Slope stability prediction method based on deformation monitoring

Abstract

The invention discloses a slope stability prediction method based on deformation monitoring, and relates to the field of slope monitoring, and the method comprises the steps: obtaining a time sequence deformation field of a target slope through employing an InSAR; generating a time sequence Dn through differential interference processing; collecting rainfall data; generating a time sequence representing regional rainfall conditions; grey correlation analysis is adopted, and the correlation degree of Dn and Pn is calculated; constructing a slope slip prediction model based on InSAR displacement and rainfall factors; if the displacement time sequence has an acceleration trend or fluctuation, marking the corresponding slope unit as an early warning analysis object; obtaining rainfall forecast data in a future period of time, and calculating a displacement predicted value di '(N + 1) of the marked side slope unit at the (N + 1) th time point in the future; and comparing the displacement predicted value di '(N + 1) of the (N + 1) th time point with a threshold value, and if di' (N + 1) exceeds the threshold value, giving a slippage early warning level. Aiming at low precision of a rainfall-induced slope slippage early warning model in the prior art, the precision is improved.

Description

Technical Field

[0001] The present application relates to the field of slope monitoring, and in particular to a slope stability prediction method based on deformation monitoring. Background Art

[0002] In recent years, extreme rainfall events caused by global climate change have occurred frequently, resulting in a significant increase in rainfall-induced slope landslide disasters, posing a serious threat to the safety of people's lives and property and the sustainable development of the social economy. Timely and accurate early warning of slope landslide disasters caused by rainfall is of great significance for minimizing disaster losses and ensuring regional safety and stability.

[0003] At present, a lot of research has been carried out on the early warning of rainfall-induced slope slip, mainly focusing on rainfall threshold analysis and physical and mechanical mechanism research. Among them, rainfall threshold analysis constructs two-dimensional critical threshold curves such as rainfall-duration or rainfall intensity-accumulated rainfall by statistically analyzing historical rainfall events and slope instability cases, which is used to determine whether the slope may slip under a given rainfall scenario. However, this type of method is generally based on statistical analysis at the regional scale, and it is difficult to consider the individual differences of the slope. The pertinence and reliability of the early warning need to be further improved.

[0004] The research on physical and mechanical mechanisms focuses on analyzing the disaster-causing mechanisms such as slope soil softening and pore water pressure changes caused by rainfall infiltration, and on this basis, a slope stability calculation model is established. This type of method requires a large number of engineering geological and hydrogeological parameters as input, and the temporal and spatial distribution information of the parameters is difficult to accurately describe, resulting in difficulty in ensuring the accuracy of model prediction. In addition, these two types of methods mainly rely on rainfall data as early warning criteria, lack effective monitoring of the actual deformation response of the slope, and are difficult to achieve quantitative slip warning and risk level classification. Summary of the invention

[0005] In response to the problem of low accuracy of rainfall-induced slope slip warning models in the prior art, the present application provides a slope stability prediction method based on deformation monitoring, which uses INSAR to obtain spatiotemporal data of slope deformation, combines regional rainfall observation data, and uses grey correlation analysis to screen out slope parts that are sensitive to rainfall response, and establishes a slope slip prediction model. At the same time, trend analysis and periodicity analysis are used to identify potential unstable slopes, thereby improving the accuracy of slope warning.

[0006] The purpose of this application is achieved through the following technical solutions.

[0007] The present application provides a slope stability prediction method based on deformation monitoring, comprising: S1, using synthetic aperture radar interferometry InSAR technology to monitor the target slope and obtain the time series deformation field of the target slope; generating a time series representing the radar line of sight displacement of the slope surface through differential interferometry processing;

[0008] D n ={d 1n ,d 2n ,......,d in ,......,d In}, where d in is the cumulative displacement of slope unit i at the nth time point; S2, multiple rain gauges are deployed on the target slope and adjacent areas to collect rainfall data; a time series P representing the regional rainfall conditions is generated n ={p 1n ,p 2n ,......,p jn ,......,p Jn}, where p jn is the rainfall recorded by rain gauge j in the nth sampling period; S3, using grey correlation analysis, calculate the displacement sequence D of the slope unit n With regional rainfall series P n According to the correlation degree, the slope units with a correlation degree with rainfall greater than a threshold are screened out, and a slope slip prediction model based on InSAR displacement and rainfall factor is constructed; S4, the displacement time series data of the slope units screened out in step S3 are extracted, and the trend analysis method is used to determine whether there is an acceleration trend in the time series, and the periodic analysis method is used to determine whether there is fluctuation in the displacement time series; if there is an acceleration trend or fluctuation in the displacement time series, the corresponding slope unit is marked as an early warning analysis object, and the process goes to step S5; otherwise, the process returns to step S1 to continue monitoring the target slope; S5, the rainfall forecast data for a period of time in the future is obtained, and the rainfall forecast data is input into the slope slip prediction model constructed in step S3 to obtain the displacement response prediction result of the slope unit marked in step S4 under the predicted rainfall scenario; the displacement prediction value d of the marked slope unit at the N+1th time point in the future is calculated. i '(N+1); S6, the displacement prediction value d of the slope unit at the N+1th time point in the future obtained in step S5 i '(N+1) is compared with the threshold value. If the displacement prediction value d i '(N+1) exceeds the threshold, then according to the displacement prediction value d i '(N+1) gives the slip warning level of the corresponding slope unit, otherwise returns to step S1.

[0009] Further, S1, the target slope is monitored by using the synthetic aperture radar interferometry technology InSAR to obtain the time series deformation field of the target slope, including: obtaining at least two synthetic aperture radar complex images of the target slope at different times; using the synthetic aperture radar complex image with the earliest acquisition time as the main image, and using the remaining synthetic aperture radar complex images acquired later than the main image as auxiliary images; using the main image as a reference, aligning the auxiliary images in the spatial coordinate system; performing cross-correlation calculation on the aligned main image and each auxiliary image to obtain a correlation coefficient representing the coherence degree of the main image and each auxiliary image. coherence coefficient γ; select the primary and auxiliary image pairs according to the coherence coefficient γ to generate an interference complex map; perform phase unwrapping on the interference complex map to obtain the radar line-of-sight displacement Δr of the target slope in the corresponding time period of the primary and auxiliary images; convert the radar line-of-sight displacement Δr into a vertical displacement Δh; repeat the phase unwrapping steps to obtain the vertical displacement map of the target slope in different time periods; according to the time interval of the primary and auxiliary image pairs, convert the vertical displacement map in each time period into a displacement time series with time as the independent variable; align the displacement time series with the geographic coordinates of the target slope to generate a time series deformation field.

[0010] Furthermore, the complex interference map is generated, including: calculating the complex value S of the registered main image 1 and the complex value S of each auxiliary image 2 The mutual correlation coefficient between them is taken as the coherence coefficient γ, and the calculation formula of the coherence coefficient γ is: Where E[*] represents mathematical expectation, |*| represents modular operation, and * represents complex conjugate operation; set the coherence coefficient threshold γ t ; When the coherence coefficient γ is greater than the threshold γ t When the corresponding main and auxiliary images are taken as the main and auxiliary image pairs, the corresponding interference complex map is calculated for each main and auxiliary image pair. The calculation formula of the interference complex map is: I = S 1 ×S 2 *, where I represents the complex interference pattern, S 1 Represents the complex value of the main image, S 2 * represents the complex conjugate value of the auxiliary image; the calculated complex interference images are integrated to obtain the complex interference image set of the main image and the auxiliary image.

[0011] Furthermore, the interferometric complex image is phase unwrapped to obtain the radar line-of-sight displacement of the target slope in the corresponding period of the primary and auxiliary images; the radar line-of-sight displacement is converted into a vertical displacement, including: using the minimum cost flow algorithm to phase unwrap each interferometric complex image in the interferometric complex image set; estimating the absolute phase value Φ of the interferometric complex image; and calculating the radar line-of-sight displacement Δr of the target slope in the corresponding period of the primary and auxiliary images according to the absolute phase value Φ of the interferometric complex image: Where λ is the wavelength of the radar signal. Using the digital elevation model DEM data, the radar line of sight displacement Δr is converted into the vertical displacement Δh: Among them, θ is the radar incident angle, and α is the surface slope angle; the vertical displacement Δh is geocoded to generate the vertical displacement map in the corresponding time period of the primary and auxiliary images.

[0012] Furthermore, the absolute phase value of the interference complex image is estimated, including: using the Goldstein filtering algorithm to perform spectrum filtering on the interference complex image to obtain the interference complex image I'; according to the filtered interference complex image I', the phase gradient between adjacent pixels is calculated, and a cost function is constructed: f(m,n)=min(|Φ(m,n)-2kπ|), wherein Φ(m,n) is the phase difference between the pixel (m,n) and its adjacent pixels, and k is an integer; the phase unwrapping problem is converted into a minimum cost flow problem, with the pixels of the interference complex image I' as nodes and the cost function value f(m,n) as the capacity of the edge, and the minimum cost maximum flow is solved to obtain the global optimal phase unwrapping result φ(m,n); according to the phase unwrapping result φ(m,n), the phase difference Φ(m,n) is calculated, and the calculation formula is: Φ(m,n)=φ(m,n)+2kπ, wherein k is an integer.

