An Improved Spatio-Temporal Kriging Interpolation Algorithm and Its Application in Slope

Slope stability analysis based on the deformation monitoring data has been commonly used to predict and warn slope disasters. However, due to breakdown of the monitoring equipment or restrictions and interferences of internal and external factors in the area, the loss of data is unavoidable during the process of the slope monitoring. This problem can be solved by the spatio-temporal Kriging interpolation algorithm. However, the subjective factors and theoretical semi-variogram of variogram models, and many parameters estimation may lead to low interpolation precision and poor calculation efficiency. In this paper, a hybrid spatio-temporal interpolation algorithm was presented by combining the improved adaptive genetic algorithm (IAGA) with the spatio-temporal Kriging interpolation method, and it was applied to monitor the deformation of HP1 slope in Yuebao open-pit mine. The results showed that the interpolation precision of the proposed method was about 1 times higher than that of the traditional spatio-temporal Kriging method and the spatio-temporal method of Gaussian process regression method. In addition, after the interpolation analysis of the missing part of the monitoring data, the maximum cumulative displacement of all monitoring points was around 35.8 mm, and the rate of the displacement deformation didn’t exceed 2 mm/d for three consecutive days. The deformation trend was basically consistent with the actual slope. It shows that the improved spatio-temporal Kriging interpolation algorithm (ISTKIA) has certain feasibility and reliability, and could provide a new idea for the research of related problems in related fields.


