A method for identifying moraine on a glacier surface based on multi-source remote sensing parameters

Abstract

The present application relates to the technical field of glacier information processing, and more particularly to a glacier 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 physically consistent data; using a cross-modal visual conversion network and a symmetric neural operator to infer the ground surface temperature under the driving of shortwave albedo, downward radiation and turbulent flux, and obtaining the thickness-presence posterior through normalization flow; constructing an uncertainty field according to information entropy and variance, selecting representative pixel unmanned aerial radar samples by means of quantum annealing, fine-tuning the model through value backflow, and then diffusing the anisotropic graph to output the spatial continuous thickness field. The present application integrates daily weather forecasts to integrate future thickness according to the law of conservation of mass, and couples with the existence index to form a breaching risk grid, thereby realizing high-precision and iterative glacial moraine dam warning.

Description

Technical Field

[0001] This invention relates to the field of glacier information processing technology, and in particular to a method for identifying glacier moraines based on multi-source remote sensing parameters. Background Technology

[0002] The front of alpine glaciers is often covered by a mixture of debris and ice (moat). The thickness of the debris controls the energy absorption and local ablation rate of the ice surface, and is also a triggering factor for the seepage and collapse of glacial moraines. Accurate knowledge of its spatial distribution and multi-year evolution is crucial for disaster early warning, regional water resource assessment, and Quaternary environmental reconstruction. Due to the rugged terrain and frigid climate of glaciers, drilling or trenching methods are costly and have limited coverage, making remote sensing technology the only feasible way to conduct macroscopic monitoring of moats. However, a single sensor can often only capture one aspect of the properties of the debris layer: the visible-near-infrared band reflects color differences, but is easily confused with coarse-grained debris and thin snow cover; thermal infrared records temperature, but is affected by diurnal variation and shading; synthetic aperture radar is sensitive to roughness and humidity, but lacks spectral discrimination. Therefore, multi-source remote sensing fusion is considered a way to overcome the current bottlenecks. Summary of the Invention

[0003] To address the numerous problems existing in the prior art, this invention provides a method for identifying glacier moraines based on multi-source remote sensing parameters. This invention physically unifies multi-domain remote sensing images, extracts cross-modal features using a bi-branch visual transformation network, and infers instantaneous temperature by combining shortwave albedo and turbulent flux inputs to a symmetric neural operator. Then, affine coupled flow outputs the thickness and posterior existence. The uncertainty field is subjected to quantum annealing to select representative pixels for field sampling. The measured thickness is used to fine-tune the model in reverse and an anisotropic diffusion is used to generate a continuous thickness field. Finally, under daily meteorological driving, the future thickness is obtained through integration, and a risk index is constructed.

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

[0005] Simultaneously 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 to generate physically consistent data.

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

[0007] An uncertainty field is constructed based on the posterior data. High-entropy and spatially correlated pixels are selected through quantum optimization and on-site sampling is carried out. The measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator. A spatially continuous thickness field is formed through anisotropic graph diffusion regularization.

[0008] The symmetric neural operator, which uses the spatial continuous thickness field and frozen weights, fuses daily weather forecasts to calculate the net radiation of pixels. The future thickness is obtained by integrating according to the mass conservation principle up to a set number of years. Based on the presence of tidal flats and the future thickness, raster data of the collapse risk index is generated.

[0009] Preferably, the cross-modal encoder employs a visual conversion network. 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. Synthetic aperture radar remote sensing images and digital elevation data are divided into a second feature sequence after being divided by the same length. The first feature sequence and the second feature sequence are processed by multi-layer self-attention modules and then concatenated in the hidden dimension to form a joint feature vector.

[0010] Preferably, the shortwave albedo is calculated based on the normalized difference between the green band reflectance and the shortwave infrared band reflectance, the downlink radiation is interpolated to a spatial resolution consistent with the physical consistency data through meteorological forecast field interpolation, and the turbulent flux is calculated based on near-surface meteorological elements using 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 configured as symmetric structures in the frequency domain according to the real and imaginary parts.

[0012] Preferably, the differentiable probability model is composed of multiple affine coupling transformation units connected in series. Each affine coupling transformation unit uses an alternating masking method to separate the input channel, and the scaling function and displacement function are implemented through a multilayer perceptron.

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

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

[0015] Preferably, the on-site active sampling collects the pixel set sequentially through the traveling sales path constructed by the UAV. The UAV is equipped with a ground wave radar to obtain radar profile and invert the thickness according to the 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.