[0013] Furthermore, a slope slip prediction model based on InSAR displacement and rainfall factor is constructed, including: the slope unit displacement time series D obtained in step S1 n ={d 1n ,d 2n ,......,d in ,......,d In} is dimensionless, and the normalized sequence of slope unit displacement time series is obtained as X = {x 1 ,x 2 ,......,x i ,......,x I},in: i=1,2,......,I;for the regional rainfall time series P obtained in step S2 n ={p 1n ,p 2n ,......,p jn ,......,p Jn} is dimensionless, and the normalized sequence of regional rainfall time series is obtained as Y = {y 1 ,y 2 ,......,y i ,......,y J},in: j=1,2,......,J; calculate the grey correlation coefficient ξ(i,j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall. The calculation formula is: Among them, ρ is the resolution coefficient, and its value range is (0, 1); the grey correlation degree r(i, j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall is calculated, and the calculation formula is: Where i represents the index time node of the slope unit displacement time series; j represents the index time node of the regional rainfall time series; k represents the comparison index used in the grey correlation calculation; N is the sequence length; the correlation threshold r is set t , the grey relational degree r(i,j) is greater than the threshold r t The slope units are screened as those with greater correlation with regional rainfall, forming the rainfall sensitive unit set R s Extract the rainfall sensitive unit set R from the slope unit displacement time series data obtained in step S1 and the regional rainfall time series data obtained in step S2 s The displacement time series data D' and the corresponding regional rainfall time series data P' of each slope unit are extracted, and the extracted displacement time series data D' and rainfall time series data P' are input into the support vector machine SVM model to build a slope slip prediction model based on InSAR displacement and rainfall factors.

[0014] Furthermore, a slope slip prediction model based on InSAR displacement and rainfall factor is constructed. The model is: Where d'(t) is the predicted displacement value of the slope unit at the tth time point in the future; L is the number of support vectors; w l is the weight coefficient corresponding to the lth support vector;

[0015] K[x(t),x l ] is the Gaussian radial basis kernel function RBF: x(t) is the input feature vector at the t-th time point in the future, including the displacement of the slope unit at the previous m time points and the regional rainfall at the previous n time points, that is: x(t) = [d(t-1), d(t-2), ..., d(tm); p(t), p(t-1), ..., p(tn)], x l is the input feature vector corresponding to the lth support vector; ||x(t)-x l || represents the Euclidean distance; σ is the parameter of the RBF kernel function, which controls the width of the kernel function; b is the bias term; ε is the model error term.

[0016] Further, S4, extracting the displacement time series data of the slope units screened out in step S3, using a trend analysis method to determine whether the time series has an acceleration trend, and using a periodic analysis method to determine whether the displacement time series has fluctuations, including: s In the above example, the displacement time series data of each slope unit are extracted to form the displacement time series data set D to be analyzed. s ; Local polynomial regression method is used to analyze the data set D s Each displacement time series data in is fitted to obtain the smooth curve f(t) corresponding to each displacement time series, and all smooth curves f(t) constitute a trend feature set F; the slope value of each smooth curve f(t) in the trend feature set F at each time point is calculated to obtain a slope sequence f'(t), and the slope sequence f'(t) is differentially operated to obtain an acceleration sequence f"(t), and all slope sequences f'(t) and acceleration sequences f"(t) constitute an acceleration feature set A; according to a preset slope threshold s 1 and acceleration threshold a 1 , the trend of the slope sequence f'(t) and acceleration sequence f"(t) of each slope unit is judged: if in a certain time period [t 1 ,t 2 ], the slope sequence f'(t) is greater than the slope threshold s 1 , and the acceleration sequence f”(t) is greater than the acceleration threshold a 1 , the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed; if the slope sequence f'(t) and the acceleration sequence f"(t) do not exceed their respective thresholds at the same time within the entire monitoring time range, the displacement time series data of the slope unit is passed to the subsequent spectrum analysis; the Fourier transform method is used to perform spectrum analysis on the displacement time series data of the slope unit that does not show an acceleration trend, and the frequency characteristics of each displacement time series are obtained, and all frequency characteristics are used to form a frequency feature set P s ; According to the preset frequency threshold f 1 , for the frequency feature set P s The frequency characteristics of each slope unit in the displacement time series are judged: if the frequency corresponding to the peak of the displacement time series spectrum is greater than the frequency threshold f 1 , then the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed; if the frequency corresponding to the peak value of the displacement time series spectrum is less than or equal to the frequency threshold f 1 , the corresponding slope unit is marked as a stable object, and the process returns to step S1 to continue monitoring the target slope.

[0017] Furthermore, the local polynomial regression method is used to classify the dataset D sFitting each displacement time series data in includes: setting the initial smoothing parameter a and the polynomial order q, defining the Gaussian weight function W; fitting the displacement time series data set D s At each time point t i , calculate the time point t j Relative to t i The time interval Δt ij , according to the Gaussian weight function W and the time interval Δt ij Calculate time point t j Relative to t i The weight value w ij , get the time point t i The weight matrix W i ; Using the weighted least squares method, at time point t i Based on the weight matrix W i For the displacement time series dataset D s Perform q-order local polynomial regression to obtain time point t i The local polynomial regression model M at i ; Based on the local polynomial regression model M i , calculate the time point t i The fitted value f at i and the fitted residual e i ; Traverse the displacement time series dataset D s For each time point in , the fitting values ​​and fitting residuals of all time points are obtained; according to the fitting residual sequence, the smoothing parameter a and the polynomial order q are optimized by the cross-validation method to obtain the optimal smoothing parameter a and the optimal polynomial order q; using the optimal smoothing parameter a and the optimal polynomial order q, the above polynomial fitting steps are repeated to obtain the final fitting values ​​of all time points, forming the smooth curve f(t) of the displacement time series.

[0018] Further, S5, obtain the displacement response prediction result of the slope unit marked in step S4 under the predicted rainfall scenario; calculate the displacement prediction value d of the marked slope unit at the N+1th time point in the future i '(N+1), including: obtaining the rainfall forecast data from the N+1th time point to the N+Mth time point in the future in the area where the target slope is located, forming a rainfall forecast sequence P'={p'(N+1),p'(N+2),......,p'(N+M)}; where M is the forecast time span; taking the rainfall forecast sequence P' as input, using the constructed slope slip prediction model based on InSAR displacement and rainfall factor, obtaining the displacement response prediction results of the marked warning analysis object at the future M time points, forming a displacement prediction sequence D i '={d i '(N+1),d i'(N+2),......,d i '(N+M)}; where d i '(n) is the predicted displacement value of slope unit i at the nth time point in the future, and the calculation formula is:

[0019] Among them, N s represents the number of support vectors, representing the number of samples used for model training; x'(n) is the input feature vector at the nth time point in the future, including the actual displacement d of slope unit i at the Nth time point i (N), the displacement prediction value d at the future n-1th time point i '(n-1) and the predicted rainfall at the nth time point in the future p'(n), that is: x'(n) = [d i (N),d i '(n-1),p'(n)],n=N+2,N+3,......,N+M, the input feature vector at the N+1th time point in the future: x'(N+1)=[d i (N),p'(N+1)], predicting the sequence D from the displacement i 'Extract the predicted displacement value d of slope unit i at the N+1th time point in the future i '(N+1), as the key monitoring object for early warning analysis.

[0020] Beneficial effects: Compared with the prior art, the advantages of this application are:

[0021] The introduction of InSAR monitoring technology to obtain high-precision, high-resolution time-series deformation fields on the slope surface overcomes the shortcomings of traditional slope monitoring methods, such as limited coverage and sparse measurement points, and provides a data basis for finely characterizing the slope deformation response characteristics. The slope radar line-of-sight displacement is obtained through differential interferometry processing, and the conversion of radar line-of-sight displacement and vertical displacement is realized in combination with the digital elevation model, eliminating the displacement estimation error caused by terrain undulations and improving the accuracy of slope deformation monitoring.

[0022] The observation data of multiple rainfall stations inside the slope are comprehensively considered to construct a time series reflecting the regional rainfall conditions. The grey correlation analysis method is used to quantitatively evaluate the correlation between the slope displacement response and the rainfall factor, screen out the slope units sensitive to rainfall, and then construct a slope slip prediction model in a targeted manner. This strategy fully explores the intrinsic relationship between slope deformation and rainfall factors, overcomes the problem of poor universality of existing rainfall warning models, and improves the accuracy of slope slip warning.

[0023] The local polynomial regression method is introduced to smooth the slope displacement time series, which effectively suppresses the interference of monitoring noise and outliers. Based on the smooth curve, trend analysis is used to determine whether the displacement time series has acceleration characteristics, and Fourier transform is used for periodic analysis to identify slope units with the risk of fluctuating instability. This method describes the dynamic evolution characteristics of slope deformation from the perspectives of deformation trend and periodicity, providing a reliable basis for identifying slope instability warning objects.

[0024] The rainfall forecast data is introduced into the slope slip prediction model to predict the future displacement response of the selected early warning analysis objects and give quantitative displacement prediction values. This strategy makes full use of the rainfall forecast information provided by the meteorological department, expands the early warning time scale, and wins precious time for timely response measures to dangerous situations. At the same time, the prediction model trained with the historical monitoring data of the slope can simulate the dynamic response of the slope displacement according to the future predicted rainfall scenario, provide quantitative early warning indicators for slip, and provide a reference for slope stability evaluation. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] The present application will be further described in the form of exemplary embodiments, which will be described in detail by the accompanying drawings. These embodiments are not restrictive, and in these embodiments, the same number represents the same structure, wherein:

[0026] Figure 1 is an exemplary flow chart of a slope stability prediction method based on deformation monitoring according to some embodiments of the present application;

