High-position loose body multi-scale deformation monitoring method based on optics and SAR (Synthetic Aperture Radar)

Abstract

The invention 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 problems of poor continuity and accuracy of the existing method, and the scheme mainly comprises the following steps: obtaining a multi-temporal optical image, an SAR image and auxiliary data, and carrying out standardized preprocessing; extracting a horizontal displacement field from the optical image through an improved frequency domain cross-correlation algorithm, and obtaining a time sequence deformation field from the SAR image by using an optimized time sequence InSAR technology; performing spatial adaptive grid fusion and Gaussian process-based time interpolation on the two deformation fields, dynamically distributing fusion weights according to deformation magnitude, topographic conditions and data quality, and generating a time-space continuous comprehensive deformation field; automatic classification of deformation modes is realized through principal component analysis and a K-means clustering algorithm, activity levels are defined in combination with conversion of the deformation modes before and after an earthquake, and a deformation type graph and an activity classification graph are output. According to the invention, the continuity and accuracy of high-position loose body monitoring are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of remote sensing monitoring, and particularly relates to a high-position loose body multi-scale deformation monitoring method based on optics and SAR. BACKGROUND

[0002] With the continuous aggravation of global warming, the retreat of glaciers in high mountain areas and the acceleration of permafrost melting process, a large amount of loose material is accumulated in high and steep slope areas to form potential disaster sources. These high-position loose bodies are extremely easy to be transformed into a disaster chain under external factors, which poses a serious threat to the downstream area. However, the traditional geological disaster investigation method faces difficulties in implementation, high cost, limited coverage and other problems in high-altitude complex terrain areas, and it is difficult to meet the needs of continuous monitoring in a large range.

[0003] With the development of remote sensing technology, optical and synthetic aperture radar (SAR) remote sensing has gradually been applied to geological disaster identification and monitoring due to its advantages of large range, periodicity and non-contact observation. At present, the deformation monitoring of high-position loose bodies mainly relies on two types of technologies: one is the pixel offset tracking (POT) method based on optical images, which can extract the horizontal displacement of the ground surface, and the precision can reach the pixel level. However, this method is obviously restricted by cloudy and rainy weather, and the cost of obtaining high-resolution images is relatively high. The other is the synthetic aperture radar interferometric (InSAR) technology, especially the time series InSAR method (such as PS-InSAR and SBAS-InSAR), which has the ability to detect millimeter-level deformation and the advantage of all-weather observation. However, in the high mountain snow and vegetation coverage area, the signal decorrelation is easy to occur, and the overlap and shadow effect caused by steep terrain also limits the integrity of the deformation information extraction. In the existing research, the technical scheme close to the present application still has obvious limitations. On the one hand, most methods rely on a single remote sensing data source, such as monitoring glacier movement only by optical images or capturing deformation only by InSAR technology, which leads to monitoring blind spots in rapidly deforming areas or complex environments. On the other hand, although existing research attempts to combine optical and SAR data, it mostly stays at the level of simple comparison or superposition, and fails to establish an effective spatio-temporal fusion model to solve the inherent differences between the two types of data in spatial coverage, temporal sampling and deformation level. How to dynamically fuse multi-source information according to the deformation stage, terrain characteristics and data quality to generate a seamless, continuous and high-credibility comprehensive deformation field is still a technical problem to be solved. In addition, the deformation feature recognition is mostly based on manual interpretation or fixed threshold segmentation, lacking automated and intelligent classification means, and it is difficult to accurately capture the nonlinear behavior in the spatio-temporal evolution of deformation. In addition, the existing methods are mostly post-analysis, and the ability to identify in advance and analyze the coupling mechanism of multiple factors is insufficient, and a reliable quantitative prediction model has not been established.

[0004] Therefore, it is urgent to develop a technology method capable of effectively fusing multi-source remote sensing data and realizing automatic identification of deformation features, so as to improve the continuity and accuracy of loose body deformation monitoring in high mountain areas. SUMMARY

[0005] The present application aims to solve the problem that the prior art cannot realize continuous high-precision monitoring and automatic identification of the whole process from creep to instability of high-position loose bodies, and proposes a multi-scale deformation monitoring method for high-position loose bodies based on optics and SAR.

[0006] The technical scheme adopted by the present application to solve the above technical problems is:

[0007] The multi-scale deformation monitoring method for high-position loose bodies based on optics and SAR comprises:

[0008] Multi-source remote sensing data of a research area are acquired, wherein the multi-source remote sensing data comprise multi-temporal optical satellite images, multi-temporal SAR satellite images and auxiliary data, and the auxiliary data comprise a digital elevation model, meteorological observation data and earthquake catalog data;

[0009] The multi-temporal optical satellite images are subjected to radiation calibration, atmospheric correction and terrain correction, the multi-temporal SAR satellite images are subjected to precise orbit correction, thermal noise removal, radiation calibration and adaptive filtering, and all auxiliary data are uniformly resampled to a consistent spatio-temporal reference to obtain a standardized multi-source remote sensing data set consistent in space and time;

[0010] The multi-temporal optical satellite images in the multi-source remote sensing data set are processed by an improved frequency domain cross-correlation algorithm, including image precise registration, pixel offset calculation, adaptive window adjustment and system error correction, to obtain a horizontal displacement field of the research area;

[0011] Based on the multi-temporal SAR satellite images in the multi-source remote sensing data set and the digital elevation model, an optimized interference pair is selected, atmospheric phase delay and orbit error correction are performed, geometric distortion is identified and processed, and time-series deformation robustness is solved to obtain a time-series deformation field of the research area;

[0012] Based on the horizontal displacement field and the time-series deformation field, spatial adaptive grid fusion, time interpolation based on Gaussian process, and dynamic allocation of fusion weights of optical and SAR data according to deformation level, terrain conditions and data quality are performed in a unified spatio-temporal framework to obtain a spatio-temporally continuous comprehensive deformation field;

[0013] Based on the comprehensive deformation field, feature extraction and dimension reduction are performed by principal component analysis, and deformation patterns are automatically classified by using an unsupervised K-means clustering algorithm, and an activity level is defined by combining the category conversion of pre-seismic and post-seismic deformation patterns to obtain a deformation type map and an activity classification map of each observation point in the research area.

[0014] Further, when the multi-temporal optical satellite images are terrain corrected, a C correction model is adopted, and an expression of the C correction model is as follows:

[0015] ;

[0016] wherein, represents a radiation brightness value after terrain correction, represents an original radiation brightness value, represents a solar zenith angle, represents an incident angle, represents a correction parameter determined by statistical regression on a terrain flat area;

[0017] When the multi-temporal SAR satellite images are adaptively filtered, a filter window size is dynamically calculated and adjusted according to a local coherence coefficient, and a calculation formula is as follows:

[0018] ;

[0019] wherein, represents the filter window size, and respectively represent a preset minimum window size and a maximum window size, represents a local coherence coefficient.

[0020] Further, a processing flow of the pixel offset calculation includes:

[0021] A two-dimensional fast Fourier transform is performed on corresponding windows of the master image and the slave image to obtain a frequency spectrum of the master image and a frequency spectrum of the slave image, a normalized cross power spectrum is calculated according to the frequency spectrum of the master image and the frequency spectrum of the slave image, an inverse Fourier transform is performed on the normalized cross power spectrum to obtain a correlation plane in a spatial domain, and a pixel-level displacement is determined by finding a peak value of the correlation plane;

[0022] In a neighborhood of a peak value position of the pixel-level displacement, a two-dimensional Gaussian surface model is used to perform a least square fitting on discrete correlation plane intensity values, and a sub-pixel peak value position is solved as a final pixel offset;

[0023] A calculation formula of the normalized cross power spectrum is as follows:

[0024] ;

[0025] wherein, represents the normalized cross power spectrum, and represent frequency domain coordinate variables, represents a frequency spectrum of the master image, represents a frequency spectrum of the slave 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 represent a neighborhood block centered at pixel and represent a neighborhood block centered at pixel represents a filter strength parameter, represents a normalization factor, represents a pixel offset of the i-th pixel before filtering.

[0037] Further, the processing procedure of the interference pair preferably comprises:

[0038] Based on the multi-temporal SAR satellite images, for all possible combinations of interference pairs, the weight of each interference pair is calculated based on its spatial baseline, time baseline and predicted coherence, and an optimized interference pair set is constructed according to the weight for time series deformation robustness solving;

[0039] The weight calculation formula of the interference pair is as follows:

[0040] ;