[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 pixel slope and pixel uncertainty. The model is iterated in an explicit difference format until the average relative difference between two adjacent thickness fields is less than a preset threshold, and then the spatial continuous thickness field is output.

[0017] Preferably, the collapse risk index is calculated based on the exponential relationship between the existence of the pixel debris layer and its thickness at future time. The thickness at 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 raster form for risk classification.

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

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

[0020] This invention achieves posterior thickness inference under energy budget constraints by using a physical driving field-differentiable probability model closed loop, overcoming the defect of empirical coefficients failing with the scene.

[0021] This invention uses quantum annealing active sampling and anisotropic graph diffusion regularization to output a spatially continuous and noise-suppressed thickness field, avoiding stripe artifacts and significantly reducing fieldwork workload. Attached Figure Description

[0022] Figure 1 This is a flowchart of the present invention;

[0023] Figure 2 This is a flowchart of the data processing in this invention;

[0024] Figure 3 This is a diagram of the fusion structure of symmetric neural operators and physical driving in this invention. Detailed Implementation

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

[0026] like Figures 1-2 As shown, a method for identifying glacier moraines based on multi-source remote sensing parameters includes the following steps:

[0027] Simultaneously 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 to 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, which differ in observation time phase, imaging mechanism, spatial resolution, and coordinate reference, to the same physical semantic space, thus laying a comparable numerical foundation for subsequent cross-modal feature encoding and coupled inference. This invention introduces a two-stage strategy of "radiative-atmospheric joint correction" and "geometric-reprojection unified processing" in this stage. Compared to traditional pipeline-style item-by-item correction methods, this approach can eliminate the superposition errors of radiation amplitude drift and spatial distortion between multi-source data in one step, significantly reduce the cross-registration errors between optical, thermal infrared, and synthetic aperture radar data, and minimize the irreversible loss of image sharpness caused by repeated resampling.

[0029] The radiometric-atmospheric joint correction stage first performs top-level radiance normalization on optical and thermal infrared remote sensing images based on orbital epochs and solar altitude angles, converting radiometer digital quantization values ​​into TOA (top atmosphere) radiance. Then, a radiative transfer model is used to remove atmospheric scattering and absorption effects from the TOA radiance, directly extrapolating it back to the surface reflectivity layer. Due to the low atmospheric pressure and thin water vapor content in high-altitude glacial regions, using empirical formulas from low-altitude areas would result in systematic overcompensation errors. This invention introduces an extrapolation module for atmospheric pressure and water vapor profile parameters in the atmospheric transfer coefficient solution, and integrates microwave radiometer data to calibrate real-time water vapor column content, making the correction coefficient more suitable for the cold and thin atmosphere. For synthetic aperture radar images, this invention adopts a radiometric-geometric integrated calibration process, combining thermal noise removal, Beta0 calibration, and slant-to-ground distance correction into a single calculation, avoiding edge fringe effects caused by multi-level resampling.

[0030] In the geometric-reprojection unified processing stage, a digital elevation model generated from laser altimetry and stereo optical imagery is used as the terrain benchmark, and terrain corrections are performed on both optical and radar imagery. For optical imagery, a polynomial image equation is used to solve for attitude errors and perform adjustment, while for radar imagery, a temporal-spatial baseline-elevation fitting method is used to compensate for terrain distortion. All data are ultimately reprojected to a unified UTM coordinate system and rasterized to a 20-meter resolution. To prevent frequency domain aliasing introduced by repeated interpolation, a cubic convolution kernel is used for optical imagery, and a sinc-Lanczos finite kernel is used for radar imagery, ensuring that high-frequency textures in the Cartesian grid still have sufficient energy. After this stage, all observations are organized into row-oriented pixel sequences and stored in X. phys The matrix features are arranged strictly according to the band number, radiance, backscattering coefficient and topographic indicators, which facilitates direct slice reading in the cross-modal encoder.

[0031] In terms of effectiveness, the combined radiation-atmosphere correction controls the differences in visible and near-infrared reflectance from multiple temporal observations on the same day to within 5%, and the thermal infrared brightness temperature drift to within one Kelvin. After unified geometric-reprojection processing, the average deviation between pixel misalignment caused by SAR image distortion in mountain radar and optical pixel boundaries is reduced to within one pixel. Taking the glacier group on the eastern slope of Gongga Mountain as the experimental area, after using the preprocessing of this invention, the mean square error between the synchronously acquired Sentinel-2MSI green reflectance and the ground-measured reflectance of the airborne spectrometer is less than half that of the traditional six-segment correction process. The SAR backscattering coefficient is highly consistent with the measured incident angle scattering response curve based on the ground angle reflector column, providing a stable baseline for subsequent estimation of debris roughness in the frequency domain.

[0032] In this embodiment, the research team selected a typical glacier front region located on the Qinghai-Tibet Plateau. First, they deployed angular reflective columns and water surface diffuse reflectors in the field to evaluate the effect of geometric and radiometric dual correction. Subsequently, they acquired Sentinel-2MSI, LandsatOLI, Sentinel-1IWSLC, and ICESat-2ATL06 data with time intervals not exceeding six hours. After applying the joint correction process of this invention, optical and thermal infrared bands were uniformly converted to surface reflectance, and SAR images underwent slant range-to-ground distance correction and calibration to σ0. Finally, X... phys The matrix is ​​input into the cross-modal encoder, which exhibits significant diagonal mutual activation of spectral-radar information in the self-attention weight heatmap of the first layer, proving that high consistency preprocessing effectively improves the coupling degree of different modes in the feature space.

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

[0034] After undergoing radiation-atmospheric joint correction and geometric reprojection, the physically consistent data are organized into a matrix X. phys The matrix column features are arranged in the following order: green reflectance, near-infrared reflectance, short-wave infrared reflectance, thermal infrared brightness temperature, dual-polarization backscattering 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 with a 20-meter resolution raster. The core computational chain then proceeds to cross-modal coding, physical coupling, and probabilistic inference.

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

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

[0037]

[0038] R G R represents the surface reflectance in the green light band. SWIR This represents the surface reflectance in the shortwave infrared band. S↓ and L↓ are obtained by interpolating the shortwave beamtop radiation and longwave downwave radiation from the meteorological model. H and LE are calculated from surface meteorological elements based on the friction velocity method. All components of the physical driving field are spatially related to X. phys Alignment is performed using the nearest neighbor time.

[0039] The symmetric neural operator is responsible for converting (z,D) phys The input is mapped to the instantaneous surface temperature field. The operator network consists of four layers of cascaded Fourier integral operators. Each layer first performs a Fast Fourier Transform to map the input to the frequency domain; then, a convolution kernel κ(f) is applied in the frequency domain. The convolution kernel satisfies the symmetric constraint κ(f) = κ(-f), ensuring that the time-domain operator matrix obtained after the inverse transform is symmetric 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 facing the band, and ε is the Stefan constant. Through calculation... The L2 norm residuals ensure that the network output is consistent with the thermal diffusion-radiation equilibrium coupling equation.

[0042] The differentiable probabilistic model employs a multi-level affine-coupled normalized flow form. The input vector is the concatenated [z, T]. The normalized flow starts from the prior distribution (log-normal-Bernoulli mixture) and maps it to the posterior distribution through an invertible affine transformation. The affine coupling block uses alternating masks to divide the channels into two groups. The first group remains unchanged, while the second group undergoes an affine transformation based on the scaling function s(·) and displacement function t(·) mapped from the first group via the perceptron, thus ensuring that the Jacobian determinant is easily computed. The model maximizes the lower bound of evidence and jointly optimizes the contrastive loss, physical residual loss, and variational loss. The final output is the posterior mean μ. d ,variance With the probability of existence p c .

[0043] The motivation for introducing posterior uncertainty in this step is to provide a quantitative basis for active sampling selection. Information entropy p c log p c With variance The pixel uncertainty field is synthesized. The spatial correlation matrix, generated by the spectrum of the convolution kernel of a symmetric neural operator, can be used to solve a quadratic disordered Boolean optimization model on a quantum annealing machine, selecting a subset of pixels that is both informative and spatially representative. This active sampling strategy significantly reduces fieldwork and substantially lowers posterior thickness field noise.

[0044] Example 1 was conducted in the extremely cold region of the Eastern Pamir Plateau. The study area used physically consistent data constructed from Sentinel-2, MODIS thermal infrared, Sentinel-1, ALOS-2, and ICESat-2 as input. Under a cross-validation framework, visualization of the self-attention weights revealed that, compared to ordinary parallel convolutional networks lacking physical field constraints, the weight matrix output by the symmetric neural operator exhibited smoother radial decay in both low-frequency and high-frequency regions, indicating that the model had intrinsically learned the stationary characteristics of the debris layer-temperature field. After testing at 500 GPS-GPR co-points, the root mean square error of the thickness decreased from three meters using traditional empirical regression to 1.6 meters.

[0045] Example 2 focuses on the southern slope glacier front, where seasonal freeze-thaw cycles are significant. Remote sensing data from April and August were input into the present invention's data link to infer the thickness difference, which was then compared with the results of pole-mounted ablation monitoring. The results show that the correlation coefficient between the model's output thickness change trend and the measured ablation amount reaches 0.89, verifying that the symmetric neural operator can still respond accurately under conditions of drastic seasonal radiation flux changes. A comparison of the same experiment without the addition of a physical driving field shows that the correlation coefficient decreases to 0.65, confirming the effect of the physical driving field on improving the model's generalization ability.

[0046] Preferably, the cross-modal encoder employs a visual conversion network. 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. Synthetic aperture radar remote sensing images and digital elevation data are divided into a second feature sequence after being divided by the same length. The first feature sequence and the second feature sequence are processed by multi-layer self-attention modules and then concatenated in the hidden dimension to form a joint feature vector.

[0047] The purpose of the cross-modal encoder is to fuse the spectral-thermal radiation information carried by optical and thermal infrared remote sensing images with the structural-roughness-topographic information carried by synthetic aperture radar remote sensing images and digital elevation data within a unified feature space. The original resolution and imaging mechanism of remote sensing images differ; direct stitching can easily lead to information occlusion due to the updating of the dominant gradient of one modality. This invention employs a visual transformation network as the cross-modal encoder, achieving scale normalization and semantic alignment of the two modalities through unified patch segmentation and a self-attention mechanism.

[0048] First, fixed-length two-dimensional non-overlapping patches are established within the optical and thermal infrared channels. The patch side length is taken as an integer multiple of a 20-meter grid cell to ensure that the patch coverage area remains consistent with the subsequent self-attention window. The original pixels are expanded according to the row master sequence and filled into the matrix, with absolute positional encoding inserted along the patch dimension. Each patch is then linearly transformed and mapped to dimension d. model The vectors are used to form the first feature sequence. The synthetic aperture radar backscattering coefficients and digital elevation data undergo the same partitioning and mapping process to form the second feature sequence. The patch size at the input end of both sequences is consistent with the linear mapping matrix, thus eliminating the scale mismatch problem caused by resolution differences at its source.

[0049] The sequence enters a multi-layered self-attention module. Each layer of the self-attention module consists of a multi-head query-key-value mapping. Query vector With key vector After performing the dot product operation, the scaling factor is applied. The weight matrix W is obtained through Softmax, and the weight matrix and value vector are... Multiplication yields the context vector C. The multi-head outputs are concatenated along the channel dimension and then linearly mapped back to the original dimension, with residuals added. Compared to the structure of convolutional kernels with fixed receptive fields, the self-attention module can automatically allocate long-range and short-range dependencies according to the scene during training, allowing spectral-thermal and microwave-terrain information to interact at different scales.

[0050] This invention introduces contrast loss during the cross-modal alignment stage. The vector output of the optical-thermal infrared branch is denoted as z. opt The vector output by the radar-terrain branch is denoted as z. sar When the batch size is B, the following definition applies:

[0051]

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

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

[0054] An additional advantage of the visual transformation network is that the self-attention weight matrix provides explicit interpretability. In experiments at the glacier front, the first layer of attention formed a long-range coupling between the backscattering coefficient of the slope area and the green reflectance of the shaded area above, reflecting the intrinsic relationship between debris layer deposition and slope aspect-light conditions. This interpretability was used in the active sampling phase to construct the space-frequency correlation matrix, allowing the quantum annealing process to incorporate prior physical information rather than relying solely on statistically independent entropy values ​​for selection.

[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 ground truth map as an evaluation set, using a traditional convolutional neural network as a control. After using the cross-modal encoder of this invention, the root mean square error of thickness on the evaluation set decreased from 3.2 meters to 1.9 meters; the active sampling points selected after information entropy-variance uncertainty mapping reduced the fieldwork workload by 60% compared to random sampling, while reducing the proportion of noise patch area in the final spatial continuous thickness field from 12% to 4%. Example 2 was conducted on the north slope of the Altun Mountains, where the experimental area had cloud shadows and deep snow cover. This invention performs global self-attention alignment of the spectral features and radar patch roughness features at the cloud shadow locations, maintaining 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 the measured reflectivity-roughness test field, the optical-thermal branch focuses on patches corresponding to high-value areas of long-wavelength radiation in the debris layer, while the radar-topography branch focuses on areas of maximum directional backscattering. The fusion of these two branches in high-dimensional space significantly improves the ability to discriminate between coarse-grained and fine-grained mixed scenes, providing sufficient conditions for the thickness-presence joint distribution of the subsequent probabilistic model. The appearance 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 zone as the global context, consistent with the actual physical law of glacial moraine distribution along topographic zonation.

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

[0058] Shortwave albedo, downflow radiation, and turbulent flux together 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 uneconomical nature of point-by-point field measurements while ensuring a precise correspondence with high spatial resolution remote sensing grids. The following explanation covers the principles, implementation details, and effects.

[0059] Shortwave albedo represents the instantaneous scattering ratio of shortwave solar radiation by the Earth's surface. Glacial debris surfaces exhibit strong absorption of visible light but show even higher scattering of shortwave infrared radiation. This invention uses the green light band reflectance R... G With shortwave infrared band reflectivity R SWIR The normalized difference is used to establish an approximately parameter-free albedo estimate:

[0060]

[0061] Where α represents shortwave albedo, which is dimensionless; R G With R SWIR These are the surface reflectances of green light and shortwave infrared after radiation-atmospheric joint correction. The formula derivation is based on the empirical relationship that the difference in scattering-absorption between the two bands is inversely proportional to the overall scattering ratio under the law of energy conservation. Compared with polynomial models that require multi-band regression, this formula only relies on two high signal-to-noise ratio bands, and can maintain stable output even under cloud cover or partial sensor failure. In the embodiment, the α calculated by the formula from the green light and shortwave infrared reflectance of the Sentinel-2 multispectral instrument is compared with the integrated sphere observations of the ground spectrometer, and the average absolute error is less than 0.03.

[0062] Downward radiation consists of a shortwave component (S↓) and a longwave component (L↓), originating from solar radiation and atmospheric thermal radiation. This invention uses numerical weather prediction products with a resolution better than 3,000 meters, projecting the forecast field onto a 20-meter grid using bilinear interpolation. To avoid overestimation of shortwave radiation due to mountain shadows, the algorithm calculates a solar incidence angle correction coefficient for each pixel after interpolation based on slope aspect-solar azimuth angle; if the corrected incidence angle is less than the horizon height, the shortwave radiation is automatically set to zero. In the Gongga Mountain south slope scenario of Example 1 in the specification, this operation reduces the estimation error of S↓ in the shadowed 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 only needs to be consistent with the interpolated value. After interpolation, the row and column coordinates of the radiation field correspond perfectly with the physical consistency data, ensuring that the symmetric neural operator can read all physical quantities at once.

[0063] Turbulent fluxes refer to sensible heat flux H and latent heat flux LE, characterizing the sensible heat and water vapor exchange between air and debris-ice surface. This invention uses the frictional velocity method for calculation. First, the near-surface temperature T is extracted from meteorological forecast products. a , specific moisture q a The surface roughness length z0 is estimated using the horizontal wind speed U and the backscattering coefficient of dual-polarized synthetic aperture radar: the backscattering is stronger in the debris rough region, corresponding to a larger z0; the backscattering is weaker in the bare ice region, corresponding to a smaller z0. Then, the friction velocity u is solved using the classical Prandtl-von Kármán relation. * The expressions for sensible heat flux and latent heat flux are as follows:

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

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

[0066] Where ρ is the air density, and c p For the specific heat capacity at constant pressure, L v For the latent heat of vaporization, T s For the surface temperature predicted by the symmetric neural operator, q s The specific humidity is saturated. Flux and temperature field form an iterative coupling: the initial iteration uses remote sensing brightness temperature to estimate T. s After obtaining the flux, input it into the neural operator to update T. s This is then fed back into the flux formula until convergence. The example shows that after two iterations, T... s The convergence error is less than 0.2℃, which meets the accuracy requirements for 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 debris-ice surface. The symmetric neural operator consists of multiple layers of Fourier integral operators, with the frequency domain convolution kernel constrained in a mirror-symmetric manner, structurally ensuring a positive definite thermal diffusion matrix. The combination of self-attention and symmetric neural operators captures multimodal spatial interactions while enforcing local energy conservation. The training loss function includes three terms: contrastive loss, energy balance residual, and variational lower bound; the weights of each term are determined through cross-validation, ensuring the model maintains physical consistency while also considering generalization.

[0068] The root mean square error (RMSE) of thickness estimation for the experimental glacier group in the Pamir Plateau, as demonstrated by this invention, is reduced by more than 30% compared to traditional empirical formulas. In areas with significant debris coverage, changes in debris thickness are sensitive to changes in albedo; the addition of shortwave albedo automatically amplifies this sensitivity, resulting in a greater relative reduction in thickness error. Furthermore, the coupled turbulent flux model of this invention incorporates radar estimation of roughness length, enabling dynamic adjustment of heat exchange at the ice-debris interface. Under sudden wind speed increases, the correlation between the peak flux output of the model and surface eddy covariance observation errors is controlled within 15%, while the error of the model without radar roughness information approaches 40%. The results of Example 2 show that the energy budget closure at three-hour resolution for a tributary glacier is improved 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 configured as symmetric structures in the frequency domain according to the real and imaginary parts.

[0070] The symmetric neural operator undertakes the core responsibility of mapping "physical driving field → surface temperature" in this invention. Its design principle is to ensure that the network structure itself embodies the symmetric positive definite characteristics of the thermal diffusion-radiation coupling equation, thereby automatically satisfying energy conservation constraints during end-to-end training without relying on external hard projection or post-processing. This operator uses multi-layer Fourier integral operators as basic units. Each layer first performs convolution operations in the frequency domain, and then returns to the spatiotemporal domain through a fast inverse transform. The convolution kernel is set with its real and imaginary parts mirrored, ensuring that the time-domain kernel obtained after the inverse transform has even function properties, thus maintaining the symmetry of the overall convolution matrix in the discrete sense. The following content elaborates on the principle, implementation, and effects in sequence.

[0071] The debris-ice-air system can be approximated at a fine scale as a three-dimensional unsteady thermal diffusion-radiation coupling process. If we neglect horizontal flow and simplify the process by assuming a linear distribution of vertical temperature, we obtain a one-dimensional diffusion-source equation:

[0072]

[0073] Where T represents the instantaneous surface temperature (in K), and k represents the effective thermal diffusivity (in m³). 2 s-1 ), ρ represents the average density of the clastic layer (unit: kg / m³). -3 ), c represents specific heat capacity (unit: J / kg) -1 K -1 ), Q net Represents the net radiation-turbulent flux source term (unit: Wm). -2 If we perform a Fourier transform on the equation in the spatial direction, we get:

[0074]

[0075] This indicates that the diffusion operator exhibits -κ|f| in the frequency domain. 2 The multiplier form of the convolution kernel is proportional to the square of the frequency and is an even function. Based on this, the present invention introduces a Fourier integral operator, directly defining the convolution kernel in the frequency domain, such that the magnitude and sign of the convolution kernel with respect to f and -f satisfy:

[0076]

[0077] Where κ(f) is the frequency domain convolution kernel, The imaginary part is represented by . This mirror constraint ensures that the time-domain kernel k(x) obtained by the inverse transform simultaneously satisfies k(x) = k(-x) (an even function) and ∫k(x)dx = 0, thus 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 allows the neural operator to naturally satisfy energy conservation when predicting temperature gradients, avoiding physically meaningless negative diffusion.

[0078] In practical applications, the symmetric neural operator is constructed by cascading four layers of Fourier integral operators. Let the input be a vector. The operation process for the l-th layer (l = 1-4) is as follows: 1. Fast Fourier Transform: Frequency domain convolution: Where ⊙ represents the inverse transformation of Hadamard multiplication 3: 4. Nonlinear mapping: u (l) =σ(v (l) +b (l) ), where σ is the linear rectified function. The convolution kernel κ of each layer... (l) After initial random sampling during training, symmetry is enforced according to the aforementioned mirroring rules; during updates, only amplitude is allowed to be adjusted synchronously on both sides and phase is adjusted in opposite directions on both sides, ensuring that symmetry is not destroyed during training.

