Optical and sar-based high-level granular multi-scale deformation monitoring method

Abstract

The application relates to the technical field of remote sensing monitoring, and discloses a high-position loose body multi-scale deformation monitoring method based on optics and SAR, and aims to solve the problem of poor continuity and accuracy of an existing method, and mainly comprises the following steps: acquiring multi-temporal optical images, SAR images and auxiliary data and performing standardized preprocessing; extracting a horizontal displacement field from the optical images through an improved frequency domain cross-correlation algorithm; acquiring a time-series deformation field from the SAR images by using optimized time-series InSAR technology; performing spatial self-adaptive grid fusion and time interpolation based on a Gaussian process on the two kinds of deformation fields, and dynamically distributing fusion weights according to deformation variable levels, terrain conditions and data quality, so as to generate a time-space continuous comprehensive deformation field; realizing automatic classification of deformation modes through principal component analysis and K-means clustering algorithm, defining activity levels by combining the conversion of pre-seismic and post-seismic deformation modes, and outputting a deformation type map and an activity classification map. The application improves the continuity and accuracy of high-position loose body monitoring.

Description

Technical Field

[0001] This invention relates to the field of remote sensing monitoring technology, specifically to a method for monitoring multi-scale deformation of high-altitude loose bodies based on optics and SAR. Background Technology

[0002] As global warming continues to intensify, the retreat of glaciers and the thawing of permafrost in high-altitude regions are accelerating, leading to the accumulation of large amounts of loose material on steep slopes, forming potential sources of disaster. These high-altitude loose masses are highly susceptible to transforming into disaster chains under external influences, posing a serious threat to downstream areas. However, traditional geological hazard investigation methods face difficulties in implementation, high costs, and limited coverage in high-altitude, complex terrain areas, making it difficult to meet the needs of large-scale continuous monitoring.

[0003] With the development of remote sensing technology, optical and synthetic aperture radar (SAR) remote sensing, due to its advantages of large-scale, periodic, and non-contact observation, has been gradually applied to geological disaster identification and monitoring. Currently, deformation monitoring of high-altitude loose bodies mainly relies on two types of technologies: one is the pixel offset tracking (POT) method based on optical images, which can extract horizontal surface displacement with pixel-level accuracy; however, this method is significantly constrained by cloud and rain conditions, and high-resolution image acquisition is costly. The other is synthetic aperture radar interferometry (InSAR) technology, especially time-series InSAR methods (such as PS-InSAR and SBAS-InSAR), which possesses millimeter-level deformation detection capabilities and all-weather observation advantages; however, signal decoherence is prone to occur in high-altitude snow and ice areas and vegetation-covered areas, and the overlapping and shadowing effects caused by steep terrain also limit the completeness of deformation information extraction. In existing research, technical solutions that are relatively close to this invention still have significant limitations. On the one hand, most methods rely on a single remote sensing data source, such as using only optical images to monitor glacier movement or relying solely on InSAR technology to capture deformation, leading to monitoring blind spots in rapidly deforming areas or complex environments. On the other hand, although existing studies have attempted to combine optical and SAR data, they mostly remain at the level of simple comparison or overlay, failing to establish an effective spatiotemporal fusion model to address the inherent differences between the two types of data in terms of spatial coverage, temporal sampling, and deformation level. How to dynamically fuse multi-source information based on deformation stage, terrain features, and data quality to generate a seamless, continuous, and highly reliable comprehensive deformation field remains an unsolved technical challenge. Furthermore, deformation feature identification is mostly based on manual interpretation or fixed threshold segmentation, lacking automated and intelligent classification methods, making it difficult to accurately capture the nonlinear behavior in the spatiotemporal evolution of deformation. In addition, existing methods are mostly post-hoc analyses, lacking the ability to proactively identify and analyze the multi-factor coupling influence mechanisms, and have not yet established reliable quantitative prediction models.

[0004] Therefore, there is an urgent need to develop a technical method that can effectively integrate multi-source remote sensing data to achieve automatic identification of deformation features, so as to improve the continuity and accuracy of loose body deformation monitoring in high-altitude areas. Summary of the Invention

[0005] This invention aims to solve the problem that existing technologies cannot achieve continuous, high-precision monitoring and automatic identification of high-altitude loose bodies throughout the entire process from creep to instability. It proposes a multi-scale deformation monitoring method for high-altitude loose bodies based on optics and SAR.

[0006] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0007] A method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR, the method comprising:

[0008] Acquire multi-source remote sensing data of the study area, including multi-temporal optical satellite imagery, multi-temporal SAR satellite imagery, and auxiliary data, including digital elevation models, meteorological observation data, and earthquake catalog data;

[0009] Radiometric calibration, atmospheric correction, and terrain correction are performed on the multi-temporal optical satellite images. Precise orbit correction, thermal noise removal, radiometric calibration, and adaptive filtering are performed on the multi-temporal SAR satellite images. All auxiliary data are uniformly resampled to a consistent spatiotemporal reference to obtain a standardized multi-source remote sensing dataset with spatiotemporal consistency.

[0010] The improved frequency domain cross-correlation algorithm is used to process multi-temporal optical satellite images in the multi-source remote sensing dataset, including image fine registration, pixel offset calculation, adaptive window adjustment and system error correction, to obtain the horizontal displacement field of the study area.

[0011] Based on the multi-temporal SAR satellite images and digital elevation models in the multi-source remote sensing dataset, the temporal deformation field of the study area is obtained through optimized interferometric pair selection, atmospheric phase delay and orbital error correction, geometric distortion identification and processing, and robust temporal deformation solution.

[0012] Based on the horizontal displacement field and temporal deformation field, spatial adaptive grid fusion, Gaussian process-based temporal interpolation, and dynamic allocation of fusion weights for optical and SAR data are performed under a unified spatiotemporal framework to obtain a spatiotemporally continuous comprehensive deformation field.

[0013] Based on the comprehensive deformation field, feature extraction and dimensionality reduction are performed through principal component analysis, and unsupervised K-means clustering algorithm is used to automatically classify deformation patterns. The activity level is defined by combining the category conversion of pre-earthquake and post-earthquake deformation patterns, and deformation type map and activity level map of each observation point in the study area are obtained.

[0014] Furthermore, when performing terrain correction on the multi-temporal optical satellite imagery, the C-correction model is adopted, and its expression is:

[0015] ;

[0016] in, This represents the radiance value after terrain correction. This represents the original radiance value. Indicates the zenith angle of the sun. Indicates the angle of incidence. This represents the correction parameters determined through statistical regression on flat terrain areas;

[0017] When performing adaptive filtering on the multi-temporal SAR satellite imagery, the size of the filtering window is dynamically calculated and adjusted based on the local coherence coefficient. The calculation formula is as follows:

[0018] ;

[0019] in, Indicates the size of the filtering window. and These represent the preset minimum and maximum window sizes, respectively. This represents the local coherence coefficient.

[0020] Furthermore, the pixel offset calculation process includes:

[0021] Two-dimensional fast Fourier transform is performed on the corresponding windows of the main image and the slave image to obtain the spectra of the main image and the slave image respectively. The normalized cross power spectrum is calculated based on the spectra of the main image and the slave image. The inverse Fourier transform is performed on the normalized cross power spectrum to obtain the spatial correlation surface. The pixel-level displacement is determined by finding the peak value of the correlation surface.

[0022] Within the neighborhood of the peak position determined by the pixel-level displacement, the discrete related surface intensity values ​​are fitted by least squares using a two-dimensional Gaussian surface model to solve for the sub-pixel peak position as the final pixel offset.

[0023] The formula for calculating the normalized cross-power spectrum is as follows:

[0024] ;

[0025] in, Represents the normalized cross-power spectrum. and Represents frequency domain coordinate variables, Represents the spectrum of the main image. Indicates the spectrum of the image. express The complex conjugate, Indicates a small quantity to prevent division by zero;

[0026] The two-dimensional Gaussian surface model is as follows:

[0027] ;

[0028] in, Indicates the relevant surface strength value. and Represents spatial domain coordinate variables, This is the sub-pixel peak position. Indicates peak amplitude. and These represent the Gaussian distributions in... and Standard deviation of direction Indicates the background noise level. This represents the natural exponential function.

[0029] Furthermore, the adaptive window adjustment process includes:

[0030] During the pixel offset calculation process, the local standard deviation of each potential window position is calculated in real time, and its deformation gradient magnitude is estimated. Based on the local standard deviation and deformation gradient magnitude, the window size for cross-correlation calculation at that position is dynamically calculated and determined. The calculation formula is as follows:

[0031] ;

[0032] in, This indicates the adaptively adjusted window size. Indicates the base window size. Indicates local standard deviation. This represents the minimum standard deviation threshold. Represents the deformation gradient mode. Indicates the maximum expected deformation. and Indicates the weighting coefficient;

[0033] The process for correcting system errors includes:

[0034] After calculating the initial displacement field, a nonlocal mean filtering algorithm is used to filter the initial displacement field. The nonlocal mean filtering algorithm is as follows:

[0035] ;

[0036] in, Indicates the filtered first... Pixel offset of each pixel , and Represented by pixels and The neighborhood block centered on, Indicates the filter strength parameter. Represents the normalization factor. Indicates the first before filtering The pixel offset of each pixel.

[0037] Furthermore, the preferred processing flow for the interference pair includes:

[0038] Based on multi-temporal SAR satellite imagery, for all possible combinations of interferometric pairs, the weight of each interferometric pair is calculated based on its spatial baseline, temporal baseline, and predicted coherence. An optimized set of interferometric pairs is then constructed based on the weights for robust solution of temporal deformation.

[0039] The weights of the interference pairs are calculated using the following formula:

[0040] ;