[0041] wherein, represents the weight of the interference pair, represents the predicted coherence, represents the vertical spatial baseline, represents the critical baseline, represents the time baseline of the interference pair, represents a feature decorrelation time constant, represents a natural exponential function;

[0042] The processing procedure of the atmospheric phase delay comprises:

[0043] Based on the digital elevation model and the interference phase map, the regression coefficient is estimated by the phase-elevation statistical relationship of the stable area, and the terrain-related atmospheric delay signal is calculated and removed by using the elevation correlation model; the elevation correlation model is as follows:

[0044] ;

[0045] wherein, represents the atmospheric delay, represents the terrain elevation, represents the regression coefficient, represents the residual turbulent atmosphere phase component;

[0046] The processing procedure of the geometric distortion identification comprises:

[0047] Based on the digital elevation model and the SAR satellite orbit parameters, the ​exponent value, when determined as a shadow area, when determined as an overlay area, and the identified shadow area and overlay area are shielded;

[0048] The exponent value is calculated as follows:

[0049] ;

[0050] wherein, represents the exponent value, represents the terrain slope, represents the slope direction, represents the radar incidence angle, represents the satellite azimuth angle.

[0051] Further, the processing flow of the time-series deformation robustness solution includes:

[0052] Based on the optimized interference pair set, the observation equation is constructed: ;

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

[0054] wherein, represents the unwrapped phase vector, represents the design matrix, represents the deformation rate vector, represents the observation noise, represents the weight matrix, represents the smoothing constraint matrix, represents the regularization parameter, represents the matrix transpose.

[0055] Further, the processing flow of the spatial adaptive grid fusion includes:

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

[0057] ;

[0058] wherein, represents the grid size, represents the grid size corresponding to the evaluation function value, represents the number of grids containing observation points, represents the total number of grids, denotes a characteristic scale, denotes a weight coefficient, denotes a natural exponential function;

[0059] After determining the optimal grid size, the optical and InSAR deformation variables are fused in each grid cell using distance inverse weight interpolation:

[0060] ;

[0061] ;

[0062] wherein, denotes a fused deformation variable, denotes a weight of an th observation point, denotes a deformation variable of an th observation point, denotes an observation accuracy of an th observation point, denotes a distance from an th observation point to the center of the grid, denotes a smoothing parameter;

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

[0064] ;

[0065] wherein, denotes a deformation variable at time , denotes a mean function, denotes a Gaussian process, denotes a covariance function, and denotes different time points, the covariance function is a Matérn kernel:

[0066] ;

[0067] wherein, denotes a signal variance, denotes a characteristic time scale;

[0068] In the dynamic allocation of the fusion weights of the optical and SAR data according to the deformation variable level, the terrain condition and the data quality, the optical data weight and the SAR data weight are calculated and determined according to the following formula:

[0069] ;

[0070] ;

[0071] wherein, denotes the optical data weight, denotes the SAR data weight, denotes the deformation rate, denotes the deformation rate threshold, denotes the terrain slope, denotes the deformation rate feature parameter, denotes the slope feature parameter, and denote the quality indicators of optical and SAR data, respectively.

[0072] Further, feature extraction and dimension reduction are performed by principal component analysis, including:

[0073] The time-series deformation data corresponding to the deformation matrix of the integrated deformation field is normalized:

[0074] ;

[0075] wherein, denotes the deformation of the th observation point at the th time sampling point, denotes the normalized deformation, denotes the mean of all time sampling deformations of the th observation point, denotes the standard deviation of all time sampling deformations of the th observation point;

[0076] A covariance matrix is constructed, and eigenvalues and eigenvectors are obtained by eigen decomposition of the covariance matrix, and the covariance matrix is as follows:

[0077] ;

[0078] wherein, denotes the covariance matrix, denotes the deformation matrix corresponding to the normalized deformation, denotes the matrix transpose, denotes the number of observation points;

[0079] The deformation matrix corresponding to the normalized deformation is projected onto the eigenvectors to obtain principal components:

[0080] ;

[0081] wherein, denotes the th principal component, denotes the th eigenvector obtained by eigen decomposition;

[0082] The cumulative variance contribution rate of the characteristic value is calculated, and the minimum value that makes the cumulative variance contribution rate reach 85% or more is selected as the principal component number, that is:

[0083] ;

[0084] Wherein, represents the first characteristic value obtained by characteristic decomposition, represents the principal component number, represents the total number of characteristic values;

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

[0086] Set the range of candidate cluster number, for each candidate number, execute K-means clustering, take the deformation variable of each observation point as a sample point, calculate the silhouette coefficient of all sample points and take the average, select the candidate number that makes the average silhouette coefficient maximum as the final cluster number, the calculation formula of the silhouette coefficient is as follows:

[0087] ;

[0088] Wherein, represents the silhouette coefficient of the i-th sample point, represents the average distance from the i-th sample point to other sample points in the same cluster, represents the average distance from the i-th sample point to all sample points in the nearest neighbor cluster.

[0089] When initializing unsupervised K-means clustering, the selection probability of each sample point is calculated, and the initial cluster center is selected according to the probability, the calculation formula of the probability is as follows:

[0090] ;

[0091] Wherein, represents the probability of the i-th sample point being selected as the next cluster center, represents the Euclidean distance from the sample point to the nearest cluster center, represents the Euclidean distance from the i-th sample point to the nearest cluster center, represents 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. represents the cumulative precipitation within the sliding time window, represents the precipitation response time lag, represents the precipitation influence coefficient;

[0108] The temperature trigger effect modeling comprises: calculating a positive accumulated temperature index for quantifying the intensity of freeze-thaw cycles:

[0109]

[0110] wherein, represents the positive accumulated temperature index, represents the daily average temperature, represents the thawing threshold, represents the summation period;

[0111] The earthquake trigger effect evaluation comprises: using a modified Newmark displacement model, in combination with the peak acceleration of ground motion, the duration and the initial stability condition of the high-risk area, to evaluate the cumulative impact of the earthquake event on the stability of the loose body, and to generate an earthquake influence factor;

[0112] The flow of the time lag correlation analysis comprises:

[0113] Based on the integrated deformation field and the meteorological observation data, the cross-correlation coefficient of the deformation time series and the trigger factor time series under different precipitation response time lags is calculated, and the optimal precipitation response time lag is determined by finding the maximum cross-correlation coefficient; the calculation formula of the cross-correlation coefficient is as follows:

[0114]

[0115] wherein, represents the cross-correlation coefficient, represents the deformation time series, represents the trigger factor time series, and respectively represent the corresponding mean values;

[0116] A multi-factor comprehensive trigger index is constructed to normalize and weightedly fuse the influences of different trigger factors under the optimal precipitation response time lag:

[0117]

[0118] wherein, represents the comprehensive trigger index, and respectively represent the normalized precipitation index and the temperature index, represents the earthquake influence factor, respectively represent the corresponding weights;

[0119] ​​​The flow of the critical state recognition comprises: based on the comprehensive deformation field, calculating the reciprocal of the deformation rate time series data of the target region in a sliding time window, and performing linear fitting on the reciprocal by using a Fukuzono-Voight model, when the slope of the straight line obtained by fitting is negative and the goodness of fitting is greater than a fitting threshold, it indicates that the deformation is accelerating and the system tends to be unstable; the Fukuzono-Voight model is as follows:

[0120] ;

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

[0122] The beneficial effects of the present application are: the high-position loose body multi-scale deformation monitoring method based on optical and SAR provided by the present application firstly combines the advantages of optical image capturing millimeter-level large deformation and SAR image detecting millimeter-level micro deformation, realizes continuous high-precision monitoring of the whole life cycle of high-position loose body from initial creep to accelerated instability, and significantly improves the integrity and accuracy of deformation information in the time and space dimensions. Secondly, the present application changes the traditional inefficient and subjective mode which relies on manual interpretation by using the automatic deformation classification technology based on machine learning, realizes objective, quantitative and efficient classification and recognition of massive time series deformation data, provides reliable and comprehensive data cornerstone for subsequent construction of accurate and advanced geological disaster warning model, and improves the prevention and control ability and decision level of high-position loose body disaster as a whole. BRIEF DESCRIPTION OF DRAWINGS

[0123] Figure 1 The flowchart of the high-position loose body multi-scale deformation monitoring method based on optical and SAR provided by the embodiment is shown;

[0124] Figure 2 The deformation rate comparison and acceleration process diagram of different regions before and after the earthquake provided by the embodiment is shown;

[0125] Figure 3 The grid resolution optimization curve diagram provided by the embodiment is shown;

[0126] Figure 4 The deformation rate and weight distribution scatter diagram provided by the embodiment is shown;