[0079] To balance deep learningability and numerical stability, this invention incorporates batch normalization into the output of each layer, ensuring that the mean of the activation values ​​is zero and the variance is one, and also applies κ... (l)An amplitude upper limit constraint is introduced to prevent high-frequency channels from amplifying indefinitely. Unlike ordinary convolutional networks, this structure can capture long-range dependencies simultaneously with frequency-domain-preferred convolution in O(n log n) complexity; symmetry 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 features output by the cross-modal encoder, and the other is a physical driving vector consisting of shortwave albedo, downlink radiation, and turbulent flux stacks. Both are mapped to the same length as the input mesh through a fully connected layer and then summed element-wise as u. (0) The output is the temperature field T. pred and with the energy balance source term Q net Calculate the physical residual:

[0081]

[0082] Neural operators with global loss:

[0083]

[0084] Joint optimization and These are derived from the cross-modal contrastive objective and the variational objective of the differentiable probabilistic model, respectively. The bidirectional gradient flow requires the symmetric neural operator to satisfy both energy conservation and provide 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 non-physical peaks. Experiments show that compared with a regular Fourier layer without symmetry constraints, the number of pixels exhibiting negative diffusion (gradient sign reversal) in the temperature gradient field is reduced by 90%, and the number of training convergence steps is shortened by about one-third. In mountainous scenarios with highly variable debris layer thickness, the model can adaptively enhance the high-frequency kernel amplitude to reduce the impact of thin, high-albedo debris layers on excessive temperature rise; while in river valleys with bare ice, the low-frequency kernel amplitude is enhanced to achieve overall smoothness of the temperature field.

[0086] Example: The front edge of the Karakoram main ridge glacier was selected, and the complete link of this invention was run based on a 20-meter grid, compared with the traditional convolutional-long short-term memory coupled model. Using 360 independent temperature-thickness-radiative flux samples as the validation set, the root mean square error of temperature for the symmetric neural operator scheme was 0.94 Kelvin, while that for the traditional scheme was more than 1 Kelvin; the coefficient of determination of the thickness posterior to the measured thickness was improved to 0.87. Due to the symmetry of the convolution kernel, the power spectral density showed that the thermal diffusion mode was symmetric and stable, and the spectral leakage phenomenon caused by high-order convolution kernels in the traditional model was not observed.

[0087] Preferably, the differentiable probability model is composed of multiple affine coupling transformation units connected in series. Each affine coupling transformation unit uses an alternating masking method to separate the input channel, and the scaling function and displacement function are implemented through a multilayer perceptron.

[0088] In the scenario of glacial moraine identification, thickness and existence are continuous-discrete mixed random variables, requiring both pixel-level point estimation and a reliable uncertainty measure. This invention uses a differentiable probabilistic model to achieve this goal. The core idea is to use an invertible transformation to map a simple prior that is easy to sample and calculate density to a complex posterior distribution; the transformation link is continuously differentiable, facilitating backpropagation with the preceding symmetric neural operator.

[0089] The model employs a normalized flow architecture composed of multiple affine coupling transformation units connected in series. Let y be the target random vector, where element one corresponds to pixel thickness and element two corresponds to the continuous representation of the existence degree of the surface after logarithmic probability transformation. Define the prior z. (0) It satisfies a log-normal-Bernoulli mixture distribution with independent dimensions. The k-th order affine coupling transformation is denoted as:

[0090]

[0091] m is a fixed binary mask, ⊙ represents element-wise multiplication, s k (·) and t k (·) represent the scaling function and the displacement function, respectively. Both are implemented using a multilayer perceptron. The perceptron input is the value of the masked channel, and the output is a vector of the same length as the transformed channel. The activation function uses a continuously differentiable rectified linear unit to ensure that the overall transformed Jacobian determinant is only related to s. k Related but not containing t k This facilitates analytical computation. The alternating masking method swaps m and 1-m in adjacent cells, ensuring that all channels are used as conditional branches at least once throughout the flow. A reversible mapping is obtained after concatenating L stages. The posterior logarithmic density is calculated by the following formula:

[0092]

[0093] in To transform the channel index, p0 is the prior density, s k,i This represents the scale output of the k-th level and i-th channel. Since the determinant is the product of the exponents of the diagonal elements, the complexity is linearly related to the number of channels and does not depend on the spatial dimension, making it suitable for batch inference of large-format remote sensing raster.

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

[0095]

[0096] Where x represents cross-modal observations, θ represents the generator network parameters, and φ represents the flow model parameters. The first term uses the conditional likelihood output by the symmetric neural operator, and the second term is the posterior-prior Kolb-Leibler divergence. Because f is invertible, the KL term can be quickly evaluated under closed-form analysis, ensuring training efficiency.

[0097] Model outputs posterior mean thickness μ h With variance and the probability of existence p c By using a closed-form gradient, the measured thickness samples obtained through active field sampling can be directly fine-tuned for φ, followed by a priori shrinkage. Due to the reversible nature of flow, thickness-existence pairs can be generated simply by performing backward reasoning on the priors during sampling, thus achieving uncertainty propagation.

[0098] The example uses the western Karakoram glacier group as the test field. A 20-meter resolution grid of 300,000 pixels was used as the training set, and 10% of the pixels were randomly selected as the test set. The ground truth thickness was obtained from a ground trench-radar joint profile. The experimental group applied a differentiable probabilistic model; the control group used a multilayer perceptron of the same size to directly regress the thickness. Results showed that the experimental group reduced the mean prediction error by 20%, and the linear correlation coefficient between the prediction variance and the absolute error reached 0.8, indicating that the variance effectively characterizes the error; the control group's variance-error correlation coefficient was less than 0.3. Further, in the active sampling phase, 1,000 pixels were selected for field measurement and feedback fine-tuning. The root mean square error of the thickness in the experimental group decreased by another 10%, while the control group's decreased by less than 5 percentage points, verifying the advantages of the flow model in incremental learning scenarios.