[0041] in, Indicates the weight of the interference pair. Indicates the coherence of the prediction. Represents the vertical spatial baseline. Indicates the critical baseline. Indicates the time baseline of the interferometer pair. This represents the decorrelation time constant of the features. Represents the natural exponential function;

[0042] The atmospheric phase delay processing procedure includes:

[0043] Based on the aforementioned digital elevation model and interferometric phase map, regression coefficients are estimated through the phase-elevation statistical relationship in stable regions. Then, an elevation correlation model is used to calculate and remove topographically related atmospheric delay signals. The elevation correlation model is as follows:

[0044] ;

[0045] in, Indicates atmospheric delay, Indicates terrain elevation. Represents the regression coefficient. Represents the residual atmospheric phase component;

[0046] The geometric distortion recognition process includes:

[0047] Based on the digital elevation model and SAR satellite orbit parameters, the elevation of each pixel is calculated. Exponential value, when When it is determined to be a shaded area, If an area is identified as an overlay area, the identified shadow and overlay areas are masked out.

[0048] The The formula for calculating the index value is as follows:

[0049] ;

[0050] in, express Index value, Indicates the slope of the terrain. Indicates slope direction. Indicates the radar incident angle. Indicates the satellite azimuth angle.

[0051] Furthermore, the processing flow for the time-series deformation robustness calculation includes:

[0052] Based on the optimized set of interferometric pairs, the observation equation is constructed as follows: ;

[0053] The deformation rate vector is solved using an iterative weighted least squares estimation method: ;

[0054] in, This represents the unwrapped phase vector. Represents the design matrix. Represents the deformation rate vector. Indicates observation noise. Represents the weight matrix. Represents the smoothing constraint matrix. Represents the regularization parameter. This indicates the matrix transpose.

[0055] Furthermore, the spatial adaptive mesh fusion processing flow includes:

[0056] The evaluation function values ​​corresponding to different grid sizes are calculated. The optimal grid size is determined by finding the optimal value of the evaluation function, thus achieving a balance between data coverage and spatial resolution. The corresponding evaluation function is as follows:

[0057] ;

[0058] in, Indicates the grid size. Indicates grid size The corresponding evaluation function value, This indicates the number of grid cells containing observation points. Indicates the total number of grid cells. Indicates the feature scale. Indicates the weighting coefficient. Represents the natural exponential function;

[0059] After determining the optimal grid size, within each grid cell, inverse-range weighted interpolation is used to fuse optical and InSAR deformations:

[0060] ;

[0061] ;

[0062] in, Represents the shape variables after fusion. Indicates the first The weight of each observation point Indicates the first The deformation of each observation point Indicates the first The observation accuracy of each observation point Indicates the first The distance from each observation point to the center of the grid. Indicates the smoothing parameter;

[0063] In the Gaussian process-based time interpolation, the deformation time evolution is expressed as:

[0064] ;

[0065] in, Indicates time The shape variable, Represents the mean function, Represents a Gaussian process. Represents the covariance function. and Represents the covariance function at different time points. For the Matérn core:

[0066] ;

[0067] in, Indicates the signal variance. Indicates the time scale of the feature;

[0068] In the dynamic allocation of fusion weights for optical and SAR data based on deformation level, terrain conditions, and data quality, the weights for optical and SAR data are calculated and determined according to the following formulas:

[0069] ;

[0070] ;

[0071] in, Indicates the weight of optical data. Indicates SAR data weights, Indicates the deformation rate. Indicates the deformation rate threshold. Indicates the slope of the terrain. The characteristic parameter representing the deformation rate, Indicates slope characteristic parameters, and These represent the quality indicators for optical and SAR data, respectively.

[0072] Furthermore, feature extraction and dimensionality reduction are performed through principal component analysis, including:

[0073] Standardize the temporal deformation data of the deformation matrix corresponding to the comprehensive deformation field:

[0074] ;

[0075] in, Indicates the first The observation point at the ... The deformation at each time sampling point Represents the deformable variable after standardization. Indicates the first The mean of the time-sampled shape variables for each observation point. Indicates the first Standard deviation of all time-sampled deformation variables at each observation point;

[0076] Construct the covariance matrix, perform eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​and eigenvectors. The covariance matrix is ​​as follows:

[0077] ;

[0078] in, Represents the covariance matrix. This represents the deformation matrix corresponding to the standardized deformation variables. Indicates matrix transpose. Indicates the number of observation points;

[0079] Projecting the deformation matrices corresponding to the standardized deformation variables onto the eigenvectors yields the principal components:

[0080] ;

[0081] in, Indicates the first Principal components, The eigenvalue obtained by eigendecomposition 1 eigenvector;

[0082] Calculate the cumulative variance contribution rate of the eigenvalues ​​and select the smallest eigenvalue that results in a cumulative variance contribution rate of 85% or higher. The value is used as the number of principal components, that is:

[0083] ;

[0084] in, The eigenvalue obtained by eigendecomposition 1 eigenvalue, Indicates the number of principal components. Indicates the total number of eigenvalues;

[0085] When using the unsupervised K-means clustering algorithm, the process for determining the optimal number of clusters includes:

[0086] Define a range for the number of candidate clusters. For each candidate, perform K-means clustering, treating the deformation of each observation point as a sample point. Calculate the silhouette coefficient of all sample points and take the average. Select the candidate that maximizes the average silhouette coefficient as the final number of clusters. The formula for calculating the silhouette coefficient is as follows:

[0087] ;

[0088] in, Indicates the first Profile coefficients for each sample point Indicates the first The average distance from each sample point to other sample points in the same cluster. Indicates the first The average distance from each sample point to all sample points in the nearest neighbor cluster.

[0089] During unsupervised K-means clustering initialization, the selection probability of each sample point is calculated, and the initial cluster centers are selected based on the probability. The formula for calculating the probability is as follows:

[0090] ;

[0091] in, Indicates the first sample points The probability of being selected as the next cluster center. Represents sample points Euclidean distance to the nearest cluster center Indicates the first sample points Euclidean distance to the nearest cluster center Indicates the number of observation points;

[0092] During the iterative update process of the unsupervised K-means clustering algorithm, the cluster centers are updated according to the weighted average of the member points:

[0093] ;

[0094] in, Indicates the first The coordinates of the center point of each cluster Indicates assignment to the first The weights for each sample point are calculated based on its measurement uncertainty and reference accuracy:

[0095] ;

[0096] in, Indicates the first Measurement uncertainty for each sample point Indicates the reference precision.

[0097] Furthermore, the activity levels include the following four categories:

[0098] Non-sensitive areas: Stable before the earthquake - Stable after the earthquake;

[0099] Self-stabilizing zone: pre-earthquake activity - post-earthquake stability;

[0100] Continuous activity zone: pre-earthquake activity - post-earthquake activity;

[0101] Earthquake activation zone: stable before earthquake - active after earthquake.

[0102] Furthermore, the method also includes:

[0103] Based on the activity grading map, meteorological observation data, and comprehensive deformation field, quantitative analysis and time-delay correlation analysis of multi-factor triggering mechanisms are performed. Critical state identification is performed based on the Fukuzono-Voight model, and early warning information is generated according to the preset graded early warning criteria.

[0104] The quantitative analysis includes: precipitation triggering effect modeling, temperature triggering effect modeling, and earthquake triggering effect assessment;

[0105] The precipitation triggering effect modeling includes: establishing a quantitative relationship between deformation acceleration and cumulative precipitation based on a sliding time window, as shown in the following model:

[0106] ;

[0107] in, express Deformation acceleration at time t, Indicates background acceleration. This represents the cumulative precipitation within the sliding time window. Indicates the time lag in precipitation response. Indicates the precipitation impact coefficient;

[0108] Modeling the temperature-triggered effect includes calculating the positive accumulated temperature index used to quantify the intensity of freeze-thaw cycles.

[0109] ;

[0110] in, Indicates positive accumulated temperature index, This represents the average daily temperature. Indicates the melting threshold. Indicates the time period for summation;

[0111] The assessment of earthquake triggering effects includes: using a modified Newmark displacement model, combined with peak ground acceleration, duration, and initial stability of high-risk areas, to assess the cumulative impact of earthquake events on the stability of loose bodies and generate earthquake impact factors;

[0112] The process of the time-delay correlation analysis includes:

[0113] Based on the comprehensive deformation field and meteorological observation data, the cross-correlation coefficient between the deformation time series and the triggering factor time series at different precipitation response lags is calculated. The optimal precipitation response lag is determined by finding the maximum cross-correlation coefficient. The formula for calculating the cross-correlation coefficient is as follows:

[0114] ;

[0115] in, Represents the cross-correlation coefficient. Indicates the deformation time sequence, Indicates the timing of triggering factors. and These represent the corresponding means;

[0116] A multi-factor comprehensive triggering index is constructed to normalize and weight the influence of different triggering factors under the optimal precipitation response time lag:

[0117] ;

[0118] in, Indicates the comprehensive trigger index, and These represent the normalized precipitation and temperature indices, respectively. Indicates the earthquake impact factor. These represent the corresponding weights;

[0119] The critical state identification process includes: based on the comprehensive deformation field, calculating the reciprocal of the deformation rate time series data of the target region within a sliding time window, and performing linear fitting using the Fukuzono-Voight model. When the slope of the fitted line is negative and the goodness of fit is greater than the goodness of fit threshold, it indicates that the deformation is accelerating and the system is tending to become unstable. The Fukuzono-Voight model is as follows:

[0120] ;

[0121] in, Indicates the deformation rate. Indicates time, and This represents the model parameters obtained by fitting deformation data within a sliding window.