[0027] Figure 2 is an exemplary flow chart of obtaining a temporal deformation field according to some embodiments of the present application;

[0028] Figure 3 is an exemplary flow chart of constructing a slope slip prediction model according to some embodiments of the present application;

[0029] Figure 4 It is an exemplary flow chart of obtaining early warning analysis slope according to some embodiments of the present application. DETAILED DESCRIPTION

[0030] The method and system provided in the embodiments of the present application are described in detail below with reference to the accompanying drawings.

[0031] Figure 1The present invention is an exemplary flow chart of a slope stability prediction method based on deformation monitoring according to some embodiments of the present invention. The slope stability prediction method based on deformation monitoring includes: using synthetic aperture radar interferometry InSAR to monitor the target slope and obtain the time series deformation field of the target slope; generating a time series D representing the radar line of sight displacement of the slope surface through differential interferometry processing; n ={d 1n ,d 2n ,......,d in ,......,d In}, where d in is the cumulative displacement of slope unit i at the nth time point.

[0032] Multiple rain gauges are deployed on the target slope and adjacent areas to collect rainfall data; a time series P representing the regional rainfall conditions is generated. n ={p 1n ,p 2n ,......,p jn ,......,p Jn}, where p jn is the rainfall recorded by rain gauge j in the nth sampling period; the displacement sequence D of the slope unit is calculated by using grey correlation analysis. n With regional rainfall series P n The correlation degree is, where I represents the total number of slope units, D n It represents the cumulative displacement of all slope units at the nth time point, P n represents the rainfall sequence recorded by all rain gauges at the nth time point, and J represents the total number of rain gauges; according to the correlation degree, the slope units with a correlation degree with rainfall greater than the threshold are screened out, and a slope slip prediction model based on InSAR displacement and rainfall factor is constructed; the displacement time series data of the screened slope units are extracted, and the trend analysis method is used to determine whether there is an acceleration trend in the time series, and the periodic analysis method is used to determine whether there is fluctuation in the displacement time series; if there is an acceleration trend or fluctuation in the displacement time series, the corresponding slope unit is marked as an early warning analysis object; otherwise, the target slope is returned to continue monitoring; the rainfall forecast data for a period of time in the future is obtained, and the rainfall forecast data is input into the constructed slope slip prediction model to obtain the displacement response prediction results of the marked slope units under the predicted rainfall scenario; the displacement prediction value d of the marked slope unit at the N+1th time point in the future is calculated i '(N+1); S6, the displacement prediction value d of the slope unit at the N+1th time point in the future will be obtained i '(N+1) is compared with the threshold value. If the displacement prediction value d i '(N+1) exceeds the threshold, then according to the displacement prediction value di '(N+1) gives the slip warning level of the corresponding slope unit, otherwise returns to continue monitoring.

[0033] Data acquisition stage: Use synthetic aperture radar interferometry InSAR to monitor the target slope and obtain the time-series deformation field of the target slope; select appropriate radar satellite platforms, such as TerraSAR-X, COSMO-SkyMed, etc., to ensure that the satellite's orbital parameters and imaging mode meet the requirements of InSAR monitoring. According to the location and size of the target slope, determine the time span and spatial coverage of the radar data. Generally, it is necessary to obtain at least 20 radar images, and the time span should cover the entire deformation process of the target slope. Preprocess the acquired radar data, including image registration, speckle noise filtering, atmospheric delay correction, etc., to improve data quality. Use InSAR technology to interfere with radar data and generate interference phase maps of the target slope at different times. Unwrap the interference phase map, convert the phase into the actual ground deformation, and obtain the time-series deformation field of the target slope.

[0034] Specifically, InSAR (Interferometric Synthetic Aperture Radar) is a technology that uses the phase information of radar signals to perform high-precision deformation monitoring of ground targets. It extracts the deformation information of ground targets in the direction of the radar line of sight by interferometric processing of radar images acquired at different times. The time series deformation field refers to the deformation of the target area over a period of time obtained using InSAR technology, which characterizes the deformation of different locations in the target area at different times. The time series deformation field usually uses time as the horizontal coordinate and the deformation as the vertical coordinate, reflecting the deformation evolution process of the target area.

[0035] Acquire at least two synthetic aperture radar complex images of the target slope at different times; select radar data from at least two different times according to the determined time span of the radar data. Perform focus processing on the selected radar data to obtain a Single Look Complex (SLC) image, i.e., a synthetic aperture radar complex image. Check the quality of the SLC image to ensure that the image is clear, has a high signal-to-noise ratio, and has no obvious missing or distortion. Sort the multiple SLC images obtained in order of acquisition time, determine the earliest acquired image as the main image, and determine the remaining images as auxiliary images. Using the main image as a reference, use an image registration algorithm such as the correlation coefficient maximization algorithm to register each auxiliary image in the spatial coordinate system. The purpose of registration is to make the position coordinates of the same target on the main and auxiliary images consistent. The registration process needs to consider the influence of factors such as orbital parameters and terrain undulations. Precise orbital data and external DEM data can be used for registration correction. To evaluate the accuracy of the registration results, it is generally required that the root mean square error of the main and auxiliary images after registration is less than 0.1 pixels. If the registration accuracy does not meet the requirements, it is necessary to adjust the registration parameters or algorithm and re-register until the accuracy requirements are met;

[0036] For each primary and secondary image pair after registration, calculate the primary image complex value S 1 and auxiliary image complex value S 2 The mutual correlation coefficient between them is taken as the coherence coefficient γ. The calculation formula of the coherence coefficient γ is: E[] represents mathematical expectation, that is, averaging the complex values. || represents modular operation, that is, calculating the modular value of a complex number. * represents complex conjugate operation, that is, negating the imaginary part of a complex number. The value range of the coherence coefficient γ is [0, 1]. The larger the γ value, the better the coherence of the main and auxiliary images. Set the coherence coefficient threshold γ t , generally takes a value of 0.3 to 0.8, which is determined according to data quality and research needs. For each primary and secondary image pair, calculate its coherence coefficient γ: if γ>γ t , then the primary and secondary image pairs are used as candidate image pairs for subsequent interference processing. t , then the main and auxiliary image pair is considered to have poor coherence and is not suitable for interference processing, so it is eliminated. For each main and auxiliary image pair selected, its interference complex map I is calculated. The calculation formula of the interference complex map I is: I = S 1 ×S 2 *, S 1 Represents the complex value of the main image, S 2 * indicates the complex conjugate value of the auxiliary image. The interferometric complex image I contains the phase difference information between the primary and auxiliary images, reflecting the ground deformation between the two observations. The interferometric complex images of all primary and auxiliary image pairs are integrated to obtain an interferometric complex atlas. The interferometric complex atlas contains the deformation information of the study area at different times and is the basic data for time series analysis.

[0037] The interferogram is phase unwrapped to obtain the radar line-of-sight displacement Δr of the target slope in the corresponding period of the main and auxiliary images; the radar line-of-sight displacement Δr is converted into the vertical displacement Δh. The spectrum filtering of the interferogram uses the Goldstein filtering algorithm, and the interferogram I is Fourier transformed: the interferogram I is divided into multiple overlapping sub-windows, each of which is N×N in size (commonly used sizes are 32×32 or 64×64). A two-dimensional Fourier transform (FFT) is performed on each sub-window to obtain the spectrum matrix F of the sub-window. Estimate the local spectrum amplitude and noise level: calculate the amplitude spectrum A and phase spectrum P of the spectrum matrix F of each sub-window. According to the amplitude spectrum A, estimate the signal-to-noise ratio (SNR) in the sub-window. Common estimation methods include: dividing the amplitude spectrum A into a signal area and a noise area, calculating the mean of the two areas respectively, and the ratio is the estimated signal-to-noise ratio. Assuming that the amplitude spectrum A follows an exponential distribution, the signal-to-noise ratio is estimated by the mean and variance of the amplitude spectrum. According to the estimated signal-to-noise ratio, the filter strength α of the sub-window is adaptively determined. The higher the signal-to-noise ratio, the smaller the filter strength. The value range of the filter strength α is [0, 1]. Construct a spectrum filter: According to the filter strength α, construct a spectrum filter H. The expression of the spectrum filter H is: H = A a , the function of the spectrum filter H is to enhance the frequency components with high signal-to-noise ratio and suppress the frequency components with low signal-to-noise ratio. Spectral filtering: multiply the spectrum filter H by the spectrum matrix F of the sub-window element by element to obtain the filtered spectrum matrix F'. F'=F×H, perform inverse Fourier transform (IFFT) on the filtered spectrum matrix F' to obtain the filtered sub-window interference complex image I sub '. Reconstruct the filtered interference complex image: The filtering results of all sub-windows I sub 'Perform overlapping averaging to obtain a complete post-filtering interference complex image I'. The purpose of overlapping averaging is to reduce the discontinuity of the sub-window edge and improve the smoothness of the filtering result.

[0038] Construct a cost function: Calculate the phase gradient between adjacent pixels: Perform phase unwrapping on the filtered interference complex image I' to obtain a continuous phase value. Calculate the phase difference ΔΦ(m,n) between each pixel (m,n) and its adjacent pixels. Commonly used adjacent pixels include four directions: up, down, left, and right. The calculation formula for the phase difference ΔΦ(m,n) is: ΔΦ(m,n) = Φ(m,n) - Φ(m',n'), where (m',n') is the adjacent pixel coordinate of pixel (m,n). Construct a cost function: For each pixel (m,n) and its adjacent pixel (m',n'), construct a cost function f(m,n).