[0099] Another practical benefit of the differentiable probability model is that it provides spatially relevant uncertainties for subsequent risk mapping. The uncertainty field formed by the weighted combination of thickness variance and existence probability plays a core role in constructing the objective function of quantum annealing optimization, ensuring that sampling points are concentrated in critical areas with high thickness variance and existence probability close to 0.5, thus significantly improving fieldwork efficiency.

[0100] An uncertainty field is constructed based on the posterior data. High-entropy and spatially correlated pixels are selected through quantum optimization and on-site sampling is carried out. The measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator. A spatially continuous thickness field is formed through anisotropic graph diffusion regularization.

[0101] The thickness expectation, variance, and existence probability given by multi-source remote sensing inversion are pixel-level random quantities, and their reliability is affected by the imbalance of training samples, radiation noise, and model extrapolation errors. To minimize this uncertainty with limited fieldwork costs, this invention proposes a three-step closed loop: "posterior-driven - quantum optimization - anisotropic diffusion." First, the posterior statistics are converted into a measurable uncertainty field. Then, a quantum annealing optimization model is used to select the most information-rich and spatially representative sampling pixels across the entire scale. Finally, the sampling results are fed back to a differentiable probability model and a symmetric neural operator, and the thickness field is smoothed using an anisotropic diffusion method, suppressing local noise while preserving details. This three-step closed loop ensures that sampling-update-regularization is completed in one step, avoiding multiple rounds of fieldwork.

[0102] The uncertainty field is constructed such that the posterior distribution of pixel thickness is output by a differentiable probability model, whose parameter is the mean thickness μ. h ,variance And the probability of existence p c The uncertainty index U needs to reflect both the magnitude of the variance and the characteristics of the existence probability being far from zero or having the highest information content at a given moment. This invention sets:

[0103]

[0104] The first term is the information entropy, and the second term is the variance. The dimensions of the two are directly added together after normalization. c Represents the probability that a pixel is a table, and is dimensionless; This represents the thickness variance, in square meters. From this, the gridded uncertainty field {U} is obtained. i}, where i is the pixel index. Experiments revealed that using only variance ignores the region where the model has the least confidence at the classification boundary when the probability is close to 0.5; using only information entropy cannot distinguish the differences in thickness fluctuation amplitude. Combining both can accommodate both continuous and discrete uncertainties.

[0105] Traditional greedy or heuristic methods for selecting sampling points in quantum annealing often only consider the uncertainty magnitude and ignore spatial redundancy. This invention incorporates spatial correlation into a Boolean quadratic disordered optimization model, solving it in one step on a quantum annealing machine. The formula is as follows:

[0106]

[0107] Among them, s i ∈{0,1} is a binary variable indicating whether a cell is selected; a i =U i For the coefficients of the first-order term, pixels with high uncertainty tend to be selected; b ij =-λC ij λ is the positive weight, C ijThe amplitude of the spectrum of the symmetric neural operator convolution kernel at pixel pair (i,j) is larger, indicating stronger temperature coupling. It is desirable that the two are not 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 to classical hill climbing or simulated annealing, the quantum tunneling mechanism can skip local minima in an exponential combinatorial space with a shorter annealing time. The target set S = {i | s} is obtained. i After setting 1, the shortest flight path for the drone is generated and it flies over the pixel center.

[0109] On-site sampling and model fine-tuning were performed using a UAV equipped with a surface-wave radar and a dual-band thermal infrared sensor. The radar profile was inverted using time-depth transformation to retrieve the measured thickness h. * Thermal infrared brightness temperature is used to verify radiation correction. Immediately after sampling, {h} * As a new supervised pair, a differentiable probabilistic model is fine-tuned with a small learning rate and a symmetric neural operator: the symmetric neural operator receives h * Update the top temperature boundary conditions; for normalized flow, use h. * Recalculate the scale-displacement perceptron weights. Thanks to the end-to-end differentiability of the model, fine-tuning does not disrupt the existing symmetric structure of the convolutional kernels and converges in only a finite number of iterations.

[0110] Anisotropic graph diffusion regularization, even after fine-tuning, still results in a new thickness field containing random noise, requiring spatial regularization. Glacier moraine thickness exhibits high smoothness along the slope and drainage channel directions, but a large gradient along the direction perpendicular to the ice streamlines. This invention constructs an undirected graph G = (V, E), where node V corresponds to a pixel, and edge weights are:

[0111]

[0112] Where H i r represents the node elevation. ij For pixel spacing. Diffusion tensor:

[0113]

[0114] The symbol η represents the uncertainty weight, which increases the diffusion coefficient in high-uncertainty directions. Explicit finite difference is used.

[0115]

[0116] Iterate until the average absolute difference between two adjacent steps is less than a threshold. (Symbol) Let λ represent the thickness at the t-th iteration, and λ be the total variational weight used to suppress staircase artifacts. This diffusion-total variational coupling preserves abrupt changes in cross-sections in high-slope regions and smooths noise in low-slope regions, ultimately outputting a spatially continuous thickness field.

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

[0118] In the Golmud River source area of ​​Qinghai, the experimental team compared the results with classical isotropic Gaussian filtering smoothing. The results showed that Gaussian filtering caused excessive blurring of the ice tongue boundary, with a boundary error exceeding five pixels. The anisotropic diffusion of the present invention prevented the thickness value from spreading across the slope direction, and the boundary error was controlled within one pixel, verifying the effectiveness of the slope-uncertainty coupled tensor design.

[0119] By combining uncertainty measurement, quantum annealing optimization, and anisotropic diffusion, this invention achieves posterior contraction and thickness field smoothing under the premise of limited field sampling costs. It makes full use of high uncertainty information and suppresses spatial redundancy. Furthermore, end-to-end fine-tuning ensures the closure of the temperature-thickness-energy balance link, significantly improving the reliability and spatial continuity of glacier moraines identification results.

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

[0121] The output of the meta-level glacial moraine identification is a continuous-discrete mixed posterior: the thickness h follows a conditional probability density p(h|x), and the existence degree c follows a Bernoulli distribution Ber(p) c To prioritize the survey of locations with the scarcest information within limited fieldwork costs, it is necessary to map this binary posterior into a single spatially measurable field. The uncertainty field proposed in this invention is precisely such a measure, unifying the uncertainty of discrete classification and the uncertainty of continuous regression to the same dimension, and indicating the steepest direction of "missing information" in the sense of gradient.

[0122] First, calculate the existence information entropy. For a 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 cThis represents the dimensionless probability that a pixel belongs to a table category. The function has a probability of p... c The maximum value is reached at p = 0.5, reflecting the location where the model lacks the most classification confidence; when p c When the entropy approaches 0 or 1, it decays to zero, and the corresponding classification has stabilized.

[0125] Secondly, the posterior variance of thickness can be directly read from the differentiable probability model. It measures the degree of dispersion in thickness estimation given observations. When pixel terrain is steep or debris spectra are mixed, An increase in the value usually indicates that the model is still uncertain in terms of regression significance.

[0126] Since the entropy and variance terms have different dimensions, directly adding them will result in one term dominating. To align the scales, this invention first performs a minimum-maximum normalization on the variance:

[0127]

[0128] in and This represents the extreme value of the entire image variance. After normalization... The range of values ​​is consistent with that of the entropy term.

[0129] The uncertainty field U is given by the following equation:

[0130]

[0131] ω H With ω σ For the preset weights, satisfy ω H +ω σ =1. This invention defaults to equal weights. When the test set shows a classification error significantly greater than the thickness error, ω can be appropriately increased through cross-validation. H Conversely, increasing ω σ .