[0122] The beneficial effects of this invention are as follows: The multi-scale deformation monitoring method for high-altitude loose bodies based on optics and SAR provided by this invention firstly combines the advantages of optical imagery capturing meter-level large deformations with SAR imagery detecting millimeter-level micro deformations, achieving continuous high-precision monitoring of high-altitude loose bodies throughout their entire lifecycle from initial creep to accelerated instability, significantly improving the completeness and accuracy of deformation information in the spatiotemporal dimensions. Secondly, this invention utilizes automated deformation classification technology based on machine learning, changing the inefficient and subjective mode that relies on manual interpretation, and achieving objective, quantitative, and efficient classification and identification of massive temporal deformation data. This provides a reliable and comprehensive data foundation for the subsequent construction of accurate and advanced geological disaster early warning models, comprehensively improving the prevention and control capabilities and decision-making level of high-altitude loose body disasters. Attached Figure Description

[0123] Figure 1 A flowchart illustrating the multi-scale deformation monitoring method for high-altitude loose bodies based on optics and SAR provided in this embodiment;

[0124] Figure 2 A schematic diagram showing the comparison and acceleration process of deformation rates in different regions before and after the earthquake, provided for an embodiment.

[0125] Figure 3 A schematic diagram of the grid resolution optimization curve provided for the embodiment;

[0126] Figure 4 A scatter plot of deformation rate and weight distribution provided for an embodiment;

[0127] Figure 5 This is a data coverage comparison diagram provided for an example.

[0128] Figure 6 A schematic diagram of the activity classification results provided for the embodiments;

[0129] Figure 7 A schematic diagram illustrating the temporal deformation evolution of the disaster event in October 2018, provided as an example.

[0130] Figure 8 A schematic diagram illustrating the time-series precipitation and temperature evolution of a disaster event in October 2018, provided as an example.

[0131] Figure 9 This is a schematic diagram illustrating the multi-factor triggering analysis of the disaster event in October 2018, provided as an example. Detailed Implementation

[0132] To overcome the limitations of a single remote sensing data source and achieve continuous and automated deformation monitoring and activity assessment of high-altitude loose bodies from initial creep to rapid instability, this invention proposes a technical solution. First, precise corrections are performed on the different physical characteristics and error sources of optical and SAR data, unifying them onto a consistent spatiotemporal reference. Then, using an improved frequency domain cross-correlation algorithm, large-scale (meter-level) horizontal deformation is extracted with high precision from optical image sequences. Optimized temporal interferometry techniques are used to precisely detect minute deformations along the radar line of sight from SAR image sequences. The two types of deformation fields, differing in spatial distribution, temporal sampling, and magnitude, are fused into a spatiotemporally continuous comprehensive deformation field through adaptive gridding and Gaussian process-based temporal interpolation. Next, the fused comprehensive deformation field is treated as a high-dimensional dataset, and dimensionality reduction is performed using principal component analysis to extract core deformation features. Subsequently, an unsupervised clustering algorithm is used to automatically group points with similar deformation evolution behaviors into the same category, thereby achieving objective quantitative classification of deformation patterns. Finally, by analyzing the transformation of cluster categories before and after the earthquake, high-risk areas activated by the earthquake are automatically identified, thus completing an intelligent assessment of the activity of loose bodies.

[0133] The technical solutions in this embodiment will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0134] Figure 1 A flowchart illustrating a method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR is shown. Please refer to [link / reference]. Figure 1 The method includes the following steps:

[0135] Step 1: Acquisition and preprocessing of multi-source remote sensing data:

[0136] Acquire multi-source remote sensing data of the study area, including multi-temporal optical satellite imagery, multi-temporal SAR satellite imagery, and auxiliary data, including digital elevation models, meteorological observation data, and earthquake catalog data;

[0137] Radiometric calibration, atmospheric correction, and terrain correction are performed on the multi-temporal optical satellite images. Precise orbit correction, thermal noise removal, radiometric calibration, and adaptive filtering are performed on the multi-temporal SAR satellite images. All auxiliary data are uniformly resampled to a consistent spatiotemporal reference to obtain a standardized multi-source remote sensing dataset with spatiotemporal consistency.

[0138] The core of this step lies in establishing a standardized multi-source remote sensing data acquisition and preprocessing process to form a high-quality basic dataset, providing support for subsequent deformation detection. By formulating a data selection strategy that conforms to the spatiotemporal evolution characteristics of high-altitude loose bodies and implementing systematic preprocessing, radiometric, geometric, and temporal differences between different sensors are eliminated, achieving consistent integration of multi-source remote sensing data.

[0139] For optical satellite imagery, images with a spatial resolution better than 10 meters and a revisit period of less than 16 days are preferred, including multi-platform images such as Planet (3.7 meters), RapidEye (5 meters), and Sentinel-2 (10 meters). Strict quality control is implemented during selection: cloud cover must be less than 30%, and the solar altitude angle must be greater than 30° to ensure image quality and minimize shadow interference. Timing selection is primarily based on seasonal transition periods, such as late spring / early summer and late autumn / early winter, when loose body activity is significant and deformation characteristics are easier to identify. Preprocessing encompasses three key stages: radiometric calibration, atmospheric correction, and topographic correction. Atmospheric correction employs an improved 6S radiative transfer model, enhancing correction accuracy by incorporating locally measured atmospheric parameters; topographic correction is based on the C-correction model to eliminate radiation distortion caused by topographic undulations, its expression being:

[0140] ;

[0141] in, This represents the radiance value after terrain correction. This represents the original radiance value. Indicates the zenith angle of the sun. Indicates the angle of incidence. This represents the correction parameters determined through statistical regression on flat terrain.

[0142] SAR satellite image processing utilizes the Sentinel-1 satellite's interferometric wide-swath mode product, simultaneously acquiring both ascent and descent data to optimize observation geometry. Data acquisition spans at least two years, ensuring coverage of the entire annual cycle and the periods before and after key events. Preprocessing includes precise orbit correction, thermal noise removal, and radiometric calibration. To enhance interferometric coherence in complex terrain areas, this embodiment introduces an adaptive filtering algorithm. The filtering window size is dynamically calculated and adjusted based on the local coherence coefficient, as shown in the following formula:

[0143] ;

[0144] in, Indicates the size of the filtering window. and These represent the preset minimum and maximum window sizes, respectively. This represents the local coherence coefficient.

[0145] The integration of auxiliary data is a crucial step in ensuring monitoring accuracy. The system collects various auxiliary data, including high-precision digital elevation models, meteorological observation data, and earthquake catalogs. For DEM data, TanDEM-X or ALOSPALSAR products with a resolution higher than 12 meters are preferred, supplemented by 30-meter SRTM data for missing areas. Meteorological data includes daily average temperature, cumulative precipitation, and snow depth, and station observations are interpolated to the entire area using an inverse distance weighting method. Earthquake data is selected based on events with a magnitude greater than 3.0 and a focal depth less than 30 kilometers for subsequent earthquake triggering effect analysis. All auxiliary data are uniformly resampled to a 30-meter spatial resolution and aligned with the remote sensing observations, thus establishing a spatiotemporally consistent foundation for multi-source remote sensing data fusion.

[0146] Step 2, Deformation extraction of optical path:

[0147] An improved frequency domain cross-correlation algorithm is used to process multi-temporal optical satellite images in the multi-source remote sensing dataset, including image fine registration, pixel offset calculation, adaptive window adjustment, and system error correction, to obtain the horizontal displacement field of the study area.

[0148] This step aims to achieve high-precision extraction of sub-pixel displacement fields from optical images through an improved frequency domain cross-correlation algorithm, with particular optimization for large-scale rapid deformation of high-altitude loose bodies. The traditional COSI-Corr algorithm is susceptible to registration errors and noise interference in complex mountainous terrain. This embodiment significantly improves the accuracy and reliability of deformation detection by introducing an adaptive window strategy, multi-scale correlation analysis, and a robust error correction mechanism.

[0149] Accurate image registration is a prerequisite for ensuring accurate deformation detection. This embodiment employs a hierarchical registration strategy. First, initial registration parameters are obtained through SIFT feature point matching, followed by fine-tuning in local areas. The feature point density is adaptively adjusted according to the terrain complexity: 20–30 feature points are extracted per square kilometer in flat areas, increasing to 50–80 per square kilometer in steep mountainous regions. A rational function model (RFM) is used for registration transformation, and the parameters are solved using least squares adjustment. To eliminate mismatches, the RANSAC algorithm is introduced for gross error removal, with an iteration threshold set at 1.5 times the mean square error. The final registration accuracy is evaluated based on the checkpoint residuals, requiring a root mean square error of less than 0.3 pixels.

[0150] In this embodiment, the pixel offset is calculated by using an improved frequency domain cross-correlation algorithm to calculate the local displacement within a sliding window. The pixel offset calculation process includes:

[0151] Two-dimensional fast Fourier transforms are performed on the corresponding windows of the main image and the slave image to obtain their spectra. Normalized cross power spectra are calculated based on the spectra of the main image and the slave image. Inverse Fourier transforms are performed on the normalized cross power spectra to obtain the spatial correlation surface. Pixel-level displacement is determined by finding the peak of the correlation surface. In the neighborhood of the peak position determined by the pixel-level displacement, the discrete correlation surface intensity values ​​are fitted by least squares using a two-dimensional Gaussian surface model to solve for the sub-pixel peak position as the final pixel offset.

[0152] The formula for calculating the normalized cross-power spectrum is as follows:

[0153] ;

[0154] in, Represents the normalized cross-power spectrum. and Represents frequency domain coordinate variables, Represents the spectrum of the main image. Indicates the spectrum of the image. express The complex conjugate, Indicates a small quantity to prevent division by zero;

[0155] The two-dimensional Gaussian surface model is as follows:

[0156] ;