[0039] The expression of the cost function f(m,n) is: Where k is an integer in the range of [-K, K], and K is the maximum integer number. The cost function f(m,n) represents the minimum cost required to untangle the phase difference ΔΦ(m,n) to 2kπ. When the phase difference ΔΦ(m,n) is close to 2kπ, the value of the cost function f(m,n) is small, indicating that the possibility of untangle to 2kπ is high.

[0040] Minimum cost flow solution: Convert the phase unwrapping problem into a minimum cost flow problem: Treat each pixel (m, n) of the interference complex graph I' as a node in the network. Establish directed edges between adjacent pixels, and the direction of the edge is the direction of the phase gradient. Use the cost function f(m, n) as the capacity (cost) of the corresponding directed edge. Construct a network flow graph: In the network flow graph, add a source point (source) and a sink point (sink). Establish directed edges between the source point and all pixel nodes, with infinite capacity and 0 cost. Establish directed edges between all pixel nodes and the sink point, with infinite capacity and 0 cost.

[0041] Solve the minimum cost maximum flow: Construct a network flow graph: According to the phase untangling problem, construct a network flow graph G = (V, E), where V is a node set and E is an edge set. The network flow graph includes source point s, sink point t and intermediate nodes, and the intermediate nodes represent the phase value of each pixel. An edge with a capacity of 1 is established between the source point s and each pixel node, indicating that each pixel can only be assigned an integer number of times. An edge with a capacity of ∞ is established between each pixel node and its adjacent pixel nodes, indicating the phase difference constraint between adjacent pixels. An edge with a capacity of 1 is established between each pixel node and the sink point t, indicating that the integer number of each pixel can only be 0 or 1. Assign a corresponding cost to each edge, indicating the cost required to select the edge. The cost can be set according to the phase gradient between adjacent pixels. The larger the gradient, the higher the cost. Initialization: Assign a distance label d(v) to each node, indicating the shortest distance from the source point s to the node v. Initially, d(s) = 0, and the distance labels of other nodes are ∞. Assign an excess flow e(v) to each node, indicating the flow surplus of node v. Initially, e(s) is the total flow of the network, and the excess flow of other nodes is 0. Create a first-in-first-out (FIFO) queue Q to store nodes with excess flow greater than 0.

[0042] Execute the Goldberg-Tarjan algorithm: Add source node s to queue Q. When queue Q is not empty, repeat the following steps: Take the head node u from queue Q. Perform a push operation on each outgoing edge (u, v) of node u: If node u has excess flow (e(u)>0), and the residual capacity of edge (u, v) is greater than 0, and d(u)=d(v)+1, then push the flow from node u to node v, update the flow and residual capacity of edge (u, v), and update the excess flow of nodes u and v.

[0043] If the excess flow of node v becomes positive and v is not in queue Q, then v is added to queue Q. If all outgoing edges of node u cannot push traffic, a relabeling operation is performed: the distance label of node u is updated to the minimum distance of all its outgoing edges to the node plus 1, that is,

[0044] d(u)=min{d(v)+1|(u,v)∈E,r(u,v)>0}. If the distance label of node u changes and u is not in queue Q, then u is added to queue Q. When queue Q is empty, the algorithm terminates and the minimum cost maximum flow is obtained.

[0045] Obtain the phase unwrapping result: According to the result of the minimum cost maximum flow, solve the corresponding integer order k for each pixel node (m, n). The integer order k can be determined by the flow from the pixel node (m, n) to the sink t in the network flow graph. If the flow is 1, then k = 1; if the flow is 0, then k = 0. Substitute the integer order k into the phase unwrapping formula to obtain the absolute phase value φ(m, n) after unwrapping: φ(m, n) = Φ(m, n) + 2kπ, where Φ(m, n) is the wrapped phase value of the pixel (m, n).

[0046] Calculate the absolute phase value: Determine the reference point: Select a pixel with a known phase value in the interference complex image as the reference point. Usually, a pixel with small phase noise and good coherence is selected. The phase value of the reference point is recorded as Φ ref Calculate the absolute phase value: For each pixel (m, n) in the interference complex image, according to the phase unwrapping result φ(m, n) and the reference point phase value Φ ref , calculate its absolute phase value Φ(m,n). The calculation formula is:

[0047] Φ(m,n)=φ(m,n)+Φ ref +2kπ, where k is an integer, is used to eliminate phase wraparound. The absolute phase value Φ(m,n) represents the actual phase value of the pixel (m,n) relative to the reference point, eliminating the effect of phase wraparound.

[0048] Calculate the radar line of sight displacement: Obtain radar parameters: Obtain the radar signal wavelength λ, which is determined by the transmission frequency of the radar system. The radar signal wavelength λ can be calculated by the following formula: λ = c / f, where c is the speed of light (approximately 3×10^8m / s) and f is the radar transmission frequency. Calculate the radar line of sight displacement: For each pixel (m, n) in the interference complex map, calculate the radar line of sight displacement Δr based on its absolute phase value Φ(m, n) and the radar signal wavelength λ. The calculation formula is: The radar line of sight displacement Δr represents the displacement of the pixel (m, n) in the radar line of sight direction, and the unit is distance (such as meter). The calculated radar line of sight displacement Δr is assigned to the corresponding pixel (m, n) to generate a radar line of sight displacement map.

[0049] Convert to vertical displacement: Obtain terrain parameters: Use the digital elevation model (DEM) to obtain the radar incident angle θ and the surface slope angle α corresponding to each pixel (m, n). The radar incident angle θ represents the angle at which the radar signal is incident on the surface, which can be calculated from the radar imaging geometry. The specific calculation method depends on the imaging mode and orbital parameters of the radar system. The surface slope angle α represents the inclination angle of the surface at the pixel (m, n), which can be calculated from the DEM data. The calculation steps are as follows: Perform the least squares plane fitting in the local window on the DEM data to obtain the plane equation coefficient at each pixel. According to the plane equation coefficient, calculate the surface slope angle α at each pixel. Specifically, obtain the digital elevation model (DEM) data of the study area, and the DEM data provides the elevation information of each pixel. Determine the size of the local window, usually select an odd-sized window such as 3×3, 5×5 or 7×7. The selection of the window size needs to be determined according to the resolution and terrain characteristics of the DEM data.

[0050] For each pixel point (m, n) in the DEM data, take the pixel as the center and select the pixel points in the local window for least square plane fitting. Assume that there are N pixels in the local window and the coordinates of each pixel point are (x i ,y i ,z i ), where i = 1, 2, ..., N. Construct the least squares plane fitting equation: z = ax + by + c, where a, b, c are the coefficients of the plane equation to be determined. Construct the error function: Represents the difference between the fitted plane and the actual elevation value. Solve the least squares problem, that is, find the values ​​of a, b, and c that minimize the error function E. This can be obtained by solving the following linear equations: Solve the linear equations to obtain the plane equation coefficients a, b, and c. Calculate the surface slope angle α: Calculate the surface slope angle α at each pixel point (m, n) based on the fitted plane equation coefficients a, b, and c. The surface slope angle α can be calculated using the following formula: Wherein, arctan is the inverse tangent function. The calculated surface slope angle α is assigned to the corresponding pixel point (m, n) to obtain a surface slope angle map.

[0051] Calculate the vertical displacement: For each pixel (m, n) in the interferometric complex image, calculate the vertical displacement Δh based on its radar line of sight displacement Δr, radar incident angle θ and ground slope angle α. The calculation formula is: The vertical displacement Δh represents the displacement of the pixel (m, n) in the vertical direction, and the unit is distance (such as meter). The calculated vertical displacement Δh is assigned to the corresponding pixel (m, n).

[0052] n), generate a vertical displacement map. Generate a vertical displacement map: Calculate the vertical displacement of all pixels in the interference complex map to obtain the vertical displacement Δh of each pixel. Assign the vertical displacement Δh to the corresponding pixel to generate a vertical displacement map. The vertical displacement map shows the vertical displacement distribution of the surface and can be used for further surface deformation analysis and application.

[0053] Generate a vertical displacement map, and use satellite orbit parameters and imaging geometry information to establish the correspondence between the interferometric complex image pixel coordinates (m, n) and the geographic coordinates (longitude, latitude). Common methods include: directly calculating the geographic coordinates corresponding to each pixel based on satellite orbit parameters and imaging models. Use ground control points (GCPs) to establish a conversion relationship between pixel coordinates and geographic coordinates, and calculate the geographic coordinates of other pixels by interpolation or fitting. According to the established correspondence between pixel coordinates and geographic coordinates, convert the pixel coordinates of the vertical displacement Δh into geographic coordinates. For each pixel (m, n), assign its vertical displacement Δh to the corresponding geographic coordinate position. Select a suitable map projection method, such as UTM projection, Gauss-Krüger projection, etc., and project the vertical displacement Δh in the geographic coordinate system to the target map coordinate system. The projection process needs to consider factors such as the earth's curvature and the geodetic datum to ensure the accuracy of the geometric position of the displacement. Generate a blank image of the same size as the interferometric complex image in the target map coordinate system. The vertical displacement Δh after projection is filled into the corresponding pixel position to form a vertical displacement map. Each pixel value in the vertical displacement map represents the vertical displacement of the position in the corresponding time period of the primary and auxiliary images, and the unit is distance (such as meters). Repeat the interference processing and phase unwrapping of each primary and auxiliary image pair to obtain the vertical displacement map of the target slope in different time periods. According to the imaging time of each primary and auxiliary image pair, the time mark of each vertical displacement map is determined. The vertical displacement maps in each time period are arranged in chronological order, and the displacement of each pixel position is extracted to form a displacement time series with time as the independent variable. The displacement time series is filtered and detrended to eliminate the influence of noise such as atmospheric delay and orbit error.