I. INTRODUCTION
With increasing demand for mineral resources during social and economic development in China, mineral enterprises continue to carry out mining activities, which results in a The associate editor coordinating the review of this manuscript and approving it for publication was Guangdeng Zong . slopes, occur various accidents and disasters under the influence of natural environment, external disturbance, as well as poor enterprise management and serious shortage of safety investment and etc., which have brought huge hidden danger to life and property security [2]- [4]. In order to ensure safety production and management in mines, scholars both at home and abroad have done a lot of research, especially researches on the slope deformation monitoring [5], [6]. The deformation monitoring data are commonly used to analyze the slope stability and forecast slope disasters.
Restricted by funding, equipment, technology and observational conditions, most of the small mining enterprises adopt traditional monitoring methods to monitor slope deformation in practical application. Some mining enterprises also use advanced devices such as measuring robots or GNSS for slope monitoring [7], [8], but the number of monitoring points and obtained monitoring data are generally hard to meet the needs of scientific research and engineering application. Some monitoring points are often damaged or missed, resulting in the loss of data, and affect the effective analysis of the slope stability. At present, it has become an active research field to solve the problems of data loss and scarcity. Domestic and foreign scholars have carried out a lot of research on remote sensing data [9], temperature [10], rainfall [11], soil moisture and temperature [12], PM2.5 [13], pollutant diffusion [14], ionosphere [15], dam deformation [16] and etc. Shen et al. [17] have comprehensively analyze the interpolation processing of the missing information of remote sensing data, and summarized four main classes: 1) spatialbased methods; 2) spectral-based methods; 3) temporal-based methods; and 4) hybrid methods. Results have been achieved in the research of the missing data interpolation, and could provide an important theoretical reference and technical support for the research of this paper.
To make up for and overcome the shortage or loss of the slope monitoring data, and accurately analyze the overall deformation trend and stability of the open-pit slope, it is necessary to interpolate the existing irregular, discontinuous, and limited monitoring data, which is useful to improve the integrity and continuity of the slope monitoring data, and support later data processing, slope stability analysis, and related scientific research. At present, there are three main classes of the research on the interpolation of the deformation monitoring data. The first classes only consider the relevance in the time domain, such as Lagrange interpolation method [18], Hermite interpolation method [19], and the spline interpolation method [20] and etc. These methods make use of the correlation of the data of the damaged points or missing points collected in different time. They can realize efficient interpolation, but it is easy to produce artifacts [21]. The slope is integral and composed of numerous points. There are a certain connections between points, and they are interdependent, rather than isolated individuals. However, these methods only consider the correlation of points in different time, and do not consider the influence of other monitoring points in the adjacent area. The second classes only consider the relevance in the spatial domain, such as the inverse distance weighting method and the spatial Kriging interpolation method [22]. Such methods utilize the correlations of the local or nonlocal information in the corrupted data itself. Although these methods consider the relationship between the points in the neighborhood, they do not study particularity of individual changes and the correlation of adjacent time data. Above all, these two classes of methods are one-sided and could affect the precision of interpolation data. According to a large number of practice and research results, the changes of the slope deformation monitoring data are related to both time attributes and spatial attributes. They have strong spatiotemporal correlations [23]. In view of this, some scholars have proposed a third kind of the spatio-temporal interpolation method which considers both the temporal correlations and the spatial correlations. Liu et al. [24] adopted an uneven spatio-temporal Kriging interpolation algorithm to analyze the landslide displacement. Lenda et al. [25] applied four different interpolation algorithms to build a surface model of the slope to analyze the progressive movement of the landslide. Wang and Zhang [26] and Wang et al. [27] used the spatio-temporal Kriging interpolation of Gaussian process regression and the spatio-temporal Kriging to interpolate the monitoring data for slope deformation analysis. The results showed that the root mean square error of the spatiotemporal Kriging interpolation was 40% higher than that of the spatial Kriging interpolation, while the average error was 55% higher, fully reflecting the spatio-temporal correlation of deformation monitoring points.
The analysis of the aforementioned spatio-temporal interpolation literatures indicates that the key to the precision of the spatio-temporal Kriging interpolation is the spatio-temporal variogram model. In practical application, the variogram model which is closest to the structural characteristics is selected artificially based on the semi-variogram of the experimental samples data. There are some limitations and shortages, such as subjective factors, theoretical semivariogram function, and many estimated parameters, which may lead to low interpolation precision and poor calculation efficiency [28]. To optimize the model and improve calculation efficiency and time convergence, scholars both at home and abroad have carried out a lot of researches in the fields of spacecraft control [29], the treatment of viral mutations [30], vertical take-off and landing of helicopters [31], time stability [32] and etc., using nonlinear algorithms and models such as the genetic algorithm, Markov chain and etc. They have made abundant results, which could provide theoretical references and thoughts for our research. On the basis of fully considering spatio-temporal correlations of monitoring points, and in order to overcome such limitations, this paper aims to propose an improved spatio-temporal Kriging interpolation algorithm (ISTKIA) to realize the interpolation of the missing data and ensure the continuity and integrity of the slope monitoring data. The method will be applied in HP1 (the number of slope subjected to landslide in the study area) slope in Yuebao open-pit mine, in hope of decreasing production costs of enterprises, improving work efficiency, providing technical support and decision basis for slope disaster prevention, mitigation and relief, as well as offering new ideas for related researches.
The rest of this paper is arranged as follows. Section 2 introduces the proposed hybrid spatio-temporal interpolation algorithm based on the improved adaptive genetic algorithm (IAGA) and the spatio-temporal Kriging interpolation method. Section 3 describes the experiments undertaken to verify the reliability and feasibility of the proposed ISTKIA by comparing the precision with the traditional spatio-temporal Kriging interpolation method and the spatio-temporal interpolation algorithm of Gaussian process VOLUME 8, 2020 regression. In Section 4, the missing data in monitoring are interpolated to analyze the change trend of the deformation monitoring points and the stability of the slope by taking the HP1 slope of Yuebao open-pit mine as an example. Finally, the paper summarizes the achievements and contributions of the research in Conclusions.

II. METHODLOGY
Ordinary Kriging is a reliable estimation method based on the continuous model of the random spatial variation for data interpolation [33], with the characteristics of unbiased estimation [34] and minimum variance of estimation error [35]. Given that Kriging model uses the preset global trend model and its approximate precision may not be optimal [36], the genetic algorithm and cross-validation were used to find out the optimal model of Kriging model and avoid the identification process falling into local optimum [37], [38]. Although interpolating missing data by spatio-temporal Kriging method fully considers the spatiotemporal correlations of monitoring points, there are still issues with variogram models, parameter estimation and a large amount of time for spatiotemporal modeling [39]. To overcome such shortcomings and limitations, this paper introduced the genetic algorithm with the advantage of finding global optimum solutions and proposed a hybrid spatio-temporal interpolation algorithm combining an improved adaptive genetic algorithm (IAGA) with the spatio-temporal Kriging interpolation method. It has great theoretical values and practical significance to make up for and overcome problems of the loss of slope monitoring data, and accurately analyze the stability and deformation trends of slopes in open-pit mines. The specific construction ideas are as follows.