[0157] in, The correlation surface intensity value represents the intensity value of a continuous correlation surface fitted by a two-dimensional Gaussian function within a local neighborhood centered on integer-order correlation peak points. It is a continuous function used to simulate and approximate the distribution of real discrete cross-correlation values. and Represents spatial domain coordinate variables, This is the sub-pixel peak position. Indicates peak amplitude. and These represent the Gaussian distributions in... and Standard deviation of direction Indicates the background noise level. This represents the natural exponential function.

[0158] In this embodiment, the adaptive window adjustment strategy can dynamically optimize the relevant window size based on local texture richness and deformation gradient. Texture richness is evaluated by calculating the local standard deviation, and deformation gradient is estimated by the displacement difference between adjacent windows. The adaptive window adjustment process includes:

[0159] During the pixel offset calculation, the local standard deviation of each potential window position is calculated in real time, and its deformation gradient magnitude is estimated. Based on the local standard deviation and deformation gradient magnitude, the window size for cross-correlation calculation at that position is dynamically calculated and determined. The calculation formula is as follows:

[0160] ;

[0161] in, This indicates the adaptively adjusted window size. Indicates the base window size. Indicates local standard deviation. This represents the minimum standard deviation threshold. Represents the deformation gradient mode. Indicates the maximum expected deformation. and Indicates the weighting coefficient;

[0162] This adaptive strategy automatically increases the window size in areas with poor texture to improve the signal-to-noise ratio, and decreases the window size in areas with large deformation gradients to maintain spatial resolution.

[0163] In this embodiment, a three-level correction strategy is adopted: the first level is track error correction, which uses a robust estimation method to fit the residual displacement surface of the stable region and automatically removes outliers exceeding three times the standard deviation; the second level targets strip noise, using a directional notch filter to suppress specific frequency components, with its bandwidth adaptively adjusted according to the strip period; the third level is random noise suppression, employing a nonlocal mean filtering algorithm to effectively reduce noise while maintaining clear deformation boundaries.

[0164] ;

[0165] in, Indicates the filtered first... Pixel offset of each pixel , and Represented by pixels and The neighborhood block centered on, Indicates the filter strength parameter. Represents the normalization factor. Indicates the first before filtering The pixel offset of each pixel.

[0166] Through the above multi-level correction, the systematic error of deformation detection is reduced to less than 0.1 pixels, and the random error is reduced by more than 60%.

[0167] Step 3: SAR path deformation extraction:

[0168] Based on the multi-temporal SAR satellite images and digital elevation models in the multi-source remote sensing dataset, the temporal deformation field of the study area is obtained through optimized interferometric pair selection, atmospheric phase delay and orbital error correction, geometric distortion identification and processing, and robust temporal deformation calculation.

[0169] This step focuses on detecting early creep and minute deformations in high-altitude loose bodies. Through optimized small baseline set interferometry (SBAS) techniques, it achieves millimeter-level accuracy in monitoring time-series deformations. Traditional SBAS methods are susceptible to interference from factors such as loss of correlation, atmospheric delay, and geometric distortion in complex mountainous terrain, limiting monitoring accuracy and coverage. Therefore, this embodiment systematically improves multiple aspects, including interferometric pair optimization, error correction, and geometric distortion processing, significantly enhancing the extraction capability and reliability of minute deformation signals.

[0170] The optimal selection of interferometric pairs is fundamental to the quality of time-series analysis. This embodiment proposes an intelligent screening mechanism based on spatiotemporal baselines and coherence prediction: the temporal baseline is controlled within 180 days to maintain continuous capture of seasonal deformation; the spatial baseline is constrained to within 40% of the critical baseline to maintain high interferometric coherence. The key innovation lies in the introduction of a coherence prediction model, which integrates land cover type, topographic slope, and seasonal factors to predict the expected coherence of each interferometric pair.

[0171] For N multi-temporal SAR satellite images, for all possible combinations of interferometric pairs, the weight of each interferometric pair is calculated based on its spatial baseline, temporal baseline and predicted coherence, and an optimized set of interferometric pairs is constructed according to the weights for time-series deformation robustness solution.

[0172] The weights of the interference pairs are calculated using the following formula:

[0173] ;

[0174] in, Indicates the weight of the interference pair. Indicates the coherence of the prediction. Represents the vertical spatial baseline. Indicates the critical baseline. Indicates the time baseline of the interferometer pair. This represents the decorrelation time constant of the features. This represents the natural exponential function.

[0175] By using this weighted selection, the number of interference pairs is controlled between 3N and 5N, which ensures redundancy while avoiding excessive computational burden.

[0176] Accurate correction of atmospheric phase delay and orbital errors is crucial for deformation extraction. This embodiment employs a multi-scale separation strategy to progressively eliminate various errors: first, the GACOS global atmospheric model is used to remove large-scale tropospheric delay; then, spatiotemporal filtering is used to separate turbulent atmospheric signals from actual deformation. Addressing the significant impact of topography on atmospheric delay, the atmospheric phase delay processing flow includes:

[0177] Based on the aforementioned digital elevation model and interferometric phase map, regression coefficients are estimated through the phase-elevation statistical relationship in stable regions. Then, an elevation correlation model is used to calculate and remove topographically related atmospheric delay signals. The elevation correlation model is as follows:

[0178] ;

[0179] in, Indicates atmospheric delay, Indicates terrain elevation. Represents the regression coefficient. This represents the residual turbulent atmospheric phase component.

[0180] The identification and processing of geometric distortions are crucial for InSAR applications in mountainous areas. This embodiment proposes a multi-scale R-index model to accurately identify distortion regions such as overlapping and shadows. The R-index comprehensively considers terrain and radar observation geometry and is defined as:

[0181] ;

[0182] in, express Index value, Indicates the slope of the terrain. Indicates slope direction. Indicates the radar incident angle. Indicates the satellite azimuth angle.

[0183] when When it is determined to be a shaded area, Areas identified as overlapping zones are masked. To improve monitoring coverage and fully leverage the complementarity of the rising and falling track data, data from another track is preferentially used to supplement areas where single-track data is affected by distortion.

[0184] In this embodiment, the robust solution of temporal deformation is obtained by using an improved singular value decomposition method to solve the deformation time series. The processing flow includes:

[0185] Based on the optimized set of interferometric pairs, the observation equation is constructed as follows: ;

[0186] The deformation rate vector is solved using an iterative weighted least squares estimation method: ;

[0187] in, This represents the unwrapped phase vector. Represents the design matrix. Represents the deformation rate vector. Indicates observation noise. Represents the weight matrix. Represents the smoothing constraint matrix. Represents the regularization parameter. This indicates the matrix transpose.

[0188] After 35 iterations of calculations, the system can automatically identify abnormal observations and reduce their weight, ultimately achieving an annual deformation rate estimation accuracy better than 2 millimeters. To facilitate geological interpretation, the radar line of sight was finally shifted to the slope direction to support subsequent stability analysis and disaster early warning.

[0189] Step 4: Spatiotemporal fusion of multi-source deformation results:

[0190] Based on the horizontal displacement field and the temporal deformation field, spatial adaptive grid fusion, Gaussian process-based temporal interpolation, and dynamic allocation of fusion weights for optical and SAR data are performed under a unified spatiotemporal framework to obtain a spatiotemporally continuous comprehensive deformation field.

[0191] The core objective of this step is to deeply fuse the large deformation information obtained from optical offset tracking with the small deformation signals detected by InSAR within a unified spatiotemporal framework, in order to overcome the differences between the two technologies in terms of spatial coverage, temporal sampling, and deformation magnitude. By constructing an adaptive fusion framework, the complementary advantages of multi-source deformation data are fully utilized to generate a complete, continuous, and high-precision deformation field, providing a reliable data foundation for subsequent activity identification.

[0192] In terms of spatial adaptive mesh fusion, an adaptive meshing strategy is adopted to dynamically determine the optimal mesh size based on the spatial distribution density of deformation points and the deformation gradient characteristics. Traditional fixed mesh methods struggle to balance detail preservation and noise suppression; this embodiment uses the maximum curvature point method to identify the inflection point location of mesh refinement.

[0193] In this embodiment, the spatial adaptive mesh fusion processing flow includes:

[0194] The evaluation function values ​​corresponding to different grid sizes are calculated. The optimal grid size is determined by finding the optimal value of the evaluation function, thus achieving a balance between data coverage and spatial resolution. The corresponding evaluation function is as follows:

[0195] ;

[0196] in, Indicates the grid size. Indicates grid size The corresponding evaluation function value, This indicates the number of grid cells containing observation points. Indicates the total number of grid cells. Indicates the feature scale. Indicates the weighting coefficient. Represents the natural exponential function;

[0197] After determining the optimal grid size, within each grid cell, inverse-range weighted interpolation is used to fuse optical and InSAR deformations:

[0198] ;

[0199] ;

[0200] in, Represents the shape variables after fusion. Indicates the first The weight of each observation point Indicates the first The deformation of each observation point Indicates the first The observation accuracy of each observation point Indicates the first The distance from each observation point to the center of the grid. This represents the smoothing parameter.

[0201] This weighting scheme takes into account both differences in observation accuracy and ensures spatial continuity.

[0202] In Gaussian process-based time interpolation, to address the inconsistency between uneven sampling of optical data due to cloud interference and periodic revisits of SAR data, a Gaussian process-based time interpolation model is established to unify irregular sampling to a standard time series. The deformation time evolution can be expressed as:

[0203] ;

[0204] in, Indicates time The shape variable, Represents the mean function, Represents a Gaussian process. Represents the covariance function. and Represents the covariance function at different time points. For the Matérn core:

[0205] ;

[0206] in, Indicates the signal variance. Indicates the characteristic time scale.