[0054] Figure 2According to the exemplary flow chart of obtaining the time series deformation field shown in some embodiments of the present application, the displacement time series is registered with the geographic coordinates of the target slope, and the time series deformation field is generated. The geographic coordinate information of the target slope area, including the latitude and longitude coordinates and elevation information of the control points, is obtained by using high-precision GPS measurement or other surveying and mapping methods. A digital terrain model (DEM) of the target slope is established to provide detailed terrain information and spatial reference. For each displacement time series image, its pixel coordinate information and corresponding displacement value are extracted. The conversion relationship between the pixel coordinates and the geographic coordinates of the displacement time series image is established using satellite orbit parameters, imaging geometry model and ground control point information. Common methods include: direct geocoding: using satellite orbit parameters and imaging geometry model to directly calculate the geographic coordinates corresponding to each pixel. Indirect geocoding: using ground control points to establish a mapping relationship between pixel coordinates and geographic coordinates, and calculating the geographic coordinates of other pixels by interpolation or fitting. Using the established correspondence between pixel coordinates and geographic coordinates, the position coordinates of each pixel in the displacement time series are registered with the geographic coordinates of the target slope. The following factors need to be considered in the registration process: Geometric distortion of the image: Geometric distortion of the image caused by factors such as satellite imaging geometry and terrain undulations needs to be corrected and rectified. Conversion of the coordinate system: Convert the coordinate system of the displacement time series with the geographic coordinate system of the target slope to ensure that both are in the same coordinate system. Registration accuracy: Improve the accuracy and reliability of the registration by selecting appropriate registration algorithms and parameters, such as least squares registration, feature point registration, etc. According to the registration results, map the displacement values ​​in the displacement time series to the geographic coordinate system of the target slope. For each geographic location on the target slope, extract the corresponding displacement time series to form a time series deformation field. The time series deformation field reflects the deformation state and deformation evolution process of the target slope at different time nodes.

[0055] Generate a time series of radar line-of-sight displacement of the slope surface through differential interferometry processing, select appropriate SAR data sources, such as TerraSAR-X, Sentinel-1, etc., to obtain multi-phase SAR data covering the target slope area. Ensure that the acquired SAR data is long enough in time span to capture the long-term evolution of slope deformation. The time interval of SAR data should be selected according to the rate of slope deformation and research needs, usually monthly, quarterly or annual time scales. Perform differential interferometry processing: Select a stable master image as the reference image for differential interferometry processing. Align each slave image with the master image to ensure that they are aligned in space. Generate a differential interferogram between each pair of master and slave images to obtain phase difference information that characterizes slope deformation.

[0056] The key steps of differential interferometry processing include: Interferogram generation: complex conjugate multiplication of the registered master and slave images to obtain the interferogram. Phase unwrapping: phase unwrapping of the interferogram to convert the phase difference into a continuous displacement. Orbit refinement and atmospheric correction: eliminate the influence of orbit error and atmospheric effect on the differential phase. Terrain removal: use the digital elevation model (DEM) to remove the terrain phase and obtain the net deformation phase.

[0057] Generate radar line-of-sight displacement time series: For each slope unit i, extract its displacement in each differential interferogram to form a displacement time series D n ={d 1n ,d 2n ,......,d in ,......,d In}. Each element d of the displacement time series in It represents the cumulative displacement of slope unit i at the nth time point relative to the main image time. The cumulative displacement refers to the total displacement from the main image time to the nth time point, reflecting the deformation of slope unit i during this time period.

[0058] Construct the displacement time series matrix of slope units: combine the displacement time series of all slope units into a matrix D n ={d 1n ,d 2n ,......,d in ,......,d In}, where I is the total number of slope units and N is the total number of time points. Matrix D n Each row represents the displacement time series of a slope unit at different time points, and each column represents the displacement distribution of all slope units at the same time point. Perform time series analysis and trend extraction: Analyze the displacement time series of each slope unit to extract characteristic information such as deformation trend and rate change. Common time series analysis methods include: Linear regression: Fit the displacement time series to obtain deformation rate and trend. Periodic analysis: Extract periodic changes in displacement time series, such as seasonal deformation. Mutation point detection: Identify mutation points in displacement time series to characterize sudden deformation events. Based on the results of time series analysis, evaluate the deformation characteristics and stability of slope units and identify potential unstable areas.

[0059] Determine the location of rain gauges: Select multiple representative locations on the target slope and its adjacent areas to place rain gauges. The following factors should be considered when placing the rain gauges: Spatial distribution: Rain gauges should be distributed as evenly as possible in the study area to capture the spatial variability of rainfall. Topographic characteristics: Select locations with relatively flat and open terrain to prevent the rain gauges from being blocked by terrain or affected by the local microenvironment. Accessibility: The location of the rain gauge should be convenient for regular maintenance and data collection. Record the precise geographic coordinates of each rain gauge so that rainfall data can be associated with the slope location. Continuously collect rainfall data from each rain gauge at set time intervals. Regularly inspect and maintain rain gauges to ensure their normal operation and data quality. Perform necessary quality control on the collected rainfall data to eliminate outliers or erroneous records. Generate a regional rainfall time series: For each sampling period n, extract the rainfall recorded by all rain gauges during that period to form a rainfall vector P n ={p 1n ,p 2n ,......,p jn ,......,p Jn}, where J is the total number of rain gauges. The rainfall vectors of all sampling periods are arranged in chronological order to generate the regional rainfall time series P n ={p 1n ,p 2n ,......,p jn ,......,p Jn}, where N is the total number of sampling periods. n Each row of represents the rainfall distribution of all rain gauges in a sampling period, and each column represents the rainfall change of a rain gauge in different sampling periods.

[0060] Figure 3 The exemplary flow chart of constructing a slope slip prediction model according to some embodiments of the present application is as follows. Grey correlation analysis is used to calculate the correlation between the slope unit displacement sequence and the regional rainfall sequence, and a slope slip prediction model based on InSAR displacement and rainfall factors is constructed. S31. Normalization of the slope unit displacement time series: The obtained slope unit displacement time series D n ={d 1n ,d 2n ,......,d in ,......,d In} is dimensionless. Calculate the minimum value min(D) and maximum value max(D) of the displacement time series. For each slope unit displacement value d i Normalize and get the normalized value, x i The formula is: i=1,2,......,I,i=1,2,......,I. The normalized slope unit displacement time series is obtained

[0061] X={x 1 ,x 2 ,......,x i ,......,x I}.

[0062] Normalization of regional rainfall time series: The obtained regional rainfall time series

[0063] P n ={p 1n ,p 2n ,......,p jn ,......,p Jn} to make it dimensionless. Calculate the minimum value min(P) and maximum value max(P) of the rainfall time series. For each time point, the rainfall p j Normalize and get the normalized value y j , the formula is: j=1,2,......,J,j=1,2,......,J。 The normalized regional rainfall time series Y={y 1 ,y 2 ,......,y i ,......,y J}.