A. THE PRINCIPLE OF THE SPATIAL KRIGING INTERPOLATION
The common spatial Kriging method is mainly based on correlation theories of variation functions and structural analysis. It takes limited known data of monitoring points as variables, and conducts linear unbiased optimal estimation for values of unsampled points. Its basic mathematical model is expressed by Eq.1: where Z(x 0 ) is the discreet value at the unknown point x 0 ; Z(x i ) is the monitored value of the unknown point at the known point x i in the neighborhood; λ i is the kriging weight and its value could be obtained by variation function; and n is the number of the known sample points. Therefore, in order to get discreet values of unknown points, the most key issue of Kriging interpolation is to get the weight coefficient λ i by constructing the variation function. Then, the estimation value γ * (h) of the variation function γ (h) could be obtained by the following Eq.2: where γ * (h) is the estimation value of the unknown point; Z(x i ) is the monitored value at the known point x i in the neighborhood of the unknown point; Z(x i + h) is the monitored value, h apart with Z(x i ); n(h) refers to the number of data pair of the sample points of h.

B. THE PRINCIPLE OF THE SPATIO-TEMPORAL KRIGING INTERPOLATION
The spatio-temporal Kriging interpolation [40], [41] extends data from spatial to spatio-temporal domain, fully taking into account spatial and temporal relevance among monitoring points [42]. The spatio-temporal position function is: , then the spatiotemporal Kriging interpolation model can be expressed as Eq.3: is the monitored value of the known sample point i in the neighborhood of the unknown point p, n(u, t) is the number of the known sample point pairs, λ i is the Kriging weight coefficient that minimizes the interpolation error, which needs to meet minimum variance unbiased estimation, as shown in Eq.4. min In the process of the spatio-temporal Kriging interpolation, the estimated values of the unknown points should be obtained. Therefore, the estimated value γ * (h s , h t ) of the vriation function γ (h s , h t ) should be calculated according to Eq.5: , is a monitoring value; h s , h t represent the spatial separation distance, respectively; n(h s , h t ) represents the number of the sample data point pairs,with the spatiotemporal separation distances of h s , h t .
Then, a commonly-used inseparable method was used to construct the variogram model of the spatio-temporal Kriging. The model is shown in Eq.6.
In Eq.6, γ st (h s , h t ) represents the spatio-temporal variation function; γ s (h s ) represents the spatial variogram; γ t (h t ) represents the time variation function; C st (0, 0), C s (0), C t (0) are the abutment values of the corresponding variograms, respectively; k 1 , k 2 and k 3 are auxiliary parameters.

C. THE IMPROVED SPATIO-TEMPORAL KRIGING INTERPOLATION ALGORITHM BASED ON IAGA
Spatial interpolation and spatio-temporal interpolation can be used to analyze the overall deformation trend of the slope. In the process of interpolation, the selected variogram model is usually subject to subjective factors, theoretical semivariograms, and parameters estimation. Then the universality and accuracy of interpolation methods will be affected to some extent.
In the computational process, the determination of spatio-temporal variation function sequence is actually a combinatorial optimization problem whose essence is NP-hard problem in optimization theory. However, due to discontinuity of the solution vector space and the difficulty in the solution neighborhood representation, it is difficult to calculate the results using classical algorithm. Then, as one of the modern optimization algorithms, the main feature of the genetic algorithm is that it can jump out of the local optimal solution with probability 1 for the nonlinear extremum problems and find the global optimal solution. It is a heuristic algorithm. The characteristic of skipping local optimum and searching for global optimum is mainly determined by crossover and mutation behaviors. References [43]- [45] have proposed an improved self-adaptive genetic algorithm for better search capacity, convergence and stability. On the basis, this paper proposes ISTKIA based on IAGA, and its calculation flow is shown in Fig.1: The basic idea of ISTKIA based on IAGA could be explained as follows:

1) CODING MODE
Generally, chromosome coding methods include floatingpoint coding, binary coding, gray coding and symbol coding. Among them, the floating-point coding can represent wide range and high precision values. It is more convenient for genetic searching in larger space, and can better represent the time-series relationship of the Kriging interpolation VOLUME 8, 2020 variogram and long-time series of slope surface. Therefore, it was selected in this paper for representation.
Firstly, according to actual situations of the mine slope, the obtained monitoring data of the slope were grouped by period and coded according to time series by the spatio-temporal Kriging variogram. Each group contained the spatial variogram function and the time variogram function in a certain period. In addition, the historical data were divided into t periods. The specific coding mode is shown in Fig.2: In Fig.2, l m1 l m2 , the m-th chromosome set, represents the spatial variogram function code l m1 (its range is 1-i, and the coded numbers refer to the adopted spatial variogram models, such as spherical models, exponential models, etc.) and the time variogram function code l m2 (its range is 1-j, and the coded numbers refer to the adopted time variogram models) in the m-th period. It should be noted that when i and j are more than single digits, l m1 l m2 could increase the length of the chromosome to indicate its coding range according to actual situations.

2) INITIAL POPULATION
Generally speaking, the initial population should be large enough to improve the probability of the optimal solution. The larger the scale, the wider the searching range, and the longer the genetic process of offspring is. In addition, the quality of the initial total group also directly affects the accuracy of the results and the efficiency of the algorithm. Therefore, in the determination of the initial variogram, the empirical model can be used firstly. On the other hand, some classical optimization algorithms, such as the improved circle algorithm based on Hamilton, can be used to optimize the initial population and obtain a better parent.

3) GENETIC CROSS AND VARIATION OF CHROMOSOMES
The chromosome cross matching of IAGA based on Kriging compile function time series will be carried out in a suitable way. That is, the order of the population should be determined according to the fitness degree of the samples. Parents with low fitness match with each other, and parents with high fitness match with each other. Then the location of the intersection will be determined according to the logistic chaotic sequence. Finally, crossover operation will be conducted for the parents who need cross pairing according to their positions.
In practice, a crossover probability value P j which is generally between 0-1 can be set first. In general, when P j = 1, it means that crossover operations have been performed on all individuals. When P j = 0, it means that no crossover operations have been implemented. A uniformly distributed random number obeying 0-1should be generated at each time. Then if the random number is greater than P j , it is necessary to perform the corresponding operations on the corresponding individual genomic set point.
Similar to the crossover operation, the mutated individual can be selected from the individuals after the crossover operation. A crossover probability value N b which is generally between 0-1 should be preset. The sample capacity is assumed to be N , and the number of variants is N b = N * b N .

4) ADAPTIVE OPTIMIZATION OF CROSSOVER PROBABILITY AND MUTATION PROBABILITY
In order to improve the quality of the offspring optimal solution and reduce the probability of the local optimal solution, the genetic mutation probability can be adaptively corrected to avoid falling into local extremum during the search in the large spatial range of the variogram time series [54]. Therefore, Eq.7 and Eq.8 were used to correct it.
In Eq.7 and Eq.8, p j 0 and p b 0 are fixed constants and often used to indicate the initial crossover probability and the mutation probability. α and β are constants with interval [0,1], f j is the maximum fitness in the cross-parent, f b is the maximum fitness of the mutant individual, f aver is the average fitness of all individuals, and f max is the maximum fitness among all individuals.

5) DETERMINATION OF THE OBJECTIVE FUNCTION
Assuming that the time-series chromosome of the spatiotemporal Kriging interpolation variogram of the open-pit mine slope surface was composed of T periods. When the variogram type was modified during a certain period, the performance of the Kriging interpolation would change accordingly, and then another set of variogram sequences would be found to optimize the interpolation effect. Then, the minimum R MSE after interpolation in the slope surface area was taken as the target by cross-validation, and the objective function was determined as Eq.9: (9) where T is the time period of the variogram sequence, n i is the number of the known points of the slope in the i-th period, and (x ij , y ij , z ij ) is the point coordinate obtained by Kriging interpolation. (x 0 ij , y 0 ij , z 0 ij ) is the original coordinate of the point obtained by Kriging interpolation. Denominator represents the number of all known monitoring points in the time period. Therefore, the mathematical meaning of the objective function is the average value of the root mean square error, that is, the smaller the objective function, the greater the individual's fitness.  [47]. In order to analyze the stability of HP1 slope and ensure the safety production and management of the mine, 40 deformation monitoring points, numbered JC01-JC40 were set up according to the actual situations of the slope, constituting 5 deformation observation lines, as shown in Fig.4. Considering the actual situation of the slope and the safety of the monitoring point arrangement worker, some platforms had few monitoring points. In order to better analyze the slope deformation trend and improve the deformation accuracy, the monitoring points were densified around the observation line near the landslide area. Leica TM30 with high accuracy and stability was used to obtain three-dimensional coordinates of each monitoring point, and the horizontal and settlement displacement were calculated.