[0207] Hyperparameters are determined by maximum likelihood estimation, enabling smooth interpolation and uncertainty quantification between observation points.

[0208] In the dynamic allocation of fusion weights for optical and SAR data based on deformation rate, terrain conditions, and data quality, the weights are dynamically adjusted according to the applicable technical conditions: optical data is prioritized when the deformation rate exceeds 100 mm / year, while InSAR is prioritized when it is below 50 mm / year; the InSAR weight is reduced when the terrain slope is greater than 35°, and the optical weight is reduced in densely vegetated areas. The weights of optical data and SAR data are calculated and determined using the following formulas:

[0209] ;

[0210] ;

[0211] in, Indicates the weight of optical data. Indicates SAR data weights, Indicates the deformation rate. Indicates the deformation rate threshold. Indicates the slope of the terrain. The characteristic parameter representing the deformation rate, Indicates slope characteristic parameters, and These represent the quality indicators for optical and SAR data, respectively.

[0212] In practical applications, the quality of the fusion results can be evaluated through cross-validation and uncertainty analysis. 10% of the observation points are randomly retained as validation samples, requiring the root mean square error between the fused and measured values ​​to be better than 1.5 times the accuracy of any single data source. Simultaneously, the uncertainty of the fusion results is estimated based on error propagation theory, generating a confidence level layer. Regions with uncertainty exceeding 50% of the deformation are marked as low confidence for cautious use in subsequent analyses. After the above processing, a comprehensive deformation field with complete spatial coverage and temporal continuity is finally formed, improving spatial coverage by over 40% and achieving a temporal resolution of 12 days, providing a solid data foundation for the accurate identification of active high-altitude loose bodies.

[0213] Step 5: Automatic identification and activity classification of deformation features based on machine learning:

[0214] Based on the comprehensive deformation field, feature extraction and dimensionality reduction are performed through principal component analysis, and unsupervised K-means clustering algorithm is used to automatically classify deformation patterns. The activity level is defined by combining the category conversion of pre-earthquake and post-earthquake deformation patterns, and deformation type map and activity level map of each observation point in the study area are obtained.

[0215] This step utilizes unsupervised machine learning methods to automatically classify and identify features in time-series deformation data, aiming to overcome the issues of strong subjectivity and low efficiency inherent in traditional manual interpretation. The core principle is to project high-dimensional time-series deformation data into a low-dimensional feature space, automatically identify different deformation patterns using cluster analysis, and define activity levels based on the transformation of deformation patterns before and after the earthquake. This provides an objective and quantitative basis for assessing the hazard of high-altitude loose bodies.

[0216] First, feature extraction and dimensionality reduction are performed using principal component analysis:

[0217] Standardize the temporal deformation data of the deformation matrix corresponding to the comprehensive deformation field:

[0218] ;

[0219] in, Indicates the first The observation point at the ... The deformation at each time sampling point Represents the deformable variable after standardization. Indicates the first The mean of the time-sampled shape variables for each observation point. Indicates the first Standard deviation of all time-sampled deformation variables at each observation point;

[0220] Construct the covariance matrix, perform eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​and eigenvectors. The covariance matrix is ​​as follows:

[0221] ;

[0222] in, Represents the covariance matrix. This represents the deformation matrix corresponding to the standardized deformation variables. Indicates matrix transpose. Indicates the number of observation points;

[0223] Projecting the deformation matrices corresponding to the standardized deformation variables onto the eigenvectors yields the principal components:

[0224] ;

[0225] in, Indicates the first Principal components, The eigenvalue obtained by eigendecomposition 1 eigenvector;

[0226] Calculate the cumulative variance contribution rate of the eigenvalues ​​and select the smallest eigenvalue that results in a cumulative variance contribution rate of 85% or higher. The value is used as the number of principal components, that is:

[0227] ;

[0228] in, The eigenvalue obtained by eigendecomposition 1 eigenvalue, Indicates the number of principal components. This represents the total number of eigenvalues.

[0229] Practice shows that the first 2-3 principal components can usually capture the main features of deformation time series. The first principal component reflects the overall deformation trend, the second principal component reflects seasonal changes, and the third principal component represents sudden changes.

[0230] Then, based on the dimensionality-reduced feature space, the K-means clustering algorithm is used to automatically identify the deformation type. Determining the number of clusters is a key parameter. In this embodiment, automatic optimization is achieved through the silhouette coefficient: a range of candidate cluster numbers is set, and K-means clustering is performed for each candidate number. The deformation of each observation point is treated as a sample point, the silhouette coefficient of all sample points is calculated, and the average value is taken. The candidate number that maximizes the average silhouette coefficient is selected as the final number of clusters. The formula for calculating the silhouette coefficient is as follows:

[0231] ;

[0232] in, Indicates the first Profile coefficients for each sample point Indicates the first The average distance from each sample point to other sample points in the same cluster. Indicates the first The average distance from each sample point to all sample points in the nearest neighbor cluster.

[0233] To improve clustering stability and convergence efficiency, this embodiment adopts an initial center selection strategy based on probability distribution. Specifically, during unsupervised K-means clustering initialization, the selection probability of each sample point is calculated, and the initial cluster centers are selected based on these probabilities. The probability calculation formula is as follows:

[0234] ;

[0235] in, Indicates the first sample points The probability of being selected as the next cluster center. Represents sample points Euclidean distance to the nearest cluster center Indicates the first sample points Euclidean distance to the nearest cluster center This indicates the number of observation points.

[0236] During the iterative update process of the unsupervised K-means clustering algorithm, the cluster centers are updated according to the weighted average of the member points:

[0237] ;

[0238] in, Indicates the first The coordinates of the center point of each cluster Indicates assignment to the first The weights for each sample point are calculated based on its measurement uncertainty and reference accuracy:

[0239] ;

[0240] in, Indicates the first Measurement uncertainty for each sample point Indicates the reference precision.

[0241] Furthermore, by analyzing the transformation of deformation patterns before and after the earthquake, this embodiment defines four types of activity characteristics: Type A (stable before earthquake - stable after earthquake) is a non-sensitive area; Type B (active before earthquake - stable after earthquake) is a self-stabilizing area; Type C (active before earthquake - active after earthquake) is a continuously active area; and Type D (stable before earthquake - active after earthquake) is a seismically activated area. Type D areas are of significant disaster prevention importance; when the deformation acceleration exceeds three standard deviations of the background value, lasts for more than 30 days, and affects an area greater than 1 square kilometer, it is identified as a high-risk active high-level loose body. This data-driven automatic identification method is more than 20 times more efficient than manual interpretation, with an identification accuracy of 92%.

[0242] In this embodiment, it also includes:

[0243] Step 6: Construction of high-level loose body advanced identification model and generation of early warning information:

[0244] Based on the activity grading map, meteorological observation data, and comprehensive deformation field, quantitative analysis and time-delay correlation analysis of multi-factor triggering mechanisms are performed. Critical state identification is performed based on the Fukuzono-Voight model, and early warning information is generated according to the preset graded early warning criteria.

[0245] This step aims to construct an advanced identification model that comprehensively considers deformation evolution patterns, external triggering factors, and environmental conditions, in order to achieve quantitative evaluation and early warning of the entire process of high-altitude loose bodies from a stable state to instability. The core innovation of this model lies in establishing a quantitative response relationship between deformation acceleration and external triggering factors, revealing the climate driving mechanism through time-lag correlation analysis, and determining graded early warning thresholds based on critical state theory, thereby advancing the warning time by 10–20 days and providing a scientific basis for disaster prevention and control decisions.

[0246] Quantitative analysis of multi-factor triggering mechanisms is the theoretical basis for advanced identification. This embodiment comprehensively considers four main triggering factors: precipitation, temperature, earthquakes, and human activities, and constructs a deformation response model. For climatic factors, a sliding time window is used to calculate the cross-correlation function between deformation rate and meteorological elements to identify the optimal response time lag.

[0247] The precipitation-triggered effect is characterized by the relationship between cumulative precipitation and deformation acceleration: a quantitative relationship between deformation acceleration and cumulative precipitation is established based on a sliding time window, as shown in the following model:

[0248] ;

[0249] in, express Deformation acceleration at time t, Indicates background acceleration. This represents the cumulative precipitation within the sliding time window. Indicates the time lag in precipitation response. This represents the precipitation impact coefficient.

[0250] The effect of temperature is mainly reflected in the weakening of the strength of loose materials by freeze-thaw cycles. It is quantified using the positive accumulated temperature index, which quantifies the intensity of freeze-thaw cycles. The calculation formula is as follows:

[0251] ;

[0252] in, Indicates positive accumulated temperature index, This represents the average daily temperature. Indicates the melting threshold (0℃). This indicates the summation period, ranging from 30 to 60 days.

[0253] The assessment of earthquake triggering effects includes: using a modified Newmark displacement model, combined with peak ground acceleration, duration, and initial stability of high-risk areas, to assess the cumulative impact of earthquake events on the stability of loose bodies and generate earthquake impact factors.

[0254] Time-delay correlation analysis is used to reveal the time delay characteristics between external factors and deformation response, and is key to achieving early warning. The process of time-delay correlation analysis includes:

[0255] Based on the comprehensive deformation field and meteorological observation data, the cross-correlation coefficient between the deformation time series and the triggering factor time series at different precipitation response lags is calculated. The optimal precipitation response lag is determined by finding the maximum cross-correlation coefficient. The formula for calculating the cross-correlation coefficient is as follows:

[0256] ;

[0257] in, Represents the cross-correlation coefficient. Indicates the deformation time sequence, Indicates the timing of triggering factors. and These represent the corresponding means;

[0258] A multi-factor comprehensive triggering index is constructed to normalize and weight the influence of different triggering factors under the optimal precipitation response time lag:

[0259] ;