[0064] Calculate the grey correlation coefficient: Calculate the grey correlation coefficient ξ(i, j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall. Select a suitable resolution coefficient ρ, with a value range of (0, 1). For each pair of slope unit i and time point j, calculate the grey correlation coefficient ξ(i, j), the formula is:

[0065]

[0066] Calculate the grey correlation degree: Calculate the grey correlation degree r(i, j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall. For each pair of slope unit i and time point j, calculate the grey correlation degree r(i, j), the formula is: Where N is the sequence length.

[0067] Screening rainfall sensitive units: Set the correlation threshold rt and select an appropriate threshold according to the actual situation. The value range of the grey correlation r(i,j) is [0, 1], where 0 means that the two sequences are completely unrelated and 1 means that the two sequences are completely related. tThe slope units i are screened as slope units with greater correlation with regional rainfall. The screened slope units constitute the rainfall sensitive unit set R s . The value range of the correlation threshold is given.

[0068] Construct a slope slip prediction model based on InSAR displacement and rainfall factor: extract the rainfall sensitive unit set R from the acquired slope unit displacement time series data. s The displacement time series data D' of each slope unit in the obtained regional rainfall time series data are extracted from the rainfall sensitive unit set R s Corresponding regional rainfall time series data P'. Ensure that the time points of the extracted displacement time series data D' and rainfall time series data P' correspond one to one. Construct a training data set based on the input form of the slope slip prediction model. For each time point t, construct an input feature vector x(t), including the displacement of the slope unit at the previous m time points and the regional rainfall at the previous n time points, that is: x(t) = [d(t-1), d(t-2), ..., d(tm); p(t), p(t-1), ..., p(tn)], and the corresponding target output value is the actual displacement d(t) of the slope unit at time point t. Combine the input feature vectors and target output values ​​of all time points into a training data set. Input the training data set into the support vector machine SVM model for training.

[0069] Select Gaussian radial basis kernel function RBF as the kernel function of SVM, the formula is:

[0070] Solve the support vector x of SVM through optimization algorithm (such as SMO algorithm) l , weight coefficient w l and bias term b to minimize the training error and model complexity. Adjust the hyperparameters of SVM, such as the penalty coefficient C and the kernel function parameter σ, to obtain the best model performance. Methods such as cross-validation can be used to select the best hyperparameter combination. Use a test data set that is not involved in training to test and evaluate the trained SVM model. For each time point t in the test data set, calculate the predicted displacement d'(t) based on the input feature vector x(t), using the formula: Compare the predicted displacement d'(t) with the actual displacement d(t), and calculate the prediction error and evaluation indicators of the model, such as root mean square error (RMSE), mean absolute error (MAE), etc. Analyze the prediction performance of the model and evaluate the generalization ability and applicability of the model. For the future time point t, construct the input feature vector x(t) based on the known slope unit displacement and regional rainfall. Input the input feature vector x(t) into the trained SVM model and calculate the predicted displacement d'(t). The predicted displacement d'(t) is used as the displacement prediction value of the slope unit at the future time point t for early warning and prevention and control decisions of slope slip.

[0071] Figure 4 is an exemplary flow chart of obtaining early warning analysis slope according to some embodiments of the present application, performing trend analysis and periodicity analysis on the displacement time series data of the slope unit, extracting the displacement time series data: from the screened rainfall sensitive unit set R s In the process of extracting the displacement time series data, each slope unit is visited one by one. For each slope unit, the corresponding displacement time series data is extracted, including: timestamp: recording the acquisition time of the displacement data. Displacement value: recording the displacement of the slope unit at the corresponding timestamp. The extracted displacement time series data are organized into a data set, and each data point contains a timestamp and a corresponding displacement value. The extracted displacement time series data are preprocessed as necessary, such as removing outliers, interpolating missing data, etc. The preprocessed displacement time series data constitutes the displacement time series data set D to be analyzed. s .

[0072] Local polynomial regression fitting: Set the initial smoothing parameter a and the polynomial order q, which control the smoothing degree and fitting order of the local polynomial regression. Define the Gaussian weight function W, which is used to calculate the weights between time points. The form of the Gaussian weight function is: Among them, Δt is the time interval between time points, and a is the smoothing parameter. s At each time point t i , perform the following operations: Calculate the time point t j Relative to t i The time interval Δt ij According to the Gaussian weight function W and the time interval Δt ij , calculate the time point t j Relative to t i The weight value w ij . All time points t j Relative to t i The weight value w ij Composed weight matrix W iUsing the weighted least squares method, at time point t i Based on the weight matrix W i For the displacement time series dataset D s Perform q-order local polynomial regression. The regression model is in the form of:

[0073] f(t)=β 0 +β 1 (tt i )+β 2 (tt i ) 2 +......+β q (tt i ) q ; where f(t) is the displacement value, t is the timestamp, β 0 ,β 1 ,......,β q is the regression coefficient. By minimizing the weighted residual sum of squares, the regression coefficient β is solved 0 ,β 1 ,......,β q , get the time point t i The local polynomial regression model M at i Based on the local polynomial regression model M i , calculate the time point t i The fitted value f at i and the fitted residual e i . Traverse the displacement time series dataset D s Repeat for each time point in to obtain the fitted values ​​and fitted residuals for all time points.

[0074] The optimal smoothing parameter a and polynomial order q are obtained through cross-validation. Data division: The displacement time series data set D s Divide into k mutually exclusive subsets (usually k=5 or k=10), denoted by D 1 ,D 2 ,......,D k. Each subset contains a roughly equal number of data points and is as evenly distributed in time as possible. Parameter grid setting: Set the value range and step size of the smoothing parameter a and the polynomial order q to form a parameter grid. For example, the value range of a can be [0.1, 0.5, 1, 2, 5], and the value range of q can be [1, 2, 3, 4]. Cross-validation loop: For each parameter combination (a, q) in the parameter grid, perform the following steps: Divide the data set Ds into k subsets, select k-1 subsets as training sets each time, and the remaining 1 subset as a validation set. Perform local polynomial regression fitting on the training set to obtain the fitting model M(a, q). Evaluate the performance of the fitting model M(a, q) on the validation set, and calculate the mean square error MSE(a, q) and the coefficient of determination R-squared(a, q). Repeat k times, each time selecting a different subset as the validation set, and obtain k performance index values. Calculate the average value of the k performance indicators as the comprehensive performance evaluation result of the parameter combination (a, q). Optimal parameter selection: Compare the comprehensive performance evaluation results of all parameter combinations (a, q). Select the parameter combination with the smallest mean square error (MSE) or the largest coefficient of determination (R-squared) as the optimal parameter (a opt ,q opt ). Model refitting: Use the optimal parameters (a opt ,q opt ), in the entire displacement time series data set D s Re-perform local polynomial regression fitting on the . The final fitting model M(a opt ,q opt ), as the basis for subsequent trend analysis. Using the optimal smoothing parameter a and the optimal polynomial order q, the entire displacement time series data set D is repeatedly s Perform local polynomial regression fitting. Obtain the final fitting values ​​of all time points to form a smooth curve f(t) of the displacement time series. All smooth curves f(t) constitute a trend feature set F for subsequent trend analysis.

[0075] Calculate the acceleration feature: For each smooth curve f(t) in the trend feature set F, perform the following operations: numerically differentiate the smooth curve f(t) and calculate the slope value at each time point. The finite difference method can be used to calculate the numerical differentiation, such as the central difference formula: Among them, f'(t i ) represents the time point t i The slope value at i+1 ) and f(t i-1 ) represent the time point t i+1 and t i-1 The displacement value at , Δt is the time interval. For the start and end time points, the one-sided difference formula can be used:

[0076] Among them, t 0 and t n Represent the start and end time points respectively. The calculated slope values ​​constitute the slope sequence f'(t). Perform differential operation on the slope sequence f'(t) to obtain the acceleration sequence f"(t): Perform first-order differential operation on the slope sequence f'(t) to calculate the slope change rate between adjacent time points. The calculation formula of the acceleration sequence f"(t) is: Among them, f”(t i ) represents the time point t i The acceleration value at i+1 ) and f'(t i ) represent the time point t i+1 and t i The slope value at , Δt is the time interval. For the starting time point, the one-sided difference formula can be used: The calculated slope sequence f'(t) and acceleration sequence f"(t) constitute the acceleration feature set A. All smooth curves in the trend feature set F are processed to obtain the complete acceleration feature set A.

[0077] Trend judgment: set the slope threshold s 1 and acceleration threshold a 1 , these thresholds are used to determine whether the displacement trend of the slope unit is abnormal. The trend of the slope sequence f'(t) and acceleration sequence f" (t) of each slope unit is determined by traversing each time point of the slope sequence f'(t) and acceleration sequence f" (t). For each time point t i , determine whether the following conditions are met at the same time: the slope sequence f'(t) at time point t i The value at is greater than the slope threshold s 1 For rock slopes, the slope threshold s 1 The value range is 1 to 5 mm / d. For soil slopes, the slope threshold s 1 The value range is 5-10 mm / d. For high-risk or important slopes, the slope threshold can be appropriately lowered, such as 0.5-2 mm / d, to improve the warning sensitivity. For low-risk or minor slopes, the slope threshold can be appropriately increased, such as 10-20 mm / d, to reduce false alarms. The acceleration sequence f”(t) at time point t i The value at is greater than the acceleration threshold a 1 For rock slopes, the acceleration threshold a 1 The value range is 0.1~0.5mm / d 2 For soil slopes, the acceleration threshold a 1 The value range is 0.5~1.0mm / d2 For high-risk or important slopes, the acceleration threshold can be appropriately lowered, such as 0.05-0.2 mm / d 2 For low-risk or minor slopes, the acceleration threshold can be appropriately increased, such as 1.0 to 2.0 mm / d. 2 , in order to reduce false alarms. If in a certain period of time [t 1 ,t 2 ], if multiple consecutive time points meet the above conditions at the same time, the time period is marked as an abnormal time period. If there is an abnormal time period, the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed for further analysis. If within the entire monitoring time range, the slope sequence f'(t) and the acceleration sequence f"(t) do not exceed their respective thresholds at the same time, the displacement trend of the slope unit is considered normal, and its displacement time series data is passed for subsequent processing. Trend judgment is performed on each slope unit to identify the early warning analysis object.

[0078] Spectral analysis is used to preprocess the displacement time series data of slope units that do not show an acceleration trend: outliers and noise in the time series data are removed to ensure the quality and continuity of the data. Detrending is performed on the time series data to eliminate linear trends or slow-changing trends in the data and stabilize it. Fourier transform is performed on the preprocessed displacement time series data: the displacement time series data in the time domain is converted into frequency features in the frequency domain. The fast Fourier transform (FFT) algorithm is used to calculate the spectrum of the displacement time series. In order to improve the calculation efficiency of FFT, the length of the displacement time series needs to be supplemented to an integer power of 2. Calculate the next integer power of 2 of the length of the displacement time series N, recorded as M. Supplement MN zero values ​​at the end of the displacement time series to make its length reach M. FFT calculation is performed on the displacement time series data after zero padding to convert the displacement time series in the time domain into frequency features in the frequency domain. The FFT algorithm efficiently calculates the discrete Fourier transform (DFT) by recursively decomposing the time series into even and odd terms and using binary inversion and butterfly operations. The time complexity of the FFT algorithm is O(N*log(N)), which is significantly more efficient than the O(N^2) of directly calculating the DFT. The result of FFT is a complex sequence, which represents the amplitude and phase information of the displacement time series at different frequencies. According to the result of FFT, the spectrum of the displacement time series is calculated. The horizontal axis of the spectrum represents the frequency, and the vertical axis represents the amplitude at the corresponding frequency. The amplitude can be obtained by calculating the modulus of the complex number in the FFT result, that is, The frequency range is where f sis the sampling frequency, that is, the inverse of the sampling interval of the time series. Analyze the spectral characteristics of the displacement time series and identify the main frequency components. Determine the frequency and amplitude corresponding to the peak in the spectrum, indicating the periodic characteristics of the displacement time series. Record the frequency characteristics of each displacement time series, including the dominant frequency, peak frequency and corresponding amplitude. The spectrum represents the amplitude and phase information of the displacement time series at different frequencies. Analyze the spectral characteristics of the displacement time series: identify the main frequency components in the spectrum and determine the dominant frequency of the displacement time series. Calculate the frequency and amplitude corresponding to the peak of the spectrum, indicating the periodic characteristics of the displacement time series. Record the frequency characteristics of each displacement time series, including the dominant frequency, peak frequency and corresponding amplitude. The frequency characteristics of all slope units constitute the frequency feature set P s , frequency feature set P s Contains the frequency characteristic information of all slope units that do not show an acceleration trend. Frequency characteristic set P s It can be used for subsequent periodic judgment and analysis.

[0079] Periodic judgment, set the frequency threshold f 1 : frequency threshold f 1 Used to determine whether the displacement time series has significant periodic characteristics. Frequency threshold f 1 The value range of can be set according to the characteristics of the slope and engineering experience, and can generally be set to 0.1 to 1.0 cycles / day. For daily periodic displacement, the frequency threshold f 1 The value can be 1.0 week / day. For seasonal displacement, the frequency threshold f 1 The value can be 0.1 to 0.3 weeks / day. s The frequency characteristics of each slope unit in the displacement time series are judged by comparing the frequency corresponding to the peak of the displacement time series spectrum with the frequency threshold f 1 If the peak frequency is greater than the frequency threshold f 1 , indicating that the displacement time series has significant periodic characteristics, the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed for further analysis. 1 , indicating that the displacement time series does not have significant periodic characteristics, the corresponding slope unit is marked as a stable object, and the process returns to step S1 to continue monitoring the target slope. According to the results of the periodicity judgment, the slope units are classified and processed: for the slope units marked as early warning analysis objects, further early warning analysis and evaluation are performed. For the slope units marked as stable objects, monitoring and data collection are continued.

[0080] In the slip warning stage, obtain rainfall forecast data: determine the forecast time span M, that is, the number of future time points that need to be predicted. Obtain rainfall forecast data from the N+1th time point to the N+Mth time point in the future in the area where the target slope is located. The obtained rainfall forecast data constitutes a rainfall forecast sequence

[0081] P'={p'(N+1),p'(N+2),......,p'(N+M)}, where p'(n) represents the predicted rainfall at the nth time point in the future.