[0132] This linear superposition possesses two important properties. First, Both the entropy-dominated region and the variance-dominated region point to the pixel 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 sorting is stable.

[0133] At the implementation level, the computational process and posterior inference share GPU tensors: for each batch of pixels, normalized stream decoding yields μ. h , p c U is then calculated using tensor addition. The entire process involves only logarithmic, exponential, and normalization operations, preserving differentiability. When interacting with the quantum annealing platform, U simply needs to be projected as the coefficients of the first-order term {a}. i}, and simultaneously utilize the convolution spectrum to generate quadratic coupling coefficients {bij}, thus a quadratic unordered Boolean optimization model can be constructed.

[0134] Example 1: After generating a U-field in the debris-rich area of ​​the West Kunlun Mountains, three sampling schemes were compared using the same fieldwork time budget: random sampling, variance-only sampling, and entropy-variance fusion. The root mean square error of the average thickness was 2.3 m for the random scheme; it decreased to 1.9 m for the variance-only scheme; and it further decreased to 1.6 m for the entropy-variance fusion scheme, while improving the classification F1 value by 0.07. This indicates that multi-source uncertainty fusion is beneficial for both continuous and discrete dual objectives.

[0135] Example 2: In the shadow-bare ice transition zone of Mount Gongga, fifty drilling points were extracted from both entropy-dominated and variance-dominated pixels. 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 spatially exhibited a mixed pattern of "band-patches," consistent with the slope-trough network, proving that it both inherits energy-topographic control and accurately exposes the model's blind spots.

[0136] The uncertainty field of this invention also participates in the construction of the anisotropic diffusion tensor. By setting the vertical diffusion coefficient to 1 + ηU and the horizontal coefficient to 1, diffusion is automatically enhanced at locations with steep slopes and high U, eliminating noise peaks; in low uncertainty regions, the system remains unchanged, preserving details. Spectral analysis shows that after introducing adaptive uncertainty diffusion, the system's maximum eigenvalue decreases, and the number of iteration convergence steps is reduced by approximately one-third.

[0137] In summary, the entropy-variance linear fusion provides an uncertainty measure that combines physical and statistical considerations: the entropy term captures the fuzziness of the classification boundary, while 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-update-smoothing, significantly improving the reliability, spatial continuity, and field efficiency of glacier moraine identification results.

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

[0139] Quantum annealing has recently been recognized as a hardware acceleration path for solving large-scale binary combinatorial optimization problems. This invention transforms the multi-source remote sensing-driven glacial moraine sampling task into a Quadratic Unconstrained Binary Optimization (QUBO) model, using quantum annealing hardware to output the optimal pixel set in one go, achieving a unified goal of "maximizing uncertainty" and "minimizing spatial redundancy." The following explanation covers four aspects: modeling principles, coefficient construction, quantum mapping, interpretation, and field results.

[0140] Modeling principle: Assume the raster to be evaluated has N cells. Define binary variables:

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

[0142] Where s i =1 indicates that the i-th pixel is selected into the field sampling set, s i =0 means not selected. The objective Hamiltonian is written as:

[0143]

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

[0145] Quantum annealing hardware searches for the ground state on the energy surface through quantum tunneling, returning {s} i To minimize H, if the coefficient of the linear term is positive, the energy minimization tendency is to choose a. i Larger pixels; if the second-order coupling coefficient is negative, it tends to avoid selecting highly correlated pixels at the same time, thereby reducing spatial redundancy.

[0146] Construction of coefficients, coefficient a of the linear term i Directly take the uncertainty field:

[0147] a i =U i

[0148] in

[0149]

[0150] H c,i The information entropy of the existence degree of the pixel table. This is for normalizing the thickness variance. Weight ω H With ω σ The quadratic coupling coefficient b is determined beforehand on the validation set and kept fixed.ij The following is obtained through the spectral similarity of the convolution kernels of symmetric neural operators:

[0151] b ij =-λC ij

[0152] Where λ>0 is the spatial penalty factor; C ij The calculation method for spectral similarity is as follows:

[0153]

[0154] F SNO This represents the amplitude matrix of the symmetric neural operator convolution kernel in the frequency domain, with values ​​normalized to [0,1]. A large spectral amplitude indicates strong temperature coupling; setting a negative sign can cause quantum annealing to disperse the sampling of pixels with strong coupling, so as to avoid repeatedly acquiring similar information.

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

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

[0157] Enter the weight format accepted by the quantum annealing machine. This invention employs the following strategy during embedding:

[0158] Weight compression: reduce a i With b ij The terms with the largest absolute values ​​are uniformly linearly scaled to the dynamic range allowed by the hardware.

[0159] Adaptive chain strength: 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 scheduling: Select the single annealing time t ann Empirical evidence suggests that the quantum temperature window with the least injection noise is at t ann The ground state hit rate is highest when the hardware limit is reached.

[0161] After the hardware completes one annealing cycle, it outputs multiple candidate solutions. The solution with the lowest energy is selected as the final sampling scheme. 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 score, retaining the top K highest s. i The corresponding pixel is sufficient.

[0162] Post-processing and closed-loop correction, high-quality measured thickness after fieldwork. Write into the training set:

[0163] Fine-tuning of differentiable probabilistic models — The input is appended to the streaming model, and backpropagation continues for 10 rounds with a small learning rate.

[0164] Symmetric neural operator correction - Update to boundary conditions, 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 real-world projects, the error has already 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 in the Pamir Plateau.

[0166] Comparison methods: random sampling, greedy maximum entropy sampling, and spatial k-means-based partitioned sampling. Metrics: Root mean square error of thickness (RMSE), classification F1 score, and field track length. 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 Partition k-means 400 60 2.0 0.80 This invention 400 58 1.7 0.85

[0169] The results show that, under the premise of the shortest track length, the proposed method reduces the thickness RMSE by 0.7 meters compared to random sampling and by 0.4 meters compared to 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 both continuous and discrete uncertainties than simple information entropy or spatial partitioning.

[0170] Example A: High-slope debris accumulation zone. The peak uncertainty is concentrated in the transition zone from the slope toe to the slope crest. Radar measurements after quantum annealing confirmed that approximately 70% of the pixels represent a mixed layer of debris and bare ice. After model update, the local H... c A decrease of 60%, A decrease of 45%.

[0171] Example B: Low-slope bare ice area. Spectral coupling matrix C ij In the low-frequency channels, the annealing results tend to be sparse and dispersed. The updated thickness error is low across multiple ice tongue peaks, while the control group, which was not treated according to this invention, still exhibits strip noise at the same number of points.

[0172] Preferably, the on-site active sampling collects the pixel set sequentially through the traveling sales path constructed by the UAV. The UAV is equipped with a ground wave radar to obtain radar profile and invert the thickness according to the 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.

[0173] The fundamental task of the active sampling stage is to introduce a small number of high-precision in-situ thickness observations without altering the macroscopic coverage advantages of remote sensing, in order to shrink the posterior distribution and correct the temperature-energy inference link. This invention places the fieldwork process within the entire differentiable pipeline: first, it relies on quantum annealing to provide a representative set of pixels; then, it uses the traveling salesman path to generate the shortest path; finally, it uses the measured thickness as a pseudo-label with a propagable gradient to write back to the differentiable probabilistic model and the symmetric neural operator. This "task-driven - path optimization - differentiable backflow" design significantly reduces the secondary errors caused by the traditional separation of fieldwork and office work.