[0260] in, Indicates the comprehensive trigger index, and These represent the normalized precipitation and temperature indices, respectively. Indicates the earthquake impact factor. These represent the corresponding weights.

[0261] Critical state identification employs a dual criterion of acceleration and displacement rate threshold. The critical state identification process includes: based on the comprehensive deformation field, calculating the reciprocal of the deformation rate time series data of the target region within a sliding time window, and performing linear fitting using the Fukuzono-Voight model. When the slope of the fitted line is negative and the goodness of fit is greater than the goodness of fit threshold, it indicates that the deformation is accelerating and the system is tending to become unstable. The Fukuzono-Voight model is as follows:

[0262] ;

[0263] in, Indicates the deformation rate. Indicates time, and This represents the model parameters obtained by fitting deformation data within a sliding window.

[0264] In practical applications, warning trigger conditions are set, and an warning is triggered when any two conditions are met:

[0265] (1) Deformation rate More than 100 mm / day;

[0266] (2) The acceleration exceeds three standard deviations of the background value;

[0267] (3) speed The linear fit has a negative slope and a goodness of fit .

[0268] The tiered early warning system integrates three dimensions—hazard, vulnerability, and risk level—and classifies early warnings into four levels:

[0269] 1) A blue alert (attention level) is issued when the deformation rate reaches 50 mm / day or the acceleration exceeds twice the standard deviation of the background value;

[0270] 2) A yellow alert (warning level) is triggered when the deformation rate exceeds 100 mm / day for three consecutive days;

[0271] 3) Orange alert (alert level) is triggered when any critical identification condition is met and the affected area exceeds 50,000 square meters;

[0272] 4) Red alert (alarm level) is triggered when two or more critical conditions are met or when experts determine that the situation is about to become unstable.

[0273] The early warning information includes the hazard level, estimated occurrence time, impact range, and response recommendations. Validated through historical case studies, the model boasts an 85% accuracy rate, a false alarm rate of less than 20%, and an average warning lead time of 15 days, significantly improving the prevention and control capabilities against high-altitude loose soil disasters.

[0274] In summary, the high-altitude loose body multi-scale deformation monitoring method based on optical and SAR provided in this embodiment constructs a complete technical chain of "data acquisition - feature extraction - fusion analysis - advanced identification." By fully leveraging the complementary advantages of optical and SAR remote sensing data, it achieves quantitative monitoring and early warning of the entire process of loose bodies in high-altitude areas from stability to instability. This not only effectively solves the key problems of existing technologies in data source fusion, feature recognition, and early warning timeliness, but also provides a new technical path for the systematic prevention and control of geological disasters in high-altitude areas. Specifically, this embodiment adopts a layered and progressive processing strategy. In the data acquisition stage, a standardized processing flow for multi-source remote sensing data is established to ensure the comparability and consistency of data from different sensors and different time phases. Subsequently, a modified pixel offset tracking algorithm is used to capture large-scale rapid deformation, while optimized temporal InSAR technology is used to detect early micro-deformation signals. These two technical paths are advanced in parallel to obtain complementary deformation information. Building upon this foundation, an adaptive mesh fusion algorithm is employed to integrate heterogeneous deformation results into a unified spatiotemporal reference frame. Subsequently, machine learning algorithms are used to automatically classify deformation features, and finally, an activity-based early warning model is established by combining external triggering factor analysis. From a system architecture perspective, this embodiment constructs a modular processing framework encompassing five core components: data input, deformation detection, data fusion, feature recognition, and early warning output. These modules transmit information through standardized data interfaces, ensuring system scalability and providing flexible technical support for practical applications.

[0275] Compared with the prior art, this embodiment has the following advantages:

[0276] (1) The fusion of multi-source remote sensing data significantly improves the integrity and continuity of deformation monitoring;

[0277] Compared to existing technologies that rely on a single remote sensing data source, resulting in monitoring blind spots, this embodiment achieves simultaneous monitoring of large deformations (meter-level) and small deformations (millimeter-level) by deeply fusing optical and SAR remote sensing data. This effectively solves the technical bottlenecks of optical technology's inability to capture early, minute deformations and SAR technology's loss of correlation in rapidly deforming areas. Through adaptive grid fusion and dynamic weight allocation strategies, this embodiment increases the spatial coverage of deformation monitoring from 60% in existing technologies to over 85%, and the temporal resolution from a monthly scale to 12 days. Experimental results show that, under complex mountainous terrain conditions, this embodiment can completely track the entire evolution of high-altitude loose bodies from initial creep and accelerated deformation to final instability, while existing single-technology methods can only capture deformation information at certain stages, easily missing the optimal early warning opportunity.

[0278] (2) Automated identification technology has greatly improved the efficiency and objectivity of deformation classification;

[0279] To address the issues of high subjectivity and low efficiency caused by existing technologies relying on manual interpretation and fixed thresholds, this embodiment constructs an automatic identification framework based on machine learning, achieving objective quantitative classification of deformation features. By using principal component analysis to reduce the dimensionality of complex temporal deformation data to 2-3 key features, and combining this with unsupervised clustering algorithms to automatically identify different deformation patterns, the processing efficiency is more than 20 times higher than traditional manual methods, reducing the processing time for a single image from several hours to less than 10 minutes. More importantly, this embodiment eliminates the uncertainty of human experience-based judgment, achieving over 98% consistency in repeated processing results for the same dataset, while traditional methods can result in up to 30% differences between different interpreters, ensuring the repeatability and comparability of monitoring results.

[0280] (3) The pre-earthquake and post-earthquake mode conversion analysis enabled the quantitative differentiation of the triggering mechanism;

[0281] Compared to existing technologies that cannot effectively distinguish between different triggering factors, this embodiment innovatively proposes a four-category activity classification method based on the transformation of deformation patterns before and after an earthquake. This method accurately identifies earthquake-induced D-type activation areas and achieves a quantitative separation of the contributions of earthquakes and climate factors to deformation. By comparing and analyzing the deformation evolution characteristics 18 months before and after the earthquake, this embodiment can accurately delineate high-risk areas directly affected by the earthquake, achieving an accuracy rate of 92%, while existing methods can only provide vague qualitative judgments. This quantitative differentiation capability enables disaster prevention departments to formulate differentiated prevention and control strategies for different triggering mechanisms, such as strengthening engineering management in earthquake-activated areas and enhancing monitoring and early warning in climate-sensitive areas, significantly improving the targeting and effectiveness of disaster prevention and control.

[0282] (4) The advanced identification model improves the timeliness of early warning to a practical level;

[0283] The multi-factor coupled early warning model established in this embodiment reveals the lag effect of precipitation (20±5 days) and temperature (10±3 days) on deformation through time-lag correlation analysis. Combined with the critical state identification criterion, the early warning lead time is extended from the existing 3-5 days to 10-20 days, gaining valuable time for emergency evacuation and engineering response. Compared with the existing technology's 85% false alarm rate and passive post-event analysis, the hierarchical early warning system of this embodiment achieves a hit rate of 85% and controls the false alarm rate to within 20%, basically meeting the requirements for practical application. In particular, in the application case after the 2017 Linzhi earthquake, this embodiment successfully issued an early warning of a large debris flow event in a river basin in southeastern Tibet 18 days in advance, while traditional monitoring methods only detected the anomaly 2 days before the disaster, fully demonstrating the significant advantages of this embodiment in early warning.

[0284] The effectiveness of the method of this invention is verified below using a river basin in southeastern Tibet as a typical research area:

[0285] (1) Overview of the study area and data acquisition:

[0286] Located on the southeastern edge of the Qinghai-Tibet Plateau, the basin ranges in altitude from 2,800 to 6,500 meters, characterized by dramatic topography, well-developed glaciers, and widespread loose deposits. Since 1950, this region has experienced seven large-scale river-blocking events. Activity significantly increased after the 6.9-magnitude earthquake in Nyingchi on November 18, 2017. On October 17 and 29, 2018, two consecutive ice avalanche-debris flow chains caused the Yarlung Tsangpo River to be cut off for dozens of hours, resulting in direct economic losses exceeding 800 million yuan. Please refer to Table 1. In accordance with the data acquisition standards of this invention, the system collected multi-source remote sensing data from October 2016 to February 2019, including 28 Planet satellite images (resolution 3.7 meters, covering key periods), 6 RapidEye images (resolution 5 meters, used for verification), 15 Sentinel-2 images (resolution 10 meters, as the main optical data source), and 42 Sentinel-1 ascending orbit images and 46 descending orbit images of SAR data, ensuring a complete observation sequence before and after the earthquake. Data quality control strictly implemented the screening criteria of cloud coverage rate below 30% and coherence greater than 0.3.

[0287] Table 1. Statistical Table of Multi-Source Remote Sensing Data Acquisition and Quality Control

[0288]

[0289] (2) Implementation and verification of optical large deformation tracking

[0290] The improved COSI-Corr algorithm of this invention was used to process optical images, focusing on the deformation characteristics changes before and after the November 2017 earthquake. Taking Sentinel-2 image pairs from November 11th and 19th, 2017 as examples, sub-pixel registration and offset calculation were performed using an adaptive window strategy (initial 64×64 pixels, final 32×32 pixels, step size 16 pixels) to obtain the east-west and north-south displacement fields. After applying polynomial surface fitting to eliminate orbital errors, the effective pixel ratio increased from the original 72% to 86%, and the noise level decreased by 65%. Figure 2 As shown, the maximum horizontal displacement detected in the upstream glacier area of ​​the basin reached 15.3 m / year, and the average movement rate of loose deposits in the main channel was 8.7 m / year. In particular, in channel II and IV, the accelerated material transport process after the earthquake was clearly identified. The displacement rate of channel IV increased from 3.2 m / year before the earthquake to 12.8 m / year after the earthquake. This acceleration process is highly consistent in time and space with the debris flow blocking the river that actually occurred on December 21, 2017, verifying the predictive ability of the method of the present invention.