[0082] Displacement prediction using slope slip prediction model: Take the rainfall forecast sequence P' as input and use the constructed slope slip prediction model based on InSAR displacement and rainfall factor to predict displacement. For the marked warning analysis object (slope unit i), predict its displacement response at the next M time points. The calculation process of displacement prediction is as follows: For the N+1th time point in the future, construct the input feature vector

[0083] x'(N+1)=[d i (N); p'(N+1)], where d i (N) is the actual displacement of slope unit i at the Nth time point, and p'(N+1) is the predicted rainfall at the N+1th time point in the future. For the nth time point in the future (n=N+2, N+3, ......, N+M), construct the input feature vector

[0084] x'(n)=[d i (N); d i '(n-1);p'(n)], where d i (N) is the actual displacement of slope unit i at the Nth time point, d i '(n-1) is the predicted displacement value at the n-1th time point in the future, and p'(n) is the predicted rainfall at the nth time point in the future.

[0085] Calculate the displacement prediction value d using the slope slip prediction model i '(n), the calculation formula is:

[0086] Among them, N s is the number of support vectors, w j is the weight of the support vector, K[·,·] is the kernel function, b is the bias term, and ε is the error term. The calculated displacement prediction values ​​constitute the displacement prediction sequence D i '={d i '(N+1),d i '(N+2),......,d i'(N+M)}, represents the displacement prediction result of slope unit i at the next M time points. Extracting key monitoring objects: From the displacement prediction sequence D i 'Extract the predicted displacement value d of slope unit i at the N+1th time point in the future i '(N+1). i '(N+1) is the key monitoring object for early warning analysis and is used for subsequent threshold comparison and early warning level judgment.

[0087] Slip warning level judgment: Set the displacement prediction value threshold to judge the slip warning level. Blue warning (mild slip): The displacement prediction value exceeds 5mm but less than 10mm. Yellow warning (moderate slip): The displacement prediction value exceeds 10mm but less than 20mm. Orange warning (severe slip): The displacement prediction value exceeds 20mm but less than 50mm. Red warning (severe slip): The displacement prediction value exceeds 50mm. Set the displacement prediction value d of slope unit i at the future N+1 time point i '(N+1) is compared with the threshold: if d i If '(N+1) exceeds the threshold, the corresponding slope unit slip warning level is given according to the size of the displacement prediction value, which can be divided into multiple levels, such as blue warning (mild slip), yellow warning (moderate slip), orange warning (severe slip), red warning (serious slip), etc. i If (N+1) does not exceed the threshold, it is considered that slope unit i has no significant sliding risk at the future N+1 time point, and returns to continue monitoring and data collection. According to the warning level, take corresponding disaster prevention measures, such as increasing monitoring frequency, launching emergency plans, evacuating dangerous areas, etc., to reduce the impact of sliding disasters.

Claims

1. A slope stability prediction method based on deformation monitoring, comprising: S1, use the synthetic aperture radar interferometry technology InSAR to monitor the target slope and obtain the time series deformation field of the target slope; Through differential interferometry processing, a time series D representing the radar line-of-sight displacement of the slope surface is generated. n ={d 1n ,d 2n ,......,d in ,......,d In }, where d in is the cumulative displacement of slope unit i at the nth time point; S2, deploy multiple rain gauges on the target slope and adjacent areas to collect rainfall data; generate a time series P representing the regional rainfall conditions n ={p 1n ,p 2n ,......,p jn ,......,p Jn }, where p jn is the rainfall recorded by rain gauge j in the nth sampling period; S3, using grey correlation analysis, calculate the displacement sequence D of the slope unit n With regional rainfall series P n According to the correlation degree, the slope units with a correlation degree greater than the threshold value are selected to construct a slope slip prediction model based on InSAR displacement and rainfall factor. S4, extract the displacement time series data of the slope unit screened out in step S3, use the trend analysis method to determine whether the time series has an acceleration trend, and use the periodic analysis method to determine whether the displacement time series has fluctuations; if the displacement time series has an acceleration trend or fluctuations, mark the corresponding slope unit as an early warning analysis object and enter step S5; otherwise, return to step S1 to continue monitoring the target slope; S5, obtain rainfall forecast data for a period of time in the future, input the rainfall forecast data into the slope slip prediction model constructed in step S3, obtain the displacement response prediction result of the slope unit marked in step S4 under the predicted rainfall scenario; calculate the displacement prediction value d of the marked slope unit at the N+1th time point in the future i '(N+1); S6, the displacement prediction value d of the slope unit at the N+1th time point in the future obtained in step S5 is i '(N+1) is compared with the threshold value. If the displacement prediction value d i '(N+1) exceeds the threshold, then according to the displacement prediction value d i '(N+1) gives the slip warning level of the corresponding slope unit, otherwise returns to step S1.

2. The slope stability prediction method based on deformation monitoring according to claim 1 is characterized by: S1, use synthetic aperture radar interferometry InSAR technology to monitor the target slope and obtain the time series deformation field of the target slope, including: Acquire at least two synthetic aperture radar complex images of the target slope at different times; The synthetic aperture radar complex image with the earliest acquisition time is used as the main image, and the remaining synthetic aperture radar complex images acquired later than the main image are used as auxiliary images; the auxiliary images are aligned in the spatial coordinate system with the main image as a reference; The registered main image and each auxiliary image are cross-correlated and a coherence coefficient γ representing the coherence degree between the main image and each auxiliary image is obtained; the main and auxiliary image pairs are selected according to the coherence coefficient γ to generate an interference complex map; Phase unwrapping is performed on the interferometric complex image to obtain the radar line-of-sight displacement Δr of the target slope in the corresponding period of the primary and secondary images; the radar line-of-sight displacement Δr is converted into a vertical displacement Δh; Repeat the phase unwrapping steps to obtain the vertical displacement map of the target slope in different time periods; according to the time interval of the primary and auxiliary image pairs, convert the vertical displacement map in each time period into a displacement time series with time as the independent variable; The displacement time series is registered with the geographic coordinates of the target slope to generate a time series deformation field.

3. The slope stability prediction method based on deformation monitoring according to claim 2 is characterized by: Generates complex interference plots, including: The mutual correlation coefficient between the registered main image complex value S1 and each auxiliary image complex value S2 is calculated as the coherence coefficient γ. The calculation formula of the coherence coefficient γ is: Among them, E[*] represents mathematical expectation, |*| represents modular operation, and * represents complex conjugate operation; Set the coherence coefficient threshold γ t ; When the coherence coefficient γ is greater than the threshold γ t When , the corresponding main and auxiliary images are taken as the main and auxiliary image pairs; For each pair of primary and secondary images, the corresponding interference complex map is calculated. The calculation formula of the interference complex map is: I=S1×S2* Where I represents the complex interference pattern, S1 represents the complex value of the primary image, and S2* represents the complex conjugate value of the secondary image; The calculated interference complex images are integrated to obtain the interference complex atlas of the main image and the auxiliary image.