[0127] Figure 5 The data coverage comparison diagram provided by the embodiment is shown;

[0128] Figure 6 The activity classification result diagram provided by the embodiment is shown;

[0129] Figure 7 A schematic diagram of the temporal deformation evolution of the disaster event in October 2018 is provided for the embodiment;

[0130] Figure 8 A schematic diagram of the temporal precipitation and temperature evolution of the disaster event in October 2018 is provided for the embodiment;

[0131] Figure 9 A schematic diagram of the multi-factor trigger analysis of the disaster event in October 2018 is provided for the embodiment. DETAILED DESCRIPTION

[0132] In order to overcome the limitations of a single remote sensing data source, realize continuous and automatic deformation monitoring and activity evaluation of high loose bodies from initial creep to rapid instability, the technical scheme of the present application is proposed. In the present application, first, according to the different physical characteristics and error sources of optical and SAR data, precise correction is carried out respectively, and they are unified to a consistent space-time reference; then, using an improved frequency domain cross-correlation algorithm, large-scale (meter level) deformation in the horizontal direction is extracted from the optical image sequence with high precision, and using an optimized time series interferometric measurement technology, small deformation in the radar line-of-sight direction is precisely detected from the SAR image sequence, and the two types of deformation fields obtained above, which are different in spatial distribution, time sampling and magnitude, are fused into a spatio-temporal continuous comprehensive deformation field through adaptive gridding and time interpolation based on Gaussian process; then, the fused comprehensive deformation field is regarded as a high-dimensional data set, and principal component analysis is used for dimension reduction to extract the core deformation features, then an unsupervised clustering algorithm is used to automatically classify points with similar deformation evolution behavior into the same category, thereby realizing objective quantitative classification of deformation modes. Finally, by analyzing the conversion of clustering categories before and after the earthquake, the high-risk areas activated by the earthquake are automatically identified, and intelligent evaluation of loose body activity is completed.

[0133] The technical scheme in the present embodiment will be described clearly and completely in combination with the drawings in the present embodiment. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments.

[0134] Figure 1 A flowchart of a multi-scale deformation monitoring method for high loose bodies based on optical and SAR is shown, please refer to Figure 1 The method comprises the following steps:

[0135] Step 1, multi-source remote sensing data acquisition and preprocessing:

[0136] Acquiring multi-source remote sensing data of the research area, the multi-source remote sensing data comprising multi-temporal optical satellite images, multi-temporal SAR satellite images and auxiliary data, the auxiliary data comprising a digital elevation model, 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, and all auxiliary data are uniformly resampled to a consistent space-time reference to obtain a standardized multi-source remote sensing data set that is consistent in space and time.

[0138] The core of this step is to establish a standardized multi-source remote sensing data acquisition and preprocessing process to form a high-quality basic data set to support subsequent deformation detection. By formulating a data selection strategy that conforms to the spatiotemporal evolution characteristics of high loose bodies and implementing systematic preprocessing, the radiometric, geometric and temporal differences between different sensors are eliminated, and the consistency of multi-source remote sensing data is achieved.

[0139] In terms of optical satellite images, images with a spatial resolution better than 10 meters and a revisit period shorter 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 performed during selection: cloud coverage must be less than 30%, and the solar elevation angle must be greater than 30°, to ensure image quality and minimize shadow interference. The selection of time phases is mainly based on seasonal transition periods, such as late spring and early summer, and late autumn and early winter, when loose body activity is significant and deformation features are more easily identified. Preprocessing includes three key steps: radiometric calibration, atmospheric correction and terrain correction. The improved 6S radiation transfer model is used for atmospheric correction to improve correction accuracy by introducing local measured atmospheric parameters; terrain correction is based on the C correction model to eliminate radiation distortion caused by terrain undulations, and its expression is:

[0140] ;

[0141] wherein, represents the radiation brightness value after terrain correction, represents the original radiation brightness value, represents the solar zenith angle, represents the incident angle, represents the correction parameter determined by statistical regression on flat terrain.

[0142] The SAR satellite image processing adopts the interferometric wide swath mode product of Sentinel-1 satellite, and the ascending and descending data are synchronously acquired to optimize the observation geometry. The data acquisition time span is not less than two years, which ensures to cover the complete annual cycle and the period before and after the key events. The preprocessing process includes precise orbit correction, thermal noise removal and radiation calibration. In order to enhance the interferometric coherence in complex terrain areas, an adaptive filtering algorithm is introduced in this embodiment, which dynamically calculates and adjusts the filter window size according to the local coherence coefficient, and the calculation formula is as follows:

[0143] ;

[0144] wherein, represents the filter window size, and respectively represent the preset minimum window size and maximum window size, represents the local coherence coefficient.

[0145] The integration of auxiliary data is an important link to ensure the monitoring accuracy. The system collects multiple types of auxiliary data such as high-precision digital elevation model, meteorological observation data and earthquake catalog. The DEM data takes TanDEM-X or ALOS PALSAR product with a resolution higher than 12 meters as the first choice, and the missing area is supplemented by 30-meter SRTM data; the meteorological data includes daily average temperature, cumulative precipitation and snow depth, etc., which are interpolated to the global domain by means of distance inverse weight method; the earthquake data selects events with magnitude greater than 3.0 and focal depth less than 30 kilometers, which are used for subsequent analysis of earthquake triggering effect. All auxiliary data are resampled to 30-meter spatial resolution and aligned with the remote sensing observation time, thereby establishing a time and space consistent multi-source remote sensing data fusion basis.

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

[0147] The multi-temporal optical satellite images in the multi-source remote sensing data set are processed by the improved frequency domain cross-correlation algorithm, including image precise 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 realize high-precision extraction of sub-pixel level displacement field in optical images by improved frequency domain cross-correlation algorithm, and is especially optimized for large-scale rapid deformation of high-bit loose body. The traditional COSI-Corr algorithm is easily disturbed by registration error and noise in complex mountainous terrain, and this embodiment introduces adaptive window strategy, multi-scale correlation analysis and robust error correction mechanism, which significantly improves the precision and reliability of deformation detection.

[0149] Image registration is the prerequisite for deformation detection. In this paper, a hierarchical registration strategy is adopted. Firstly, the initial registration parameters are obtained by matching SIFT feature points, and then the local area is optimized. The density of feature points is adjusted adaptively according to the complexity of the terrain: 20-30 feature points per square kilometer in flat areas, and 50-80 in steep mountainous areas. Rational function model (RFM) is selected for registration transformation, and the parameters are solved by least squares adjustment. In order to exclude false matching, RANSAC algorithm is introduced to remove gross errors, and the iteration threshold is set to 1.5 times the mean error. Finally, the root mean square error of the check points is used to evaluate the registration accuracy, which is required to be 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 in a sliding window. The processing flow of pixel offset calculation includes:

[0151] The two-dimensional fast Fourier transform is performed on the corresponding windows of the master image and the slave image to obtain the frequency spectrum of the master image and the slave image. The normalized cross-power spectrum is calculated according to the frequency spectrum of the master image and the slave image. The inverse Fourier transform is performed on the normalized cross-power spectrum to obtain the correlation surface in the spatial domain. The pixel-level displacement is determined by finding the peak value of the correlation surface. In the neighborhood of the peak value position determined by the pixel-level displacement, the least squares fitting is performed on the discrete correlation surface intensity values by using the two-dimensional Gaussian surface model to solve the sub-pixel peak position as the final pixel offset.

[0152] The calculation formula of the normalized cross-power spectrum is as follows:

[0153] ;

[0154] Wherein, represents the normalized cross-power spectrum, and represents the frequency domain coordinate variable, represents the frequency spectrum of the master image, represents the frequency spectrum of the slave image, represents the complex conjugate of represents a small amount to prevent division by zero;

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

[0156] ;

[0157] Wherein, represents the correlation surface intensity value, i.e. the intensity value of the continuous correlation surface fitted by the two-dimensional Gaussian function in the local neighborhood centered on the integer-level correlation peak point, which is a continuous function for simulating and approximating the distribution of the real discrete cross-correlation value, and denotes a spatial domain coordinate variable, denotes a sub-pixel peak position, denotes a peak amplitude, denotes a Gaussian distribution, denotes a standard deviation of the Gaussian distribution in denotes a standard deviation of the Gaussian distribution in denotes a background noise level, denotes a natural exponential function.

[0158] In this embodiment, the adaptive window adjustment strategy can dynamically optimize the relevant window size according to the local texture richness and the deformation gradient. The texture richness is evaluated by calculating the local standard deviation, and the deformation gradient is estimated by the difference of the displacement of adjacent windows. The processing flow of adaptive window adjustment includes:

[0159] In the pixel offset calculation process, the local standard deviation of each potential window position is calculated in real time, and the deformation gradient module thereof is estimated. The window size for the intercorrelation calculation at this position is dynamically calculated and determined according to the local standard deviation and the deformation gradient module, and the calculation formula is as follows:

[0160] ;

[0161] wherein, denotes the window size after adaptive adjustment, denotes the reference window size, denotes the local standard deviation, denotes the minimum standard deviation threshold, denotes the deformation gradient module, denotes the maximum expected deformation, denotes a natural exponential function. denotes a weight coefficient;

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

[0163] In this embodiment, a three-level correction strategy is adopted: the first level is track error correction, which fits the residual displacement surface of the stable area through a robust estimation method, and automatically removes outliers exceeding 3 times the standard deviation; the second level is for strip noise, which uses a directional notch filter to suppress specific frequency components, and the bandwidth is adaptively adjusted according to the strip period; the third level is random noise suppression, which uses a non-local mean filter algorithm to effectively reduce noise while maintaining clear deformation boundaries:

[0164] ;

[0165] wherein, denotes the pixel offset of the i th pixel point after filtering, denotes the pixel offset of the i th pixel point before filtering.​ , and represent a neighborhood block centered at pixel and , represents a filter strength parameter, represents a normalization factor, represents a pixel offset of the pixel before filtering.

[0166] Through the above multi-stage correction, the system error of deformation detection is reduced to below 0.1 pixel, and the random error is reduced by more than 60%.

[0167] Step 3, deformation extraction of SAR path:

[0168] Based on the multi-temporal SAR satellite images and digital elevation models in the multi-source remote sensing data set, the time series deformation field of the study area is obtained through optimized interference pair optimization, atmospheric phase delay and orbit error correction, geometric distortion identification and processing, and time series deformation robustness solution.

[0169] This step focuses on the detection of early creep and micro-deformation of high-level loose bodies. Through the optimized small baseline set interferometric measurement technology, millimeter-level precision time series deformation monitoring is realized. Traditional SBAS method is easily disturbed by decorrelation, atmospheric delay and geometric distortion in high mountain complex terrain area, which limits the monitoring accuracy and coverage. Therefore, the embodiment improves the system from multiple aspects such as interference pair optimization, error correction and geometric distortion processing, which significantly improves the extraction ability and reliability of small deformation signals.

[0170] Among them, the optimization of interference pairs is the basis of time series analysis quality. This embodiment proposes an intelligent screening mechanism based on space-time baseline and coherence prediction: the time baseline is controlled within 180 days to maintain the continuous capture ability of seasonal deformation; the spatial baseline is constrained within 40% of the critical baseline to maintain high interference coherence. The key innovation is to introduce a coherence prediction model, which comprehensively considers the factors of land cover type, terrain slope and seasonality to predict the expected coherence of each interference pair.

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

[0172] The weight calculation formula of the interference pair is as follows:

[0173] ;

[0174] wherein, denotes the weight of an interferometric pair, denotes the predicted coherence, denotes the vertical spatial baseline, denotes the critical baseline, denotes the temporal baseline of an interferometric pair, denotes the characteristic decorrelation time constant, denotes the natural exponential function.

[0175] Through such a weighted selection, the number of interferometric pairs is controlled between 3N~5N, both ensuring the redundancy and avoiding excessive computational burden.

[0176] The precise correction of atmospheric phase delay and orbit error is crucial for deformation extraction. This embodiment adopts a multi-scale separation strategy to gradually eliminate various errors: first, the global atmospheric model GACOS is used to remove the large-scale tropospheric delay, and then the turbulence atmospheric signal and the true deformation are separated by spatial-temporal filtering. In view of the significant influence of terrain on atmospheric delay, the processing flow of atmospheric phase delay includes:

[0177] Based on the digital elevation model and the interferometric phase map, the regression coefficient is estimated through the phase-elevation statistical relationship of the stable area, and the terrain-related atmospheric delay signal is calculated and removed by using the elevation correlation model; the elevation correlation model is as follows:

[0178] ;

[0179] wherein, denotes the atmospheric delay, denotes the terrain elevation, denotes the regression coefficient, denotes the residual turbulence atmospheric phase component.

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

[0181] ;

[0182] wherein, denotes the index value, denotes the terrain slope, denotes the aspect, denotes the radar incidence angle, denotes the satellite azimuth.

[0183] When , it is determined as a shadow area, and when , it is determined as a layover area.The shadow area and the overlap area are shielded. In order to improve the monitoring coverage and give full play to the complementarity of the ascending and descending track data, the data of the other track is used to supplement the data of the single track in the area affected by the distortion.

[0184] The robustness of the time sequence deformation in the embodiment is solved by using an improved singular value decomposition method to solve the deformation time sequence. The processing flow includes:

[0185] Based on the optimized interference pair set, the observation equation is constructed: ;

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

[0187] wherein, represents the unwrapped phase vector, represents the design matrix, represents the deformation rate vector, represents the observation noise, represents the weight matrix, represents the smoothing constraint matrix, represents the regularization parameter, represents the matrix transpose.

[0188] After 35 iterations of calculation, the abnormal observations can be automatically identified and the weight thereof can be reduced, and finally the estimation accuracy of the annual average deformation rate is better than 2 mm. In order to facilitate geological interpretation, the radar line-of-sight deformation is finally converted to the along-slope direction, supporting subsequent stability analysis and disaster warning.

[0189] Step 4, spatio-temporal fusion of multi-source deformation results:

[0190] Based on the horizontal displacement field and the time sequence deformation field, spatial adaptive grid fusion, time interpolation based on Gaussian process, and dynamic allocation of fusion weights of optical and SAR data according to the deformation level, terrain conditions and data quality are performed in a unified spatio-temporal framework to obtain a spatio-temporally continuous comprehensive deformation field.

[0191] The core goal of this step is to deeply fuse the large deformation information obtained by tracking optical offset and the small deformation signal detected by InSAR in a unified spatio-temporal framework to overcome the differences in spatial coverage, time sampling and deformation level between the two types of technologies. 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, which provides a reliable data basis for subsequent activity identification.

[0192] In the aspect of spatial adaptive grid fusion, an adaptive grid strategy is adopted to dynamically determine the optimal grid size according to the spatial distribution density of deformation points and deformation gradient characteristics. The traditional fixed grid method is difficult to balance detail preservation and noise suppression, and the present embodiment identifies the inflection point position of grid refinement through the maximum curvature point method.

[0193] In the present embodiment, the processing flow of spatial adaptive grid fusion 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, the balance between data coverage and spatial resolution is achieved, and the corresponding evaluation function is as follows:

[0195] ;

[0196] Among them, represents the grid size, represents the grid size corresponding to the evaluation function value, represents the number of grids containing observation points, represents the total number of grids, represents the characteristic scale, represents the weight coefficient, represents the natural exponential function;

[0197] After determining the optimal grid size, in each grid unit, the optical and InSAR deformation variables are fused by using the distance inverse weight interpolation:

[0198] ;

[0199] ;

[0200] Among them, represents the fused deformation variable, represents the weight of the th observation point, represents the deformation variable of the th observation point, represents the observation accuracy of the th observation point, represents the distance from the th observation point to the grid center, represents the smoothing parameter.

[0201] This weighting scheme not only considers the difference in observation accuracy, but also ensures spatial continuity.

[0202] In the aspect of time interpolation based on Gaussian process, a time interpolation model based on Gaussian process is established to unify irregular sampling to standard time series in view of the inconsistency between uneven sampling caused by cloud interference of optical data and regular revisit of SAR data. The deformation time evolution can be expressed as:

[0203] ;

[0204] where, denotes the deformation at time , denotes the mean function, denotes the Gaussian process, denotes the covariance function, and denotes different time points, the covariance function is Matérn kernel:

[0205] ;

[0206] where, denotes the signal variance, denotes the characteristic time scale.

[0207] The hyperparameters are determined by maximum likelihood estimation to realize smooth interpolation between observation points and uncertainty quantification.

[0208] In the aspect of dynamically assigning fusion weights of optical and SAR data according to deformation level, terrain condition and data quality, the weights are dynamically adjusted according to technical application conditions: when the deformation rate exceeds 100 mm / year, optical data is preferred, and when it is lower than 50 mm / year, InSAR is mainly used; when the terrain slope is greater than 35°, the weight of InSAR is reduced, and in the dense vegetation area, the weight of optical data is reduced. The optical data weight and SAR data weight are calculated and determined according to the following formula:

[0209] ;

[0210] ;

[0211] where, denotes the optical data weight, denotes the SAR data weight, denotes the deformation rate, denotes the deformation rate threshold, denotes the terrain slope, denotes the deformation rate characteristic parameter, denotes the slope characteristic parameter, and denote the quality indicators of optical and SAR data respectively.

[0212] ​In practical application, the quality of the fusion result can also be evaluated by cross-validation and uncertainty analysis. 10% of the observation points are randomly reserved as validation samples, and the root mean square error between the fusion value and the measured value is required to be better than 1.5 times the accuracy of any single data source. At the same time, based on the error propagation theory, the uncertainty of the fusion result is estimated, and the confidence layer is generated. The area with uncertainty exceeding 50% of the deformation is marked as low confidence for cautious use in subsequent analysis. After the above processing, the final comprehensive deformation field with complete spatial coverage and continuous time is formed, the spatial coverage rate is improved by more than 40%, and the time resolution is 12 days, which provides a solid data foundation for the accurate identification of active high loose bodies.

[0213] Step 5, automatic recognition of deformation features and activity classification based on machine learning:

[0214] Based on the comprehensive deformation field, feature extraction and dimension reduction are performed by principal component analysis, and unsupervised K-means clustering algorithm is used for automatic classification of deformation mode. Combined with the category conversion of pre-seismic and post-seismic deformation mode, the activity level is defined, and the deformation type map and activity classification map of each observation point in the study area are obtained.

[0215] This step uses unsupervised machine learning method to automatically classify and identify time series deformation data, aiming to overcome the strong subjectivity and low efficiency in traditional manual interpretation. The core is to project high-dimensional time series deformation data into low-dimensional feature space, automatically identify different deformation modes by clustering analysis, and define the activity level according to the transformation of pre-seismic and post-seismic deformation mode, so as to provide objective and quantitative basis for risk assessment of high loose bodies.

[0216] First, feature extraction and dimension reduction are performed by principal component analysis:

[0217] The time series deformation data corresponding to the deformation matrix of the comprehensive deformation field are standardized:

[0218] ;

[0219] Among them, represents the deformation of the th observation point at the th time sampling point, represents the standardized deformation, represents the mean of all time sampling deformation of the th observation point, represents the standard deviation of all time sampling deformation of the th observation point;

[0220] The covariance matrix is constructed, the covariance matrix is decomposed to get the eigenvalue and eigenvector, and the covariance matrix is as follows:

[0221] ;

[0222] wherein, denotes the covariance matrix, denotes the deformation matrix corresponding to the normalized deformation variable, denotes the matrix transpose, denotes the number of observation points;

[0223] Project the deformation matrix corresponding to the normalized deformation variable onto the eigenvectors to obtain the principal components:

[0224] ;

[0225] wherein, denotes the th principal component, denotes the th eigenvector obtained by eigen decomposition;

[0226] Calculate the cumulative variance contribution rate of the eigenvalues, and select the minimum value that makes the cumulative variance contribution rate reach more than 85% as the number of principal components, i.e.:

[0227] ;

[0228] wherein, denotes the th eigenvalue obtained by eigen decomposition, denotes the number of principal components, denotes 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 embodies seasonal changes, and the third principal component represents sudden changes.

[0230] Then, based on the reduced feature space, the K-means clustering algorithm is used to automatically identify deformation types. The determination of the number of clusters is a key parameter, and the embodiment automatically optimizes it through the Silhouette Coefficient: set a range of candidate cluster numbers, perform K-means clustering for each candidate number, take the deformation variable of each observation point as a sample point, calculate the Silhouette Coefficient of all sample points and take the average, and select the candidate number that makes the average Silhouette Coefficient the largest as the final number of clusters. The calculation formula of the Silhouette Coefficient is as follows:

[0231] ;

[0232] wherein, denotes the 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] Further, by analyzing the conversion of deformation patterns before and after the earthquake, the 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-stable area; type C (active before earthquake-active after earthquake) is a continuous active area; and type D (stable before earthquake-active after earthquake) is a seismic activation area. Among them, the D-type area has important disaster prevention significance. When the deformation acceleration exceeds 3 times the standard deviation of the background value, the duration is more than 30 days, and the influence range is greater than 1 square kilometer, it is determined as a high-risk active high loose body. The data-driven automatic identification method has improved the processing efficiency by more than 20 times compared with manual interpretation, and the identification accuracy is 92%.

[0242] In the embodiment, further comprising:

[0243] Step 6, construction of a high loose body early identification model and generation of early warning information:

[0244] Based on the activity classification map, meteorological observation data and comprehensive deformation field, quantitative analysis of multi-factor triggering mechanism and time lag correlation analysis are carried out, critical state recognition is carried out based on Fukuzono-Voight model, and early warning information is generated according to preset classification early warning criteria.

[0245] This step aims to build an early identification model that comprehensively considers deformation evolution law, external triggering factors and environmental conditions, to realize quantitative evaluation and early warning of high loose body from stable state to instability. The core innovation of the model is to establish a quantitative response relationship between deformation acceleration and external triggering factors, to reveal the climate driving mechanism through time lag correlation analysis, and to determine the classification early warning threshold based on the critical state theory, so as to advance the early warning time by 10-20 days, and provide scientific basis for disaster prevention and control decision-making.

[0246] Quantitative analysis of multi-factor triggering mechanism is the theoretical basis of early identification. The embodiment considers four main triggering factors of precipitation, temperature, earthquake and human activity, and constructs a deformation response model. For climate factors, the cross-correlation function of deformation rate and meteorological elements is calculated by using a sliding time window to identify the best response time lag.

[0247] Among them, the triggering effect of precipitation is described by the relationship between cumulative precipitation and deformation acceleration: based on the sliding time window, a quantitative relationship between deformation acceleration and cumulative precipitation is established, and the model is as follows:

[0248] ;

[0249] Among them, represents the deformation acceleration at time t, represents the background acceleration, represents the background acceleration, represents the cumulative precipitation within the sliding time window, represents the precipitation response time lag, represents the precipitation influence coefficient.

[0250] The temperature influence is mainly reflected in the weakening of the loose body strength by freeze-thaw cycles. The positive accumulated temperature index is used to quantify the strength of freeze-thaw cycles, and the calculation formula is as follows:

[0251] ;

[0252] wherein, represents the positive accumulated temperature index, represents the daily average temperature, represents the melting threshold (0℃), represents the summation period, which is 30-60 days.

[0253] The earthquake triggering effect evaluation includes: using the modified Newmark displacement model, combining the peak ground acceleration, duration and initial stability condition of high-risk areas, evaluating the cumulative influence of earthquake events on the stability of loose bodies, and generating an earthquake influence factor.

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

[0255] Based on the comprehensive deformation field and meteorological observation data, the cross-correlation coefficient of deformation time series and trigger factor time series under different precipitation response time lags is calculated, and the best precipitation response time lag is determined by finding the maximum cross-correlation coefficient; the calculation formula of the cross-correlation coefficient is as follows:

[0256] ;

[0257] wherein, represents the cross-correlation coefficient, represents the deformation time series, represents the trigger factor time series, and respectively represent the corresponding mean value;

[0258] A multi-factor comprehensive trigger index is constructed to normalize and weightedly fuse the influence of different trigger factors under the best precipitation response time lag:

[0259] ;

[0260] wherein, represents the comprehensive trigger index, and respectively represent the normalized precipitation index and temperature index, represents the earthquake influence 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 a danger level, an estimated occurrence time, an influence range, and a disposal suggestion. Through historical case backtracking verification, the model hit rate reaches 85%, the false alarm rate is lower than 20%, the average early warning advance period is 15 days, and the prevention and control capability of the high loose body disaster is significantly improved.

[0274] In summary, the high loose body multi-scale deformation monitoring method based on optical and SAR provided by the embodiment, by constructing a complete technical chain of "data acquisition-feature extraction-fusion analysis-advance identification", by giving full play to the complementary advantages of optical and SAR remote sensing data, realizing the quantitative monitoring and early warning of the whole process from stability to instability of loose body in high mountain area, not only effectively solves the key problems of data source fusion, feature recognition and early warning timeliness of the prior art, but also provides a new technical path for the systematic prevention and control of geological disasters in high mountain area. Specifically, the embodiment adopts a hierarchical progressive processing strategy. In the data acquisition stage, a standardized processing flow of multi-source remote sensing data is established to ensure the comparability and consistency of different sensors and different temporal data. Then, the improved pixel offset tracking algorithm is used to capture large-scale rapid deformation, and the optimized time series InSAR technology is used to detect early micro-deformation signals. The two technical paths are advanced in parallel to obtain complementary deformation information. On this basis, the adaptive grid fusion algorithm is used to integrate the heterogeneous deformation results into a unified space-time reference system, and then the machine learning algorithm is used to realize the automatic classification of deformation features, and finally the activity advance identification model is established combined with external trigger factor analysis. From the perspective of system architecture, the embodiment constructs a modular processing framework, covering five core components of data input, deformation detection, data fusion, feature recognition and early warning output. These modules realize information transmission through standardized data interfaces, not only ensuring the scalability of the system, but also providing flexible technical support for practical application.

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

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