[0291] (3) Implementation of time-series InSAR small deformation detection

[0292] Optimized SBAS-InSAR technology was used to process Sentinel-1 data from before the earthquake (October 2016 to November 17, 2017) and after the earthquake (November 19, 2017 to February 2019). An R-index model was used to identify and remove geometrically distorted regions covering 23% of the study area, ensuring the reliability of the deformation results. 156 interferometric pairs were constructed before the earthquake, and 198 interferometric pairs were constructed after the earthquake. The time baseline was controlled within 180 days, and the vertical baseline was limited to within 40% of the critical baseline. After atmospheric correction and orbital error removal, a deformation rate field with millimeter-level accuracy was obtained, with an average deformation rate of 12.5 mm / year before the earthquake and increasing to 35.8 mm / year afterward. Creep signals of multiple newly formed cracks were detected in the bedrock slope area in the middle and lower reaches of the watershed, with deformation rates between 5 and 15 mm / year. These minute deformations are undetectable by optical methods, highlighting the necessity of multi-source data fusion. By combining the data from the lifting and lowering rails, the LOS deformation was decomposed into three-dimensional deformation components. It was found that the vertical settlement was mainly concentrated in the glacial deposit area, with a maximum settlement rate of -85 mm / year.

[0293] (4) Spatiotemporal fusion of multi-source data and deformation field reconstruction

[0294] Please see Figure 3 The adaptive mesh fusion strategy of this invention determines the optimal mesh size to be 85 meters using the maximum curvature point method. At this resolution, both local deformation details are preserved and noise is effectively suppressed. Please refer to [link / reference]. Figure 4 During spatial fusion, the weights of optical and InSAR data are dynamically adjusted based on deformation gradient and data quality. In rapid deformation regions with deformation rates greater than 100 mm / year, the weight of optical data reaches 0.7; in slow deformation regions with deformation rates less than 50 mm / year, the InSAR weight increases to 0.8. Gaussian interpolation is used in the time dimension to unify irregular sampling into a standard 12-day time series, with interpolation uncertainty controlled within 1.5 times the original measurement error. Please refer to [link / reference]. Figure 5 The spatial coverage of the fused deformation field increased from 62% to 87% of the single InSAR and from 54% to 87% of the single optical field, achieving near-complete coverage of deformation information in the study area, with only a small amount of data gaps in extremely steep cliffs and permanent shadow areas.

[0295] (5) Automatic identification of deformation features and activity classification

[0296] The machine learning framework of this invention was applied to automatically classify the fused temporal deformation data. Principal component analysis was used to extract the first two principal components (cumulative variance contribution rate 87.6%). The first principal component reflects the long-term deformation trend, and the second principal component reflects seasonal fluctuations (see Tables 2 and 3). Based on silhouette coefficient optimization, the number of clusters was determined to be 2, automatically classifying the 13,847 valid deformation points into two categories: stable (Type I, accounting for 68%) and active (Type II, accounting for 32%). Further, based on the pre- and post-earthquake category transformations, four types of activity were identified: Type A stable zone (8,234 points, 59.5%), Type B self-stabilizing zone (1,126 points, 8.1%), Type C persistently active zone (2,308 points, 16.7%), and Type D earthquake-activated zone (2,179 points, 15.7%) (see Tables 2 and 3). Figure 6 The Type D region, covering an area of ​​4.2 square kilometers, is mainly distributed on loose deposits within 60 kilometers of the epicenter. These areas experienced three large-scale instability events within 18 months after the earthquake, confirming the accuracy of the classification results. The entire automatic identification process took only 45 minutes, while traditional manual interpretation requires 3-5 working days.

[0297] Table 2 Statistical Table of Deformation Characteristics of Different Activity Types

[0298]

[0299] Table 3 Deformation Characteristic Analysis Table

[0300]

[0301] Note: ± represents the standard deviation; the acceleration coefficient is the ratio of the deformation rate after the earthquake to the deformation rate before the earthquake; the deformation direction dominance refers to the direction in which the principal component contribution rate is >60%.

[0302] (6) Application of advanced identification model and verification of early warning effect

[0303] A multi-factor advance identification model for the study area was constructed, and the optimal precipitation response time lags for meteorological factors were determined through sliding correlation analysis: 18 days for temperature and 22 days for precipitation. (Please refer to...) Figures 7 to 9 , Figure 8 The bar chart represents daily precipitation, and the dashed line represents daily average temperature. Before the two disaster events in October 2018, the deformation-type D region experienced a 15-day period of high temperatures with daily average temperatures above 10°C and a heavy rainfall event with accumulated precipitation reaching 186 mm. The trigger index exceeded the warning threshold on October 5th. The Fukuzono-Voight model was used to analyze the critical state of key points. When the reciprocal of the rate approached zero and the goodness of fit R² > 0.85, the system triggered an orange alert on October 12th, with a predicted instability window of October 15th-20th. The actual disaster occurred on October 17th, with a warning lead time of 5 days and a prediction error of only 2 days. Through retrospective testing of all deformation events between 2016 and 2019, the method of this invention successfully issued early warnings in 9 out of 11 events (hit rate 81.8%), with 2 false alarms (false alarm rate 18.2%), and an average early warning lead time of 12.5 days, which is a significant improvement compared to the 3-day early warning period of traditional monitoring methods. This provided sufficient time for the safe evacuation of more than 3,800 residents in 17 downstream villages, fully verifying the practical value of this embodiment in the early identification of active high-altitude loose bodies.

Claims

1. A method for monitoring multi-scale deformation of high-altitude loose bodies based on optics and SAR, characterized in that, The method includes: Acquire multi-source remote sensing data of the study area, including multi-temporal optical satellite imagery, multi-temporal SAR satellite imagery, and auxiliary data, including digital elevation models, meteorological observation data, and earthquake catalog data; Radiometric calibration, atmospheric correction, and terrain correction are performed on the multi-temporal optical satellite images. Precise orbit correction, thermal noise removal, radiometric calibration, and adaptive filtering are performed on the multi-temporal SAR satellite images. All auxiliary data are uniformly resampled to a consistent spatiotemporal reference to obtain a standardized multi-source remote sensing dataset with spatiotemporal consistency. The improved frequency domain cross-correlation algorithm is used to process multi-temporal optical satellite images in the multi-source remote sensing dataset, including image fine registration, pixel offset calculation, adaptive window adjustment and system error correction, to obtain the horizontal displacement field of the study area. Based on the multi-temporal SAR satellite images and digital elevation models in the multi-source remote sensing dataset, the temporal deformation field of the study area is obtained through optimized interferometric pair selection, atmospheric phase delay and orbital error correction, geometric distortion identification and processing, and robust temporal deformation solution. Based on the horizontal displacement field and temporal deformation field, spatial adaptive grid fusion, Gaussian process-based temporal interpolation, and dynamic allocation of fusion weights for optical and SAR data are performed under a unified spatiotemporal framework to obtain a spatiotemporally continuous comprehensive deformation field. Based on the comprehensive deformation field, feature extraction and dimensionality reduction are performed through principal component analysis, and unsupervised K-means clustering algorithm is used to automatically classify deformation patterns. The activity level is defined by combining the category conversion of pre-earthquake and post-earthquake deformation patterns, and deformation type map and activity level map of each observation point in the study area are obtained.

2. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, When performing terrain correction on the multi-temporal optical satellite imagery, the C-correction model is used, and its expression is: ； in, This represents the radiance value after terrain correction. This represents the original radiance value. Indicates the zenith angle of the sun. Indicates the angle of incidence. This represents the correction parameters determined through statistical regression on flat terrain areas; When performing adaptive filtering on the multi-temporal SAR satellite imagery, the size of the filtering window is dynamically calculated and adjusted based on the local coherence coefficient. The calculation formula is as follows: ； in, Indicates the size of the filtering window. and These represent the preset minimum and maximum window sizes, respectively. This represents the local coherence coefficient.

3. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, The pixel offset calculation process includes: Two-dimensional fast Fourier transform is performed on the corresponding windows of the main image and the slave image to obtain the spectra of the main image and the slave image respectively. The normalized cross power spectrum is calculated based on the spectra of the main image and the slave image. The inverse Fourier transform is performed on the normalized cross power spectrum to obtain the spatial correlation surface. The pixel-level displacement is determined by finding the peak value of the correlation surface. Within the neighborhood of the peak position determined by the pixel-level displacement, the discrete related surface intensity values ​​are fitted by least squares using a two-dimensional Gaussian surface model to solve for the sub-pixel peak position as the final pixel offset. The formula for calculating the normalized cross-power spectrum is as follows: ； in, Represents the normalized cross-power spectrum. and Represents frequency domain coordinate variables, Represents the spectrum of the main image. Indicates the spectrum of the image. express The complex conjugate, Indicates a small quantity to prevent division by zero; The two-dimensional Gaussian surface model is as follows: ； in, Indicates the relevant surface strength value. and Represents spatial domain coordinate variables, This is the sub-pixel peak position. Indicates peak amplitude. and These represent the Gaussian distributions in... and Standard deviation of direction Indicates the background noise level. This represents the natural exponential function.

4. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, The adaptive window adjustment process includes: During the pixel offset calculation process, the local standard deviation of each potential window position is calculated in real time, and its deformation gradient magnitude is estimated. Based on the local standard deviation and deformation gradient magnitude, the window size for cross-correlation calculation at that position is dynamically calculated and determined. The calculation formula is as follows: ； in, This indicates the adaptively adjusted window size. Indicates the base window size. Indicates local standard deviation. This represents the minimum standard deviation threshold. Represents the deformation gradient mode. Indicates the maximum expected deformation. and Indicates the weighting coefficient; The process for correcting system errors includes: After calculating the initial displacement field, a nonlocal mean filtering algorithm is used to filter the initial displacement field. The nonlocal mean filtering algorithm is as follows: ； in, Indicates the filtered first... Pixel offset of each pixel , and Represented by pixels and The neighborhood block centered on, Indicates the filter strength parameter. Represents the normalization factor. Indicates the first before filtering The pixel offset of each pixel.

5. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, The preferred processing procedure for interference pairs includes: Based on multi-temporal SAR satellite imagery, for all possible combinations of interferometric pairs, the weight of each interferometric pair is calculated based on its spatial baseline, temporal baseline, and predicted coherence. An optimized set of interferometric pairs is then constructed based on the weights for robust solution of temporal deformation. The weights of the interference pairs are calculated using the following formula: ； in, Indicates the weight of the interference pair. Indicates the coherence of the prediction. Represents the vertical spatial baseline. Indicates the critical baseline. Indicates the time baseline of the interferometer pair. This represents the decorrelation time constant of the features. Represents the natural exponential function; The atmospheric phase delay processing procedure includes: Based on the aforementioned digital elevation model and interferometric phase map, regression coefficients are estimated through the phase-elevation statistical relationship in stable regions. Then, an elevation correlation model is used to calculate and remove topographically related atmospheric delay signals. The elevation correlation model is as follows: ； in, Indicates atmospheric delay, Indicates terrain elevation. Represents the regression coefficient. Represents the residual atmospheric phase component; The geometric distortion recognition process includes: Based on the digital elevation model and SAR satellite orbit parameters, the elevation of each pixel is calculated. Exponential value, when When it is determined to be a shaded area, If an area is identified as an overlay area, the identified shadow and overlay areas are masked out. The The formula for calculating the index value is as follows: ； in, express Index value, Indicates the slope of the terrain. Indicates slope direction. Indicates the radar incident angle. Indicates the satellite azimuth angle.

6. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 5, characterized in that, The processing flow for the time-series deformation robustness solution includes: Based on the optimized set of interferometric pairs, the observation equation is constructed as follows: ; The deformation rate vector is solved using an iterative weighted least squares estimation method: ; in, This represents the unwrapped phase vector. Represents the design matrix. Represents the deformation rate vector. Indicates observation noise. Represents the weight matrix. Represents the smoothing constraint matrix. Represents the regularization parameter. This indicates the matrix transpose.

7. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, The spatial adaptive mesh fusion processing flow includes: The evaluation function values ​​corresponding to different grid sizes are calculated. The optimal grid size is determined by finding the optimal value of the evaluation function, thus achieving a balance between data coverage and spatial resolution. The corresponding evaluation function is as follows: ； in, Indicates the grid size. Indicates grid size The corresponding evaluation function value, This indicates the number of grid cells containing observation points. Indicates the total number of grid cells. Indicates the feature scale. Indicates the weighting coefficient. Represents the natural exponential function; After determining the optimal grid size, within each grid cell, inverse-range weighted interpolation is used to fuse optical and InSAR deformations: ； ； in, Represents the shape variables after fusion. Indicates the first The weight of each observation point Indicates the first The deformation of each observation point Indicates the first The observation accuracy of each observation point Indicates the first The distance from each observation point to the center of the grid. Indicates the smoothing parameter; In the Gaussian process-based time interpolation, the deformation time evolution is expressed as: ； in, Indicates time The shape variable, Represents the mean function, Represents a Gaussian process. Represents the covariance function. and Represents the covariance function at different time points. For the Matérn core: ； in, Indicates the signal variance. Indicates the time scale of the feature; In the dynamic allocation of fusion weights for optical and SAR data based on deformation level, terrain conditions, and data quality, the weights for optical and SAR data are calculated and determined according to the following formulas: ； ； in, Indicates the weight of optical data. Indicates SAR data weights, Indicates the deformation rate. Indicates the deformation rate threshold. Indicates the slope of the terrain. The characteristic parameter representing the deformation rate, Indicates slope characteristic parameters, and These represent the quality indicators for optical and SAR data, respectively.

8. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, Feature extraction and dimensionality reduction are performed using principal component analysis, including: Standardize the temporal deformation data of the deformation matrix corresponding to the comprehensive deformation field: ； in, Indicates the first The observation point at the ... The deformation at each time sampling point Represents the deformable variable after standardization. Indicates the first The mean of the time-sampled shape variables for each observation point. Indicates the first Standard deviation of all time-sampled deformation variables at each observation point; Construct the covariance matrix, perform eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​and eigenvectors. The covariance matrix is ​​as follows: ； in, Represents the covariance matrix. This represents the deformation matrix corresponding to the standardized deformation variables. Indicates matrix transpose. Indicates the number of observation points; Projecting the deformation matrices corresponding to the standardized deformation variables onto the eigenvectors yields the principal components: ； in, Indicates the first Principal components, The eigenvalue obtained by eigendecomposition 1 eigenvector; Calculate the cumulative variance contribution rate of the eigenvalues ​​and select the smallest eigenvalue that results in a cumulative variance contribution rate of 85% or higher. The value is used as the number of principal components, that is: ； in, The eigenvalue obtained by eigendecomposition 1 eigenvalue, Indicates the number of principal components. Indicates the total number of eigenvalues; When using the unsupervised K-means clustering algorithm, the process for determining the optimal number of clusters includes: Define a range for the number of candidate clusters. For each candidate, perform K-means clustering, treating the deformation of each observation point as a sample point. Calculate the silhouette coefficient of all sample points and take the average. Select the candidate that maximizes the average silhouette coefficient as the final number of clusters. The formula for calculating the silhouette coefficient is as follows: ； in, Indicates the first Profile coefficients for each sample point Indicates the first The average distance from each sample point to other sample points in the same cluster. Indicates the first The average distance from each sample point to all sample points in the nearest neighbor cluster; During unsupervised K-means clustering initialization, the selection probability of each sample point is calculated, and the initial cluster centers are selected based on the probability. The formula for calculating the probability is as follows: ； in, Indicates the first sample points The probability of being selected as the next cluster center. Represents sample points Euclidean distance to the nearest cluster center Indicates the first sample points Euclidean distance to the nearest cluster center Indicates the number of observation points; During the iterative update process of the unsupervised K-means clustering algorithm, the cluster centers are updated according to the weighted average of the member points: ； in, Indicates the first The coordinates of the center point of each cluster Indicates assignment to the first The weights for each sample point are calculated based on its measurement uncertainty and reference accuracy: ； in, Indicates the first Measurement uncertainty for each sample point Indicates the reference precision.

9. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to claim 1, characterized in that, The activity levels include the following four categories: Non-sensitive areas: Stable before the earthquake - Stable after the earthquake; Self-stabilizing zone: pre-earthquake activity - post-earthquake stability; Continuous activity zone: pre-earthquake activity - post-earthquake activity; Earthquake activation zone: stable before earthquake - active after earthquake.

10. The method for multi-scale deformation monitoring of high-altitude loose bodies based on optics and SAR according to any one of claims 1 to 9, characterized in that, The method further includes: Based on the activity grading map, meteorological observation data, and comprehensive deformation field, quantitative analysis and time-delay correlation analysis of multi-factor triggering mechanisms are performed. Critical state identification is performed based on the Fukuzono-Voight model, and early warning information is generated according to the preset graded early warning criteria. The quantitative analysis includes: precipitation triggering effect modeling, temperature triggering effect modeling, and earthquake triggering effect assessment; The precipitation triggering effect modeling includes: establishing a quantitative relationship between deformation acceleration and cumulative precipitation based on a sliding time window, as shown in the following model: ； in, express Deformation acceleration at time t, Indicates background acceleration. This represents the cumulative precipitation within the sliding time window. Indicates the time lag in precipitation response. Indicates the precipitation impact coefficient; Modeling the temperature-triggered effect includes calculating the positive accumulated temperature index used to quantify the intensity of freeze-thaw cycles. ； in, Indicates positive accumulated temperature index, This represents the average daily temperature. Indicates the melting threshold. Indicates the time period for summation; The assessment of earthquake triggering effects includes: using a modified Newmark displacement model, combined with peak ground acceleration, duration, and initial stability of high-risk areas, to assess the cumulative impact of earthquake events on the stability of loose bodies and generate earthquake impact factors; The process of the time-delay correlation analysis includes: Based on the comprehensive deformation field and meteorological observation data, the cross-correlation coefficient between the deformation time series and the triggering factor time series at different precipitation response lags is calculated. The optimal precipitation response lag is determined by finding the maximum cross-correlation coefficient. The formula for calculating the cross-correlation coefficient is as follows: ； in, Represents the cross-correlation coefficient. Indicates the deformation time sequence, Indicates the timing of triggering factors. and These represent the corresponding means; A multi-factor comprehensive triggering index is constructed to normalize and weight the influence of different triggering factors under the optimal precipitation response time lag: ； in, Indicates the comprehensive trigger index, and These represent the normalized precipitation and temperature indices, respectively. Indicates the earthquake impact factor. These represent the corresponding weights; The critical state identification process includes: based on the comprehensive deformation field, calculating the reciprocal of the deformation rate time series data of the target region within a sliding time window, and performing linear fitting using the Fukuzono-Voight model. When the slope of the fitted line is negative and the goodness of fit is greater than the goodness of fit threshold, it indicates that the deformation is accelerating and the system is tending to become unstable. The Fukuzono-Voight model is as follows: ； in, Indicates the deformation rate. Indicates time, and This represents the model parameters obtained by fitting deformation data within a sliding window.

Images