B. ANALYSIS OF RESULTS
In order to compare and analyze the accuracy of the interpolation more clearly, the cross-validation method was used to evaluate its accuracy and test the quality of the determined variogram model. The basic ideas are as follows.
First, the observed values Z (x i )(i = 1, 2, . . . , n) were removed from the monitoring data column, then the variogram model was established by using the remaining observation values, and the estimated value Z * (x i ) of the removed monitoring point was predicted. Meanwhile, Z * (x i ) values were returned to the original data column, and the other points were selected to repeat the above operations until the estimated values of all the monitoring points were calculated. Finally, the error d Zi between the original data and the estimated value was calculated to analyze the consistency of the data. The R MSE was obtained according to Eq.10 to test the accuracy of the interpolation. Meanwhile, the correlation degree was evaluated by the coefficient R 2 . The calculation expression is shown in Eq.11. When R 2 was approximately 1, the correlation between the two values was stronger. The flow of calculation steps is shown in Fig.5.
In Eq.11, S res and S total are the sum of the squared residuals and the total sum of squares, respectively. The calculation functions are shown in Eq.12 and Eq.13.
In Eq.13,Z is the average of the sample points. In order to verify the feasibility of ISTKIA based on IAGA, the collected deformation monitoring data were divided into three periods (December 23, 2016, May 31, 2017 and February 07, 2018). The missing data were interpolated by the spatial Kriging method, the spatio-temporal Kriging method, the spatio-temporal method of Gaussian process regression and the ISTKIA based on IAGA. Four interpolation errors were obtained as shown in Table 1 according to cross-validation. In addition, the residual distribution of the ISTKIA is shown in Fig. 6.
Therefore, according to the cross-validation results of Table 1 and the residual distribution of Fig.5, the conclusions can be drawn: 1) The ISTKIA based on IAGA overcomes limitations of the variogram model subjected to subjective factors, VOLUME 8, 2020  the theoretical semi-variogram and many estimated parameters. The interpolation accuracy is better than traditional spatio-temporal Kriging method and the spatio-temporal interpolation method of Gaussian process regression, increasing about 1 time. It shows that it has certain universality and reliability. 2) With the increasing amount of the slope monitoring data and the relatively large deformation change, the ISTKIA based on IAGA fully embodies the advantages of calculation adaptability in the optimization of the variogram model and the estimated parameters.
It has better operation efficiency and higher interpolation precision.
3) The calculation accuracy of spatio-temporal interpolation method is significantly higher than that of the spatial interpolation method, regardless of the period of time, showing that the monitoring points not only have the spatial correlation, but also have strong spatiotemporal correlation. The interpolation results, consistent with the actual situations, fit the state of the slope monitoring points well. It is beneficial to the analysis of the whole deformation trend of the slope. 90724 VOLUME 8, 2020  It shows that the proposed ISTKIA has high stability and the interpolation results are feasible, and improves the accurate analysis of the change trend of monitoring points and the slope stability.