4. The slope stability prediction method based on deformation monitoring according to claim 3 is characterized by: Phase unwrapping is performed on the interference complex image to obtain the radar line-of-sight displacement of the target slope in the corresponding period of the primary and auxiliary images; the radar line-of-sight displacement is converted into vertical displacement, including: The minimum cost flow algorithm is used to perform phase unwrapping on each complex interference image in the complex interference image set; the absolute phase value Φ of the complex interference image is estimated; According to the absolute phase value Φ of the interference complex image, the radar line of sight displacement Δr of the target slope in the corresponding period of the main and auxiliary images is calculated: Where, λ is the wavelength of radar signal; Using the digital elevation model DEM data, the radar line of sight displacement Δr is converted into the vertical displacement Δh: Among them, θ is the radar incident angle, α is the surface slope angle; The vertical displacement Δh is geocoded to generate a vertical displacement map within the corresponding time period of the primary and auxiliary images.

5. The slope stability prediction method based on deformation monitoring according to claim 4 is characterized in that: Estimate the absolute phase value of the complex interferogram, including: The Goldstein filtering algorithm is used to perform spectrum filtering on the interference complex image to obtain the interference complex image I'; According to the filtered interference complex image I', the phase gradient between adjacent pixels is calculated and the cost function is constructed: f(m,n)=min(Φ(m,n)-2kπ) Where Φ(m,n) is the phase difference between pixel (m,n) and its adjacent pixel, and k is an integer; The phase unwrapping problem is transformed into a minimum cost flow problem. The pixels of the interference complex graph I' are taken as nodes, the cost function value f(m,n) is taken as the capacity of the edge, and the minimum cost maximum flow is solved to obtain the global optimal phase unwrapping result φ(m,n); According to the phase unwrapping result φ(m,n), the phase difference Φ(m,n) is calculated using the following formula: Φ(m,n)=φ(m,n)+2kπ Wherein, k is an integer.

6. The method for predicting slope stability based on deformation monitoring according to any one of claims 2 to 5, characterized in that: Construct a slope slip prediction model based on InSAR displacement and rainfall factor, including: The slope unit displacement time series D obtained in step S1 n ={d 1n ,d 2n ,......,d in ,......,d In } is dimensionless, and the normalized sequence of slope unit displacement time series is obtained as X = {x1, x2, ..., x i ,......,x I },in: The regional rainfall time series P obtained in step S2 n ={p 1n ,p 2n ,......,p jn ,......,p Jn } is dimensionless, and the normalized sequence of regional rainfall time series is obtained as Y = {y1, y2, ..., y i ,......,y J },in: Calculate the grey correlation coefficient ξ(i, j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall. The calculation formula is: Among them, ρ is the resolution coefficient, and its value range is (0, 1); Calculate the grey correlation degree r(i, j) between the normalized sequence X of slope unit displacement and the normalized sequence Y of regional rainfall. The calculation formula is: Where i represents the index time node of the slope unit displacement time series; j represents the index time node of the regional rainfall time series; k represents the comparison index used in the calculation of grey correlation; N is the sequence length; Set the correlation threshold r t , the grey relational degree r(i,j) is greater than the threshold r t The slope units are screened as those with greater correlation with regional rainfall, forming the rainfall sensitive unit set R s ; From the slope unit displacement time series data obtained in step S1 and the regional rainfall time series data obtained in step S2, the rainfall sensitive unit set R is extracted. s The displacement time series data D' and the corresponding regional rainfall time series data P' of each slope unit are extracted, and the extracted displacement time series data D' and rainfall time series data P' are input into the support vector machine SVM model to build a slope slip prediction model based on InSAR displacement and rainfall factors.

7. The slope stability prediction method based on deformation monitoring according to claim 6 is characterized by: A slope slip prediction model based on InSAR displacement and rainfall factor is constructed. The model is: Where d'(t) is the predicted displacement value of the slope unit at the tth time point in the future; L is the number of support vectors; w l is the weight coefficient corresponding to the lth support vector; K[x(t),x l ] is the Gaussian radial basis kernel function RBF: x(t) is the input feature vector at the tth time point in the future, including the displacement of the slope unit at the previous m time points and the regional rainfall at the previous n time points, that is: x(t)=[d(t-1),d(t-2),...,d(tm);p(t),p(t-1),...,p(tn)] x l is the input feature vector corresponding to the lth support vector; ‖x(t)-x l ‖ represents the Euclidean distance; σ is the parameter of the RBF kernel function, which controls the width of the kernel function; b is the bias term; ε is the model error term.

8. The slope stability prediction method based on deformation monitoring according to claim 6 is characterized by: S4, extracting the displacement time series data of the slope unit screened out in step S3, using a trend analysis method to determine whether the time series has an acceleration trend, and using a periodic analysis method to determine whether the displacement time series has fluctuations, including: The rainfall sensitive unit set R selected in step S3 s In the above example, the displacement time series data of each slope unit are extracted to form the displacement time series data set D to be analyzed. s ; The local polynomial regression method is used to train the dataset D s Fit each displacement time series data in to obtain the smooth curve f(t) corresponding to each displacement time series, and form a trend feature set F with all smooth curves f(t); Calculate the slope value of each smooth curve f(t) in the trend feature set F at each time point to obtain the slope sequence f'(t), and perform differential operation on the slope sequence f'(t) to obtain the acceleration sequence f"(t). All slope sequences f'(t) and acceleration sequences f"(t) constitute the acceleration feature set A; According to the preset slope threshold s1 and acceleration threshold a1, the slope sequence f'(t) and acceleration sequence f"(t) of each slope unit are trend judged: If within a certain time period [t1, t2], the slope sequence f'(t) is greater than the slope threshold s1, and the acceleration sequence f"(t) is greater than the acceleration threshold a1, the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed; If the slope sequence f'(t) and the acceleration sequence f"(t) do not exceed their respective thresholds simultaneously within the entire monitoring time range, the displacement time series data of the slope unit will be passed to the subsequent spectrum analysis; The Fourier transform method is used to perform spectrum analysis on the displacement time series data of the slope units that do not show an acceleration trend, and the frequency characteristics of each displacement time series are obtained. All frequency characteristics are combined into a frequency feature set P s ; According to the preset frequency threshold f1, the frequency feature set P s The frequency characteristics of each slope unit in are judged: If the frequency corresponding to the peak value of the displacement time series spectrum is greater than the frequency threshold f1, the corresponding slope unit is marked as an early warning analysis object, and step S5 is executed; If the frequency corresponding to the peak value of the displacement time series spectrum is less than or equal to the frequency threshold f1, the corresponding slope unit is marked as a stable object, and the process returns to step S1 to continue monitoring the target slope.

9. The method for predicting slope stability based on deformation monitoring according to claim 8 is characterized in that: The local polynomial regression method is used to train the dataset D s Each displacement time series data in is fitted, including: Set the initial smoothing parameter a and the polynomial order q, and define the Gaussian weight function W; For the displacement time series dataset D s At each time point t i , calculate the time point t j Relative to t i The time interval Δt ij , according to the Gaussian weight function W and the time interval Δt ij Calculate time point t j Relative to t i The weight value w ij , get the time point t i The weight matrix W i ; Using the weighted least squares method, at time point t i Based on the weight matrix W i For the displacement time series dataset D s Perform q-order local polynomial regression to obtain time point t i The local polynomial regression model M at i ; Based on the local polynomial regression model M i , calculate the time point t i The fitted value f at i and the fitted residual e i ; Traverse the displacement time series dataset D s For each time point in , get the fitting values ​​and fitting residuals of all time points; According to the fitting residual sequence, the smoothing parameter a and the polynomial order q are optimized by cross-validation method to obtain the optimal smoothing parameter a and the optimal polynomial order q; Using the optimal smoothing parameter a and the optimal polynomial order q, the above polynomial fitting steps are repeated to obtain the final fitting values ​​of all time points, forming the smooth curve f(t) of the displacement time series.

10. The method for predicting slope stability based on deformation monitoring according to claim 9 is characterized in that: S5, obtaining the displacement response prediction result of the slope unit marked in step S4 under the predicted rainfall scenario; calculating the displacement prediction value d of the marked slope unit at the N+1th time point in the future i '(N+1), including: Obtain the rainfall forecast data from the N+1th time point to the N+Mth time point in the future in the area where the target slope is located, and form a rainfall forecast sequence P'={p'(N+1),p'(N+2),......,p'(N+M)}; where M is the forecast time span; Taking the rainfall forecast sequence P' as input, the slope slip prediction model based on InSAR displacement and rainfall factor is constructed to obtain the displacement response prediction results of the marked warning analysis object at the next M time points to form the displacement prediction sequence D i '={d i '(N+1),d i '(N+2),......,d i '(N+M)}; where d i '(n) is the predicted displacement value of slope unit i at the nth time point in the future, and the calculation formula is: Among them, N s represents the number of support vectors, representing the number of samples used for model training; x'(n) is the input feature vector at the nth time point in the future, including the actual displacement d of slope unit i at the Nth time point i (N), the displacement prediction value d at the future n-1th time point i '(n-1) and the predicted rainfall p'(n) at the nth time point in the future, that is: x'(n)=[d i (N),d i '(n-1),p'(n)],n=N+2,N+3,......,N+M The input feature vector at the N+1th time point in the future: x'(N+1)=[d i (N),p'(N+1)] Predicting sequence D from displacement i 'Extract the predicted displacement value d of slope unit i at the N+1th time point in the future i '(N+1), as the key monitoring object for early warning analysis.

Images