[0174] The UAV platform utilizes a multi-rotor aircraft with a high thrust-to-payload ratio, equipped with two types of sensors: a ground-penetrating radar (GPR) and a dual-band thermal infrared sensor (TIR). The GPR acquires the two-way time delay from debris to the ice surface to the ice substrate; the TIR synchronously records surface brightness temperature, serving as the measured boundary condition for the symmetric neural operator. Flight path planning employs a graph-theoretic Traveling Salesman Problem (TSP) model, where nodes are the pixel center coordinates output by quantum annealing, and edge weights are geodetic distances. After obtaining the Hamiltonian shortest loop through integer linear programming, a flight command set for the UAV is generated. The commands include latitude and longitude, a preset altitude, and hovering time, ensuring 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 debris-ice layers. The two-way time delay Δt is automatically acquired by the radar controller from the time difference between the echoes from the first to the third interface, and output in real time by the onboard computing unit after image stabilization filtering. The thickness d is calculated using the formula:

[0176]

[0177] Where v is the velocity of the electromagnetic wave in the measured medium; ε is the effective dielectric constant of the conventional debris-ice mixture. r The speed estimation formula was obtained through pre-calibration experiments. (c is the speed of light in vacuum) This is written into the firmware. This approach avoids multiple manual parameter adjustments in the field, improving the consistency of thickness calculations. To reduce diffuse reflection mismatch at the ice-debris interface, the radar antenna attitude controller maintains a small tilt angle relative to the ground surface normal during flight, correcting for delays based on real-time attitude angles.

[0178] The TIR sensor simultaneously acquires brightness temperatures in two atmospheric window bands, and after radiometric calibration, writes them into the temperature stack in real time for boundary backfeeding of the symmetric neural operator convolutional layer. If the TIR brightness temperature deviates from the remotely sensed temperature of the same pixel by more than a threshold, the system automatically marks the pixel to enter the high-weight queue of the next round of quantum annealing, achieving closed-loop bootstrapping.

[0179] After the sampled data is returned to the office server, two types of fine-tuning are triggered. The first type: differentiable probabilistic models... For the new supervised pair, 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: The symmetric neural operator treats the measured thickness as the depth boundary, recalculates the surface temperature-energy residual, and applies a small gradient descent to the convolution kernel amplitude to make the local diffusion-radiation coupling more closely match the field conditions.

[0180] The thickness grid obtained after fine-tuning often exhibits sudden drops or rises near the sampling points. To suppress interpolation ringing, this invention introduces anisotropic map diffusion regularization. The main diagonal element of the diffusion tensor is selected as a weighted sum of aspect and uncertainty, while the lateral component retains its baseline value. Under the explicit difference scheme, the iteration step size is adaptively determined by the spectral radius, and the iteration termination condition is that the mean absolute difference of pixels is less than a given threshold. Due to the orientation dependence of the diffusion tensor, diffusion in high aspect regions is suppressed, and ice tongue-sidewall details are preserved; diffusion increases in high uncertainty regions, and noise peaks tend to be degenerate.

[0181] Example: A 28-kilometer-long UAV flight path was deployed in the high-risk section of the Nyainqêntanglha Mountains' ice-dammed landslide lake. Quantum annealing output 960 sampling pixels, and after TSP optimization, the flight path length was compressed to 23 kilometers, 18% shorter than the manual zoning scheme. A single UAV could complete the task with only two battery replacements. Radar inversion thickness was verified by trenching, with an average absolute error of 0.7 meters; after fine-tuning, the root mean square error of the thickness decreased from 2 meters at the baseline to 1.4 meters. After 30 iterations of diffusion regularization, the proportion of noise patches decreased from 10% to 3%, and the error between the ice edge curve and the digital elevation model of the high-resolution stereo image pair was less than one pixel.

[0182] Compared to the classic "fixed grid drilling-Kriging interpolation" process, this invention increases the number of verification points by an order of magnitude and reduces the error by nearly half while maintaining the same fieldwork time. Compared to using only airborne optical-radar remote sensing inversion, it avoids systematic underestimation in areas with strong debris-ice mixture distribution. The results verify the efficiency and reliability of the quantum annealing-TSP-anisotropic diffusion closed loop: the annealing model ensures that the sampling points are optimal in both information and spatial dimensions; path planning saves flight paths; radar-thermal infrared dual-sensor complementary scale output; fine-tuning and diffusion smoothing allow field information to be accurately integrated into the global posterior in the form of gradients, achieving iterative convergence of 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 pixel slope and pixel uncertainty. The model is iterated in an explicit difference format until the average relative difference between two adjacent thickness fields is less than a preset threshold, and then the spatial continuous thickness field is output.

[0184] After active sampling is completed and the measured thickness is fed back into the model, the pixel-level thickness field still exhibits three types of local defects: the first type is random noise peaks, commonly found at locations where debris roughness changes abruptly; the second type is banded oscillations, originating from systematic errors in the directionality of radar inversion paths; and the third type is step-like fractures, mostly occurring in slope breaks where ice-debris thickness changes abruptly and the model is insufficiently fitted. Directly using isotropic smoothing (such as Gaussian filtering or conventional total variation) would lead to excessive blurring of the actual ice edge and ice surface micro-topography, affecting subsequent melting evolution and risk classification. This invention proposes an anisotropic graph diffusion regularization model to perform post-processing on the thickness field, which weakens noise while preserving spatial details controlled by topography.

[0185] The principle originates 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, which is dimensionless. To achieve direction-dependent smoothing, the second term suppresses high-frequency oscillations when the gradient magnitude is less than the noise threshold. The implementation steps are as follows:

[0188] I. Graph structure construction: For each pixel v on a uniform 20-meter resolution grid. i Create a graph node, with adjacency determined by eight-neighbor connections. Edge weights:

[0189]

[0190] H i With H j These are the elevations of the two nodes, respectively; r ij This represents the Euclidean distance between pixels. This suppresses diffusion between pixels with large elevation differences or large distances.

[0191] II. Definition of the diffusion tensor for node v i Construct the diagonal tensor:

[0192]

[0193] U i U is the uncertainty field value, which is dimensionless; η is the uncertainty amplification factor. The longitudinal component is introduced into U. iTo allow high-uncertainty pixels to release noise more quickly in the longitudinal slope direction, while maintaining the baseline value of the lateral component, excessive lateral blurring can be avoided. If the slope S... i If the value exceeds the set threshold, then multiply by the scaling factor k = S. i / S max Suppress the uphill diffusion on steep slopes.

[0194] III. Explicit difference discretization, the five-point difference scheme is written as follows:

[0195]

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

[0197]

[0198] A safety factor of 0.45 is typically used based on experience.

[0199] IV. Iteration Termination Criterion: Define 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 is, stop the iteration and output d. # =d (t+1) Threshold ε stop Take 10 -3 It can remove noise while maintaining the details of the bevel.

[0202] In an example, in the glacier tongue region of the Qilian Mountains, unnormalized thickness fields exhibit striped noise along the radar track direction. After thirty steps of anisotropic diffusion, the peak noise amplitude was reduced by 90%, and the glacier tongue boundary error remained within one pixel. In contrast, processing with an isotropic Gaussian filter at half maximum width (FWHM) of three pixels also reduced noise, but the boundary blur width expanded to five pixels. The method of this invention achieves a trade-off between boundary sharpness and noise suppression.

[0203] Another experiment was conducted at the glacier front with intense debris cover in the eastern Pamir Plateau. After incorporating slope suppression, the thickness gradient of the high-slope debris deposit wall was preserved; however, when η was set to zero, pseudo-smoothing occurred along the high-slope direction, increasing the thickness shaving error by 30%. This demonstrates that the slope-uncertainty coupled diffusion tensor is effective for complex terrain.

[0204] The symmetric neural operator, which uses the spatial continuous thickness field and frozen weights, fuses daily weather forecasts to calculate the net radiation of pixels. The future thickness is obtained by integrating according to the mass conservation principle up to a set number of years. Based on the presence of tidal flats and the future thickness, raster data of the collapse risk index is generated.