[0277] Compared with the monitoring blind area problem caused by the single remote sensing data source in the prior art, the embodiment realizes the synchronous monitoring of large deformation (meter level) and small deformation (millimeter level) by deeply fusing optical and SAR remote sensing data, effectively solving the technical bottlenecks that the optical technology cannot capture early micro-deformation and the SAR technology is not correlated in the rapid deformation area. Through the adaptive grid fusion and dynamic weight distribution strategy, the spatial coverage of the deformation monitoring is improved from 60% in the prior art to more than 85%, and the time resolution is improved from the monthly scale to 12 days. The test results show that, under the complex mountainous terrain conditions, the embodiment can completely track the whole process evolution of the high-position loose body from the initial creep, accelerated deformation to the final instability, while the existing single technology method can only capture the deformation information of part of the stage, and the best early warning opportunity is easily missed.

[0278] (2) The automatic recognition technology greatly improves the efficiency and objectivity of deformation classification;

[0279] In view of the strong subjectivity and low efficiency problem caused by the dependence of the prior art on manual interpretation and fixed threshold, the automatic recognition framework based on machine learning constructed in the embodiment realizes the objective quantitative classification of deformation features. Through principal component analysis, the complex time series deformation data is reduced to 2-3 key features, and different deformation modes are automatically recognized by combining an unsupervised clustering algorithm, the processing efficiency is improved by more than 20 times compared with the traditional manual method, and the processing time of a single image is shortened from several hours to within 10 minutes. More importantly, the embodiment eliminates the uncertainty of human experience judgment, the consistency of repeated processing results of the same data set reaches more than 98%, while the results of different interpreters using the traditional method can differ by 30%, ensuring the repeatability and comparability of the monitoring results.

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

[0281] Compared with the defect that the prior art cannot effectively distinguish different triggering factors, the embodiment innovatively proposes a four-type activity classification method based on the pre-earthquake and post-earthquake deformation mode conversion, accurately identifies the D-type activation area induced by earthquakes, and realizes the 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, the embodiment can accurately delineate the high-risk area directly affected by the earthquake, with an identification accuracy of 92%, while the existing method can only give a vague qualitative judgment. This quantitative differentiation capability enables the disaster prevention department to develop differentiated prevention and control strategies for different triggering mechanisms, such as strengthening engineering management in earthquake-activated areas and strengthening monitoring and early warning in climate-sensitive areas, significantly improving the pertinence and effectiveness of disaster prevention and control.

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

[0283] The multi-factor coupling early identification model established in this embodiment reveals the lagging effect of precipitation (20±5 days) and temperature (10±3 days) on deformation through time delay correlation analysis, and in combination with the critical state identification criterion, the early warning period is extended from 3-5 days of the prior art to 10-20 days, which saves valuable time for emergency evacuation and engineering disposal. Compared with the 85% false alarm rate and passive situation of post-analysis of the prior art, the hit rate of the hierarchical early warning system of this embodiment reaches 85%, the false alarm rate is controlled within 20%, and the practicality requirement is basically reached. In particular, in the application case after the 2017 Nyingchi earthquake, this embodiment successfully warned 18 days in advance of a large-scale debris flow event in a certain basin in Southeast Tibet, while the traditional monitoring method found the anomaly only 2 days before the disaster occurred, fully proving the significant advantage of this embodiment in early identification.

[0284] The effectiveness of the method of the present application is verified below in a certain basin in Southeast Tibet as a typical research area:

[0285] (1) Research area overview and data acquisition:

[0286] The basin is located in the southeast margin of the Qinghai-Tibet Plateau, with an elevation range of 2800-6500 meters, a rugged terrain, developed glaciers, and widely distributed loose deposits. Since 1950, there have been 7 large-scale river blocking events in this area, among which the activity was significantly enhanced after the 2017 Nyingchi M6.9 earthquake on November 18, 2018, and two consecutive ice avalanche-debris flow disaster chains occurred on October 17 and 29, causing the Brahmaputra River to be blocked for dozens of hours, with direct economic losses exceeding 800 million yuan. Please refer to Table 1, according to the data acquisition standard of the present application, the system collects multi-source remote sensing data from October 2016 to February 2019, including 28 scenes of Planet satellite images (resolution 3.7 meters, covering the key period), 6 scenes of RapidEye images (resolution 5 meters, for verification), 15 scenes of Sentinel-2 images (resolution 10 meters, as the main optical data source), and 42 scenes of Sentinel-1 ascending track SAR data and 46 scenes of Sentinel-1 descending track SAR data, ensuring complete observation sequences before and after the earthquake, and strict data quality control is performed with cloud coverage less than 30% and coherence greater than 0.3.

[0287] Table 1 Multi-source remote sensing data acquisition and quality control statistics

[0288]

[0289] (2) Optical large deformation tracking implementation and verification

[0290] The optical image is processed by using the improved COSI-Corr algorithm, and the deformation characteristic changes before and after the earthquake in November 2017 are analyzed. Taking the Sentinel-2 image pair on November 11, 2017 and November 19, 2017 as an example, sub-pixel registration and offset calculation are carried out through an adaptive window strategy (initial 64*64 pixels, final 32*32 pixels, step 16 pixels), and east-west and north-south displacement fields are obtained. After removing the orbit error by applying a polynomial surface fitting, the effective pixel ratio is increased from the original 72% to 86%, and the noise level is reduced by 65%. As shown in Figure 2 , the maximum horizontal displacement of 15.3 meters / year is detected in the glacier area upstream of the basin, and the average moving rate of loose deposits in the main channel is 8.7 meters / year. In particular, in the channel II and IV regions, the post-earthquake accelerated material transport process is clearly identified, and the displacement rate of channel IV increases from 3.2 meters / year before the earthquake to 12.8 meters / year after the earthquake, which is highly consistent with the actual occurrence of the debris flow blocking river event on December 21, 2017 in space and time, verifying the prediction ability of the method.

[0291] (3) Time series InSAR small deformation detection implementation

[0292] The optimized SBAS-InSAR technology is used to process the pre-earthquake (October 2016 to November 17, 2017) and post-earthquake (November 19, 2017 to February 2019) Sentinel-1 data, and the R index model is used to identify and remove the geometric distortion area accounting for 23% of the study area, ensuring the reliability of the deformation results. In the pre-earthquake stage, 156 interference pairs are constructed, and in the post-earthquake stage, 198 interference pairs are constructed, with a time baseline controlled within 180 days and a vertical baseline limited within 40% of the critical baseline. After atmospheric correction and orbit error removal, a millimeter-level precision deformation rate field is obtained, with an average deformation rate of 12.5 mm / year before the earthquake and an increase to 35.8 mm / year after the earthquake. In the bedrock slope area in the middle and lower reaches of the basin, the creep signals of multiple new cracks are detected, with a deformation rate of 5-15 mm / year. These small deformations cannot be detected by optical methods, demonstrating the necessity of multi-source data fusion. Through joint analysis of ascending and descending track data, the LOS deformation is decomposed into three-dimensional deformation components, and it is found that the vertical subsidence is mainly concentrated in the moraine accumulation area, with a maximum subsidence rate of -85 mm / year.

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

[0294] Please refer to Figure 3 , through the adaptive grid fusion strategy of the application, the maximum curvature point method is used to determine the optimal grid size of 85 meters, which not only retains the local deformation details but also effectively suppresses the noise. Please refer toFigure 4 In the spatial fusion process, the weights of optical and InSAR data are dynamically adjusted according to the deformation gradient and data quality. In the rapid deformation area with a deformation rate greater than 100 mm / year, the weight of optical data reaches 0.7; in the slow deformation area with a deformation rate less than 50 mm / year, the InSAR weight is increased to 0.8. In the time dimension, Gaussian process interpolation is used to unify the irregular sampling to the standard time series with an interval of 12 days, and the interpolation uncertainty is controlled within 1.5 times of the original measurement error. Please refer to Figure 5 The spatial coverage of the fused deformation field is improved from 62% of single InSAR to 87%, and from 54% of single optical to 87%, achieving quasi-complete coverage of deformation information in the study area, with only a small amount of data blank in the extremely steep cliff and permanent shadow area.

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