IV. SLOPE STABILITY AND DEFORMATION TRENDS ANALYSIS OF ISTKIA BASED ON IAGA
The accurate and complete deformation monitoring data are key to the analysis of slope stability and deformation trends, which is also the premise of safety production and management in mines. To solve the problem of data missing in actual monitoring, taking the HP1 slope as the research object, this paper interpolated the missing data of each monitoring point by using ISTKIA based on IAGA to ensure the continuity and integrity. However, considering many deformation monitoring points and large amount of monitoring data in the monitoring of HP1 slope as well as the paper length, only the missing coordinate values of JC09 and JC12 were calculated, and some of them were listed in Table 2. The data in the bold is the interpolation data. On this basis, the stability and the deformation trend of the HP1 slope were analyzed.
The results could provide an important technical basis and decision support for the disaster prevention, reduction and relief of mine slopes. Generally speaking, the commonly-used method to judge slope stability is to see whether the displacement of the monitoring point is more than its set warning threshold or control VOLUME 8, 2020  threshold in the deformation monitoring analysis. When the accumulative displacement of the monitoring point is greater than its early warning threshold or control threshold, there is a high probability of slope instability. Therefore, setting a reasonable and effective early warning threshold or control threshold is the key to slope stability judgement. According to the above analyses, historical data and actual situations of the HP1 slope of the Yuebao open-pit mine, as well as relevant specifications and references [48], [49], the alarm threshold was set as 30mm, and the control threshold was set as 35mm. The rate of daily average displacement shall not exceed 2mm/d for three consecutive days. When the displacement monitoring data are less than the above parameters, the probability of slope instability is small. With a large number of deformation monitoring points, JC09 and JC12 monitoring points were taken as research objects to study the stability of the slope (as shown in Fig.7 and Fig.8).
According to the cumulative displacement-time curves of JC09 and JC12 in all directions, it can be seen that the displacement of each monitoring point shows an increasing trend in the early period (Dec. 2016 to Aug. 2017). The maximum deformation could be observed at JC12 monitoring point, and the maximum cumulative displacement is about 35.8mm, reaching the control threshold; JC09 comes second, with the maximum cumulative displacement of 31.5mm which exceeds the alarm threshold. The deformation of other monitoring points is less than the alarm threshold of 30mm, and the displacement deformation does not exceed 2mm/d for three consecutive days.
Suffering from continuous rainfall, the deformation of the monitoring points in some areas gradually increases from Dec. 2016 to Aug. 2017. There is a possibility of further slope instability, which needs primary focus, heightened vigilance and strengthened prevention and control. Under such conditions, on the premise of studying the influence of the rainfall on the slope stability [50]- [52], the corresponding prevention and control measures were taken to change physical and mechanical parameters such as friction angle, cohesion and groundwater level caused by the surface rainwater infiltration. With further slope treatments, the displacement of each monitoring point on the HP1 slope basically slightly fluctuated up and down around the maximum cumulative displacement since Sep. 2017.
To sum up, although the cumulative displacement of the monitoring points on the HP1 slope was large in the early period, most of them failed to reach the alarm threshold. The cumulative displacement of all monitoring points didn't further expand after taking the corresponding treatment measures in the mine. The slope remained stable.

V. CONCLUSIONS
The analysis of slope stability through deformation monitoring data is an essential way to forecast and warn slope disasters. However, the scarce monitoring data in practical deformation monitoring couldn't meet the needs of scientific research and engineering applications and may bring deviations in slope stability analysis, which may further lead to mistakes in slope construction schemes and slope disaster prevention and control. To overcome the problems during the calculation by traditional method, such as varogram model and parameter estimation and etc., an improved spatio-temporal Kriging interpolation algorithm (ISTKIA) was developed in this paper with full consideration of spatiotemporal correlations among monitoring points. And it was successfully used to interpolate the missing monitoring data and analyze the HP1 slope stability. The main conclusions and contributions are as follows: 1) According to the calculation results, the ISTKIA based on IAGA is more universal and reliable in the time of optimal variogram and the calculated R MSE and R 2 were significantly higher than those of the traditional spatio-temporal Kriging, increasing by about 1 time. It better overcomes limitations of traditional spatiotemporal Kriging methods, and improves the interpolation precision and calculation efficiency. 2) According to the deformation monitoring data and the interpolation results, the maximum cumulative displacement of the HP1 slope was about 35.8 mm at JC12, and the displacement deformation of all monitoring points didn't exceed 2mm/d for three consecutive days. Although the displacement reached the control threshold, all monitoring points were basically stable and weren't further enlarged through proper treatments. The deformation trend of each monitoring point was basically consistent with the actual slope. It shows that the ISTKIA based on IAGA is feasible and practical, and can be adopted to guide safety production and mines management. 3) To some extent, the ISTKIA makes up for and overcomes the missing, discontinuity and unevenness of the slope monitoring data, as well as ensures the integrity of the monitoring data and the effective analysis of the monitoring points change trend. It provides an important technical basis and decision support for the disaster prevention, reduction and relief of mine slopes, and a new idea for the research of the related problems in related fields. 4) In practical deformation monitoring, due to economic and external factors, the monitoring points are usually composed of dozens of points, which can not accurately describe the slope surface. The ISTKIA based on IAGA solves the problem of scarce and missing monitoring data of the slope surface. It can effectively analyze the whole spatio-temporal evolution law of the slope, improve work efficiency of employees and reduce the production cost of enterprises. 5) With the increasing amount of the monitoring data, the initial calculation data should be increased in the later data processing, in order to optimize the calculation model, ensure interpolation precision and improve accuracy and reliability of the slope stability analysis.