[0205] After diffusion through anisotropic maps, the spatially continuous thickness field d0(i) has eliminated random noise and preserved bend details at the grid scale, and can be used as the initial condition for thermal-mass evolution. Subsequent steps utilize a symmetric neural operator with frozen weights. Daily weather drivers are mapped to pixel net radiance, and future thickness is obtained by integral calculation based on mass conservation. A breakdown risk index grid is then constructed based on the presence of moraines. This process adheres to the ice-debris energy budget equation while avoiding regional extrapolation errors of traditional empirical daily factors. Net radiance is calculated for each pixel i and daily time step t:

[0206]

[0207] Where, α i Shortwave albedo; and ε represents the shortwave and longwave radiation under meteorological conditions; ε is the longwave emissivity of the debris-ice mixture surface; σ is the Stefan constant. The instantaneous surface temperature is the output of the symmetric neural operator.

[0208] The solution is obtained from near-surface meteorological elements using the friction velocity method. The convolution kernel spectrum of the symmetric neural operator has been frozen in the previous steps and will not be updated further to ensure that the physical consistency constraints are not violated by subsequent fine-tuning.

[0209] when When melting occurs, the latent heat consumed by melting corresponds to the change in thickness:

[0210]

[0211] ρ represents the equivalent density of the debris layer, L f This represents the latent heat of fusion of ice, where Δt is the length of a day. If... This indicates net freezing, with the thickness remaining constant, because the debris-ice mixture layer loses heat primarily through radiation at low temperatures and has no significant source of thickening. Let the daily cycle index t = 1, ..., T (T corresponding 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 it can resume calculation from the breakpoint if unexpected events occur during long-period integration.

[0214] The probability p of the existence of the surface c,i Derived from the final posterior of a differentiable probability model. The risk index is constructed based on the debris-ice dam collapse mechanism: the thinner and more abundant the debris layer, the greater the likelihood of thermal channels forming and initiating a collapse. Index form:

[0215]

[0216] satisfy When p approaches 1 c,i → Approaching 0 at 0. Using a continuous index is superior to hard threshold classification; color bands can be directly rendered in GIS and overlaid with administrative boundaries to achieve hierarchical inspection.

[0217] Preferably, the collapse risk index is calculated based on the exponential relationship between the existence of the pixel debris layer and its thickness at future time. The thickness at 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 raster form for risk classification.

[0218] Glacial till dam instability often begins with a "penetration-thinning" process of the debris layer: as the surface till layer thins year by year, internal thermal conductivity increases, pore ice is exposed, and hydrothermal channels are established to the backwater surface of the dam, triggering erosion and rapid collapse. Risk assessment must therefore consider both the presence and thinness of debris. This invention, after multi-year thermal-mass integration, maps both types of information into a single-value raster—the collapse risk index R—and directly serves GIS early warning and classification.

[0219] Future moment thickness d T Using the mass conservation equation:

[0220]

[0221] Where, d t Indicates the current thickness of the pixel, in meters (m). Represents net radiative-turbulent flux, in W / m³. -2 ρ represents the equivalent density of the debris-ice mixture, in kg / m³. -3 L f The latent heat of fusion of ice is expressed in J / kg. -1 For a full year T1 step, a five-year T5 step, or a ten-year T... 10 By integrating daily, the target annual thickness d is obtained. T .

[0222] The existence degree p of the surface c , derived from the posterior of a differentiable probability model, with a numerical interval [0,1]. p c →1 indicates that the table definitely exists, p c →0 indicates that it definitely does not exist.

[0223] Exponential mapping, based on vulnerability theory, treats thin-layer debris with high prevalence as a danger zone and employs an exponential decay function:

[0224]

[0225] Where d minIt is a minimal thickness constant to prevent division by zero. Function property: when d T Decrease or p c As p increases, R rapidly approaches 1; c →0 or d T When it is very large, R approaches 0.

[0226] In practical applications, data alignment is achieved by using the same 20-meter resolution and UTM coordinate system for the thickness and existence grids, with consistent cell indices. Parallel computing is also achieved by using tensor frameworks for one-time vectorization operations, with millions of pixels completed in milliseconds on the GPU.

[0227] Threshold grading is implemented, with three default levels: Green (R < 0.3): routine inspection; Orange (0.3 ≤ R < 0.6): seasonal inspection; Red (R ≥ 0.6): on-site special team investigation. The thresholds can be fine-tuned through historical accident review.

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

[0229] The actual application effects of this invention are shown in Table 2:

[0230] Table 2

[0231] index Experience-based daily living model This invention's integral-exponential link Thickness RMSE (m) 2.1 1.4 Historical Breakthrough Point Hitting Rate 53％ 81％ Calculation duration (10 years) 2h 12min

[0232] The results show that the combination of thermal-mass integral and physical driving field significantly reduces thickness error; exponential mapping amplifies the weight of thin-layer high-presence pixels, greatly improving the hit rate; GPU batch processing enables ten-year rolling forecasts to be completed in just over ten minutes, which can meet the annual update requirements.

[0233] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for identifying glacier moraines based on multi-source remote sensing parameters, characterized in that, Includes the following steps: Simultaneously 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 to generate physically consistent data. The physical consistency data is input into a cross-modal encoder, combined with the physical driving field composed of shortwave albedo, downlink radiation and turbulent flux, and the surface temperature is inferred using a symmetric neural operator. The temperature and joint feature vector are used as conditions to input a differentiable probability model to obtain the posterior data of the till thickness and till presence. The cross-modal encoder uses a visual conversion network. 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. Synthetic aperture radar remote sensing images and digital elevation data are divided into a second feature sequence after being divided by the same length. The first feature sequence and the second feature sequence are processed by a multi-layer self-attention module and then spliced ​​together in the hidden dimension to form a joint feature vector. The symmetric neural operator is composed of multiple Fourier integral operators connected in series, and the convolution kernels of each layer are set as symmetric structures in the frequency domain according to the real and imaginary parts; the differentiable probability model is composed of multiple affine coupling transformation units connected in series, and each affine coupling transformation unit uses an alternating masking method to separate the input channel, and the scaling function and the displacement function are implemented through a multilayer perceptron. An uncertainty field is constructed based on the posterior data. High-entropy and spatially correlated pixels are selected through quantum optimization and then sampled in-situ. The measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator. A spatially continuous thickness field is formed through anisotropic graph diffusion regularization. The uncertainty field is obtained by linearly combining the information entropy of the pixel surface existence degree and the pixel thickness variance according to preset weights. Quantum optimization uses quantum annealing to solve the quadratic disordered Boolean optimization model. The pixel uncertainty is set as the coefficient of the first term, and the spectral similarity of the convolution kernel of the symmetric neural operator is set as the coupling coefficient of the quadratic term to obtain a set of pixels with high entropy and strong spatial correlation. 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 pixel slope and pixel uncertainty. It is iterated in an explicit difference format until the average relative difference between two adjacent thickness fields is less than a preset threshold, and the output is a spatially continuous thickness field. The symmetric neural operator, which uses the spatial continuous thickness field and frozen weights, fuses daily weather forecasts to calculate the net radiation of pixels. The future thickness is obtained by integrating according to the mass conservation principle up to a set number of years. Based on the presence of tidal flats and the future thickness, raster data of the collapse risk index is generated.

2. The method according to claim 1, characterized in that, Shortwave albedo is calculated based on the normalized difference between the reflectance of the green band and the reflectance of the shortwave infrared band. Downward radiation is interpolated to a spatial resolution consistent with the physical consistency data through meteorological forecast field interpolation. Turbulent flux is calculated based on near-surface meteorological elements using the friction velocity method.

3. The method according to claim 1, characterized in that, Active sampling on-site uses a traveling salesman path constructed by a drone to collect a set of pixels sequentially. The drone is equipped with a ground-wave radar to obtain radar profiles and invert 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.

4. The method according to claim 1, characterized in that, The collapse risk index is calculated based on the exponential relationship between the existence of the cell debris layer and its thickness at future time. The thickness at 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 raster form for risk classification.

Images