[0296] The machine learning framework of the application is applied to automatically classify the fused time series deformation data. The first two principal components (cumulative variance contribution rate 87.6%) are extracted by principal component analysis. The first principal component reflects the long-term deformation trend, and the second principal component reflects the seasonal fluctuation, please refer to Table 2 and Table 3. Based on the contour coefficient optimization, the number of clusters is determined as 2, and 13,847 effective deformation points are automatically divided into two categories: stable type (type I, accounting for 68%) and active type (type II, accounting for 32%). Further according to the category conversion before and after the earthquake, four activity types are identified: A type stable area 8,234 points (59.5%), B type self-stable area 1,126 points (8.1%), C type sustained activity area 2,308 points (16.7%), and D type earthquake activated area 2,179 points (15.7%), please refer to Figure 6 The D type area is mainly distributed in the loose accumulation body within 60 kilometers from the epicenter, with an area of 4.2 square kilometers. These areas have experienced 3 large-scale instability events within 18 months after the earthquake, confirming the accuracy of the classification results. The entire automatic identification process takes only 45 minutes, while the traditional manual interpretation takes 3-5 working days.

[0297] Table 2 Deformation characteristic statistical table of different activity types

[0298]

[0299] Table 3 Deformation characteristic analysis table

[0300]

[0301] Note: ± indicates standard deviation; acceleration coefficient is the ratio of post-seismic and pre-seismic deformation rate; deformation direction dominance refers to the direction with a principal component contribution rate greater than 60%.

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

[0303] A multi-factor early identification model for the study area was constructed, and the best precipitation response time lag of meteorological factors was determined through sliding correlation analysis: 18 days for temperature and 22 days for precipitation. Please refer to Figure 7 to Figure 9 , Figure 8 The column chart shows the daily precipitation, and the dotted line shows the daily average temperature. Before the two disaster events in October 2018, the deformation D-shaped area experienced a high-temperature period of 15 consecutive days with daily average temperature above 10°C and a strong rainfall process with cumulative precipitation reaching 186 mm, triggering the index to exceed the warning threshold on October 5. The Fukuzono-Voight model was applied to analyze the critical state of the key points, and when the rate inverse approached zero and the goodness of fit R² was greater than 0.85, the system triggered an orange warning on October 12, predicting the instability time window to be October 15-20. The actual disaster occurred on October 17, with a 5-day early warning period and a 2-day time prediction error. Through backtracking tests on all deformation events from 2016 to 2019, the method successfully warned 9 out of 11 times (hit rate 81.8%), with 2 false alarms (false alarm rate 18.2%). The average early warning period was 12.5 days, which was significantly improved compared to the traditional monitoring method of 3 days. This provided sufficient time for the safe relocation of more than 3800 residents in the downstream 17 villages, fully verifying the practical value of the present embodiment in the early identification of active high loose bodies.

Claims

1. A method for high-relief loose material multiscale deformation monitoring based on optical and SAR, characterized in that, The method comprises: acquiring multi-source remote sensing data of a study area, the multi-source remote sensing data comprising multi-temporal optical satellite images, multi-temporal SAR satellite images and auxiliary data, the auxiliary data comprising a digital elevation model, meteorological observation data and earthquake catalog data; performing radiation calibration, atmospheric correction and terrain correction on the multi-temporal optical satellite images, performing precise orbit correction, thermal noise removal, radiation calibration and adaptive filtering on the multi-temporal SAR satellite images, and uniformly resampling all auxiliary data to a consistent space-time reference to obtain a standardized multi-source remote sensing data set that is consistent in space and time; processing the multi-temporal optical satellite images in the multi-source remote sensing data set by means of an improved frequency domain cross-correlation algorithm, including image fine registration, pixel offset calculation, adaptive window adjustment and system error correction, to obtain a horizontal displacement field of the study area; based on the multi-temporal SAR satellite images and the digital elevation model in the multi-source remote sensing data set, obtaining a time-series deformation field of the study area by means of optimized interferometric pair selection, atmospheric phase delay and orbit error correction, geometric distortion identification and processing, and time-series deformation robustness solving; based on the horizontal displacement field and the time-series deformation field, performing spatial adaptive grid fusion, Gaussian process-based time interpolation and dynamic allocation of optical and SAR data fusion weights according to deformation level, terrain conditions and data quality in a unified space-time framework to obtain a spatio-temporally continuous comprehensive deformation field; based on the comprehensive deformation field, performing feature extraction and dimension reduction by principal component analysis, and automatically classifying deformation patterns by means of unsupervised K-means clustering algorithm, and defining activity level by combining the category conversion of pre-seismic and post-seismic deformation patterns to obtain a deformation type map and an activity classification map of each observation point in the study area.

2. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, When performing terrain correction on the multi-temporal optical satellite images, a C correction model is used, and its expression is as follows: ； wherein, denotes a terrain-corrected radiance value, denotes an original radiance value, denotes a solar zenith angle, denotes an incidence angle, denotes a correction parameter determined by statistical regression on terrain flat areas; When performing adaptive filtering on the multi-temporal SAR satellite images, the filter window size is dynamically calculated and adjusted according to the local coherence coefficient, and the calculation formula is as follows: ； wherein, denotes a filter window size, and denotes a preset minimum window size and a preset maximum window size, respectively, denotes a local coherence coefficient.

3. The optical and SAR based high lying loose material multi-scale deformation monitoring method according to claim 1, characterized in that, The processing flow of the pixel offset calculation comprises: performing two-dimensional fast Fourier transform on the corresponding windows of the master image and the slave image to obtain the frequency spectrum of the master image and the slave image, calculating the normalized cross-power spectrum according to the frequency spectrum of the master image and the slave image, performing inverse Fourier transform on the normalized cross-power spectrum to obtain a correlation surface in the spatial domain, and determining the pixel-level displacement by finding the peak value of the correlation surface; in the neighborhood of the peak value position determined in the pixel-level displacement, a two-dimensional Gaussian surface model is used to perform least squares fitting on the discrete correlation surface intensity values to solve the sub-pixel peak value position as the final pixel offset; the calculation formula of the normalized cross-power spectrum is as follows: ； wherein, denotes the normalized cross-power spectrum, and denotes the frequency domain coordinate variable, denotes the frequency spectrum of the primary image, denotes the frequency spectrum of the secondary image, denotes the complex conjugate of denotes a small quantity that prevents division by zero; the two-dimensional Gaussian surface model is as follows: ； wherein, represents the relevant face intensity value, and represents the spatial domain coordinate variable, is the sub-pixel peak position, represents the peak amplitude, and respectively represent the standard deviation of the Gaussian distribution in and directions, represents the background noise level, represents the natural exponential function.

4. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, the processing flow of the adaptive window adjustment comprises: in the pixel offset calculation process, the local standard deviation of each potential window position is calculated in real time, and its deformation gradient module is estimated, and the window size for cross-correlation calculation at this position is dynamically calculated and determined according to the local standard deviation and the deformation gradient module, and the calculation formula is as follows: ； wherein, represents the window size after adaptive adjustment, represents the reference window size, represents the local standard deviation, represents the minimum standard deviation threshold, represents the deformation gradient modulus, represents the maximum expected deformation, and represents the weight coefficient; the processing flow of the system error correction comprises: After the initial displacement field is calculated, the initial displacement field is filtered using a non-local mean filtering algorithm, the non-local mean filtering algorithm being as follows: ； wherein, represents a pixel offset of the filtered th pixel point, , and respectively represent a neighborhood block centered at the pixel and , represents a filter strength parameter, represents a normalization factor, represents a pixel offset of the unfiltered th pixel point.

5. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, The preferred processing procedure for the interference pairs includes: Based on multi-temporal SAR satellite images, for all possible combinations of interference pairs, the weight of each interference pair is calculated based on its spatial baseline, time baseline and predicted coherence, and an optimized interference pair set is constructed according to the weight for time series deformation robustness solving; The weight calculation formula of the interference pair is as follows: ； wherein, denotes a weight of an interference pair, denotes a predicted coherence, denotes a vertical spatial baseline, denotes a critical baseline, denotes a time baseline of an interference pair, denotes a characteristic decorrelation time constant, denotes a natural exponential function; The processing procedure of the atmospheric phase delay includes: Based on the digital elevation model and the interference phase map, the regression coefficient is estimated through the phase-elevation statistical relationship of the stable area, and the terrain-related atmospheric delay signal is calculated and removed using the elevation correlation model; the elevation correlation model is as follows: ； wherein, denotes the atmospheric delay, denotes the terrain elevation, denotes the regression coefficients, denotes the residual turbulent atmosphere phase component; The processing procedure of the geometric distortion identification includes: Based on the digital elevation model and the SAR satellite orbit parameters, the exponential value of each pixel is calculated When the exponential value is less than a predetermined threshold, the pixel is determined as a shadow area, and when the exponential value is greater than the predetermined threshold, the pixel is determined as an overlap area. When the exponential value is less than a predetermined threshold, the pixel is determined as a shadow area, and when the exponential value is greater than the predetermined threshold, the pixel is determined as an overlap area. The identified shadow areas and overlap areas are shielded. The The formula for calculating the index value is as follows: ； wherein, denotes exponent value, denotes terrain slope, denotes aspect, denotes radar incidence angle, denotes satellite azimuth.

6. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 5, characterized in that, The processing procedure of the time series deformation robustness solving includes: Based on the optimized set of interference pairs, the observation equation is constructed: ; The deformation rate vector is solved by using an iterative weighted least square estimation method: ; wherein, denotes the unwrapped phase vector, denotes the design matrix, denotes the deformation rate vector, denotes the observation noise, denotes the weight matrix, denotes the smoothing constraint matrix, denotes the regularization parameter, denotes the matrix transpose.

7. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, The processing procedure of the spatial adaptive grid fusion includes: The evaluation function value corresponding to different grid sizes is calculated, the optimal grid size is determined by finding the optimal value of the evaluation function, the balance between data coverage and spatial resolution is achieved, and the corresponding evaluation function is as follows: ； wherein denotes the grid size, denotes the grid size the corresponding evaluation function value, denotes the number of grids containing observation points, denotes the total number of grids, denotes the characteristic scale, denotes the weight coefficient, denotes the natural exponential function; After the optimal grid size is determined, in each grid unit, the optical and InSAR deformation variables are fused by adopting distance inverse weight interpolation: ； ； wherein, denotes the deformation after fusion, denotes the weight of the th observation point, denotes the deformation of the th observation point, denotes the observation accuracy of the th observation point, denotes the distance of the th observation point to the center of the grid, denotes the smoothing parameter; In the Gaussian process-based time interpolation, the deformation time evolution is represented as: ； where, denotes the deformation at time denotes the mean function, denotes the Gaussian process, denotes the covariance function, and denotes the covariance function at different time points, is a Matérn kernel:​ ； wherein, denotes the signal variance, denotes the characteristic time scale; In dynamically assigning the fusion weights of optical and SAR data according to the deformation variable level, terrain conditions and data quality, the optical data weight and the SAR data weight are calculated and determined according to the following formula: ； ； wherein, denotes the optical data weight, denotes the SAR data weight, denotes the deformation rate, denotes the deformation rate threshold, denotes the terrain slope, denotes the deformation rate characteristic parameter, denotes the slope characteristic parameter, and denote the quality indicators of the optical and SAR data, respectively.

8. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, Feature extraction and dimension reduction are performed by principal component analysis, including: The time series deformation data corresponding to the deformation matrix of the comprehensive deformation field are standardized: ； wherein, denotes the deformation of the th observation point at the th time sample, denotes the normalized deformation, denotes the mean of all time sample deformations of the th observation point, denotes the standard deviation of all time sample deformations of the th observation point; A covariance matrix is constructed, the covariance matrix is decomposed to obtain eigenvalues and eigenvectors, and the covariance matrix is as follows: ； wherein, denotes the covariance matrix, denotes the deformation matrix corresponding to the normalized deformation variable, denotes the matrix transpose, denotes the number of observation points; The deformation matrix corresponding to the standardized deformation variable is projected onto the eigenvectors to obtain the principal components: ； wherein represents the first principal component, represents the first eigenvector obtained by eigen decomposition. 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: ； wherein, denotes the first eigenvalue obtained by eigen decomposition, denotes the number of principal components, denotes the total number of eigenvalues; When the unsupervised K-means clustering algorithm is used, the procedure for determining the optimal number of clusters includes: The range of the candidate number of clusters is set, K-means clustering is performed for each candidate number, the deformation variable of each observation point is taken as a sample point, the silhouette coefficient of all sample points is calculated and the average value is obtained, the candidate number that makes the average silhouette coefficient maximum is selected as the final number of clusters, and the calculation formula of the silhouette coefficient is as follows: ； wherein, denotes the silhouette coefficient of the th sample point, denotes the average distance of the th sample point to other sample points of the same cluster, denotes the average distance of the th sample point to all sample points of the nearest neighbor cluster; In the initialization of the unsupervised K-means clustering, the selection probability of each sample point is calculated, and the initial cluster center is selected according to the probability, and the calculation formula of the probability is as follows: ； wherein, denotes the probability that the th sample point is selected as the next cluster center, denotes the Euclidean distance of the sample point to the nearest cluster center, denotes the Euclidean distance of the th sample point to the nearest cluster center, denotes the number of observations; In the iteration and updating process of the unsupervised K-means clustering algorithm, the cluster center is updated according to the weighted average of the member points: ； wherein, denotes the center point coordinate of the th cluster, denotes the weight assigned to the th sample point, the weight being calculated for each sample point according to its measurement uncertainty and the reference accuracy: ； wherein, represents the measurement uncertainty of the th sample point, represents the reference precision.

9. The optical and SAR-based high-relief geomorphological mass movement monitoring method according to claim 1, characterized in that, The activity level includes the following four categories: Non-sensitive area: stable before earthquake-stable after earthquake; Self-stable area: active before earthquake-stable after earthquake; Continuously active area: active before earthquake-active after earthquake; Seismic activation area: stable before earthquake-active after earthquake.

10. The optical and SAR-based high-relief geomorphological deformation monitoring method according to any one of claims 1 to 9, characterized in that, The method further includes: Based on the activity classification map, meteorological observation data and the comprehensive deformation field, quantitative analysis of multi-factor triggering mechanism and time lag correlation analysis are performed, critical state recognition is performed based on Fukuzono-Voight model, and warning information is generated according to preset classification warning criteria; The quantitative analysis includes: precipitation triggering effect modeling, temperature triggering effect modeling and earthquake triggering effect evaluation; The precipitation triggering effect modeling includes: based on a sliding time window, a quantitative relationship between deformation acceleration and cumulative precipitation is established, and the model is as follows: ； wherein, denotes denotes the deformation acceleration at the time instant, denotes the background acceleration, denotes the accumulated precipitation within the sliding time window, denotes the precipitation response time lag, denotes the precipitation influence coefficient; The temperature triggering effect modeling includes: calculating a positive accumulated temperature index for quantifying the intensity of freeze-thaw cycles: ； wherein, denotes the positive accumulated temperature index, denotes the daily mean temperature, denotes the thawing threshold, denotes the summation period; The earthquake triggering effect evaluation includes: using a modified Newmark displacement model, combining peak ground acceleration, duration and initial stability conditions of high-risk areas, evaluating the cumulative impact of earthquake events on the stability of loose bodies, and generating a seismic influence factor; The process of the time lag correlation analysis includes: Based on the comprehensive deformation field and meteorological observation data, the cross-correlation coefficient of deformation time series and triggering factor time series under different precipitation response time lags is calculated, and the best precipitation response time lag is determined by finding the maximum cross-correlation coefficient; the calculation formula of the cross-correlation coefficient is as follows: ； wherein, denotes the cross-correlation coefficient, denotes the deformation time series, denotes the trigger factor time series, and denotes the corresponding mean; A multi-factor comprehensive triggering index is constructed, the influence of different triggering factors under the best precipitation response time lag is normalized and weighted fusion: ； wherein, denotes the integrated trigger index, and denote the normalized precipitation index and temperature index, respectively, denotes the earthquake impact factor, denote the corresponding weights, respectively. The process of the critical state recognition includes: based on the comprehensive deformation field, the reciprocal of the deformation rate time series data of the target area is calculated in the sliding time window, and linear fitting is performed using the Fukuzono-Voight model; when the slope of the fitted straight line is negative and the goodness of fit is greater than the goodness threshold, it indicates that the deformation is accelerating and the system is tending to be unstable; the Fukuzono-Voight model is as follows: ； wherein, denotes the rate of deformation, denotes time, and denotes the model parameters obtained by fitting the deformation data within the sliding window.

Images