Data Assimilation of PSO-Kalman Filter and InSAR/GNSS for High Precision Monitoring and Stability Evaluation of Key Dam Body

1. Introduction

The tailings reservoirs are places for storing and managing tailings or other industrial waste discharged after ore separation and processing in each mine [1]. It is an artificial debris flow with a risk of dam collapse. Once the key dam body collapses, it will cause extremely serious damage to the environment and residential facilities downstream of the tailings pond. According to statistics, the probability of tailings dam break is relatively high, with an average of two to five large tailings ponds worldwide breaking each year [2]. There are more than 8000 tailings reservoirs in China, ranking first in the world, and most of them are located in the valleys along the rivers in mountainous areas, presenting a typical distribution characteristic of the river basin. Once accidents such as leakage and dam break occur, they are easy to flow into low-lying valleys and rivers, causing secondary environmental pollution along the river. Therefore, the monitoring of the surface deformation of the tailings dam body and the stability evaluation and analysis of the tailings pond are indispensable in China.

In recent years, many scholars have carried out a lot of experimental research on the safety of tailings dam. Yang Chunhe [3] conducted a systematic study on the disaster mechanism and prevention and control methods of the key dam body of tailings ponds. His team discussed the stability of the key dam body through theoretical analysis and numerical simulation. To obtain sufficient data on the stability of the key dam body, a monitoring system that integrates different methods is needed. The representative achievements in the monitoring and analysis of the stability of tailings dam bodies are as follows: Risk monitoring of tailings dam bodies, Liu Pei et al. [4] used the deep semantic segmentation method to process GF1 satellite images, and comprehensively analyzed the stability of the tailings dam body and the hazards caused by dam failure. Yan Shiyong et al. [5] conducted wide-area monitoring of the stability of tailings dam groups by using the improved maximum likelihood estimation and the DS-InSAR technology with joint distributed targets, achieving certain results. Wang Guangjin et al.[6] introduced the structural importance coefficient using the BowTie model to analyze the influence of the water table height on the tailings dam body. Their research results can provide reliable technical support for the design and stability evaluation of the dam body. In addition to on-site monitoring, Zhao et al. conducted physical experiments, using digital image technology to simulate the changes of saturation line and dam deformation under rainfall conditions in tailings reservoir. The study showed that the height of the saturation line has a significant impact on the stability of the dam body, and is of great importance for the monitoring of key dam bodies of the tailings pond. Meanwhile, a large number of deformation data were obtained by monitoring. Hu Jun et al. [7] used the cuckoo algorithm to optimize the SVM model, which solved the problem of excessively long fitting accuracy and poor robustness of the SVM model. The final experiment showed that the fitting time of this model was reduced and the robustness was improved. In view of the current research results, many scholars have mainly conducted studies on the key dam body monitoring and stability of tailings ponds from the perspectives of single monitoring techniques, data analysis, and numerical simulation. There are few ways to assimilate data using different monitoring technologies and models. This paper proposes to adopt the collaborative monitoring of InSAR and GNSS technologies for data assimilation and prediction.

Interferometric Synthetic Aperture Radar,InSAR,a new type of microwave remote sensing technology for the Earth developed in recent years. It can meet the requirements of high-precision, large-scale, and round-the-clock monitoring at a low cost. However, due to its own technical limitations, it is susceptible to atmospheric delay errors, temporal-spatial decoherence noise, and changes in the ground environment, etc. [8,9], resulting in the resolution in the time domain not reaching very high requirements and including certain gross errors in the measurement results. It is suitable for obtaining displacement and deformation information of the surface of tailings dam bodies with large areas. The development of Global Navigation Satellite System(GNSS) has provided high-precision continuous monitoring information for surface deformation monitoring. Research institutions can obtain high-time-resolution and high-precision deformation information by setting up GNSS continuous monitoring devices at ground points in the monitoring area. However, since GNSS is a monitoring technology based on a single base station, its spatial resolution and monitoring cost are limited by the number and networking form of GNSS devices. It is quite difficult to realize monitoring in large areas like tailings dams [10]. GNSS can continuously and rapidly provide punctiform time-series deformation information with millimeter-level accuracy. However, due to the high cost of its monitoring equipment, it is difficult to set up a large number of monitoring points over a wide area, making it hard to obtain high spatial resolution dam deformation data and failing to meet the requirements of large-scale monitoring fields. In contrast, InSAR can acquire deformation fields over a wide area, but it is susceptible to factors such as long time spans, large deformation amplitudes, and atmospheric interference, resulting in the inability to obtain high time-series dam deformation data and having the defect of insufficient time resolution. The Long Short Term Memory ( LSTM ) was used to predict the fused deformation data and the height of the saturation line of the tailings dam at the same period. Based on the predicted settlement data and the height of the saturation line, a Principal Component Analysis (PCA) - Random Forest (RF) evaluation model was constructed. This paper innovatively and initialtively proposes joint evaluation of the time-series stability of the tailings dam by using the predicted settlement data and the height of the saturation line, providing reliable technical support for the long-term safe operation of the critical dam body.

2. Materials and Methods

2.1. Study Area and Data Acquisition

The research area of this paper is the Houjiachong tailings reservoir of a mining company in Anhui Province. The dam not only ensures the safety of the tailings reservoir but is also a key dam body for the municipal transportation and road network. Therefore, it is extremely necessary to conduct precise and stable monitoring of it. This mining company is the first joint-stock enterprise established by the mining industry in Anhui Province that obtained the mining rights through public bidding and established through attracting investment. It mainly engages in the mining, processing and sales of iron ore products. Houjiachong tailings reservoir is located in Longqiao Village, Huangtun Township, Lujiang County, Anhui Province. It’s a fourth-class tailings reservoir. The tailings dam in the reservoir area is the Lucheng Highway, and the north side of the tailings reservoir is the backwater pool.

According to the "Safety Technical Regulations for Tailings Dams", 6 GNSS online displacement observation piles W1-W6 were arranged at the top line of the slope outside the dam, and 8 monitoring holes for the saturation line J1-J8 were set inside the dam. The observation data from each observation point are transmitted to the internal safety online monitoring system of Houjiachong tailings reservoir in real time. The orthographic images of tailings dam and the distribution of monitoring points are shown in Figure 1.

The data used in this paper mainly consist of the following three types: (1)30 sets of GNSS longitudinal deformation monitoring data of 6 displacement monitoring points W1-W6 within the study area. The data were collected regularly from the tailings dam safety online monitoring system, with a time span from December 19, 2023 to December 15, 2024; (2) 30 sets of InSAR data of the 6 displacement monitoring points within the study area. The SAR data originated from Sentinel-1 satellite images and were processed by SBAS-InSAR, with the same time span as the GNSS data;(3) 30 sets of saturation line buried depth monitoring data of 8 saturation line monitoring points J1-J8 within the study area. The data were collected from the tailings dam safety online monitoring system, with the same time span as the GNSS data. Among them, since the InSAR data obtains the LOS deformation of the radar line of sight, the GNSS monitoring data analysis shows that the key dam body of the tailings reservoir has almost no deformation in the horizontal direction. Therefore, in this paper, the LOS-directional deformation variable has been converted into a vertical-directional deformation variable by Formula (1) [11], which is regarded as the settlement result of the study area.

$$\mathrm{\Delta }H=d_{LOS}/\cos \theta$$

(1)

In this formula: $\mathrm{\Delta }H$ represents the longitudinal deformation, ${\mathrm{d}}_{L\mathrm{os}}$ represents the deformation along the radar line of sigh and $\theta$ represents the radar incidence angle.

The data processing flow is shown in Figure 2.

Figure 2 shows the data processing flow, that is, SBAS-InSAR is used to obtain time series deformation data in the study area, and PSO-Kalman filter model is used to fuse InSAR and GNSS data to improve the accuracy of InSAR data. The fusion data and the saturation line data are used to train the LSTM neural network model and the tailings dam stability evaluation model respectively. The PCA model is used for the dimensionality reduction of the existing multi-dimensional feature data, and the results are used for the training of the RF classifier. Finally, the LSTM is used to predict the time series data of the monitoring points, and the results are input into the RF model to obtain the time series stability of the dam.

2.2. Principle of SBAS-InSAR Technology

SBAS-InSAR is a new time-series surface deformation monitoring method based on InSAR proposed in recent years, which can improve the problem of decorrelation of conventional InSAR technology in temporal baseline and improve the accuracy and time continuity of InSAR monitoring [12]. Mainly by setting a threshold for the spatio-temporal baselines, all the obtained SAR images are combined in multiple ways to generate multiple small baseline sets and obtain interferograms. Then, the least square method is used within each set, and the singular value decomposition method is used between the sets. The joint solution is used to obtain high-precision time series deformation values [13]. Figure 3. shows the data processing flow of the SBAS-InSAR technology.

Under the assumption that M SAR images of the same area are obtained and arranged according to time, one of the images is selected as the master image. Based on the interferometric combination principle, the threshold is set to form N interferograms, which satisfies the following formula [14] :

$${\displaystyle \frac{M}{2}}\le N\le {\displaystyle \frac{M(M-1)}{2}}$$

(2)

For the interference images generated at the $t_a$ and $t_b$ ($ta<tb$) moments, take the $i（i=1，2，3\dots N）$ images among them. After removing the influence of the flat-earth effect and the terrain phase, and without considering some interference phase factors such as decoherence, atmospheric refraction and DEM error in the image [15], then there is :

$$\delta {\phi }_i=\phi \left(t_b\right)-\phi \left(t_a\right)\approx {\displaystyle \frac{4\mathsf{\pi }}{\lambda }}\left[d\left(t_b\right)-d\left(t_a\right)\right]$$

(3)

In formula (3), $\lambda$represents center wavelength; $d(t_b)$,$d(t_b)$represents the Line-of-sight cumulative deformation at time $t_b$,$t_a$; Taking the cumulative deformation of the line of sight at time $t_1$ as the reference value, that is, considering $d(t_1)=0$. Then, when$i=1,2,\mathrm{3...}M$, assuming the temporal deformation is $\phi (t_i)$, which is the associated phase part, then there is $\phi \left(t_i\right)\approx {\displaystyle \frac{4\mathsf{\pi }}{\lambda }}d\left(t_i\right)$.

Suppose $\phi \left(t_i\right)$ is the unwrapped phase, after the unwrapped phase error correction and the use of far-field or control point for calibration, then that is:

$${\phi }^T=\left[\phi \left(t_2\right),\phi \left(t_3\right),\dots ,\phi \left(t_M\right)\right]$$

(4)

In formula (4), ${\phi }^T$represents the unknown phase values associated with deformation, which is the $M-1$dimensional vector, that is:

$$\delta {\phi }^T=\left[\delta {\phi }_1,\delta {\phi }_2,\dots ,\delta {\phi }_N\right]$$

(5)

In formula (5),$\delta {\phi }^T$represents an N-dimensional vector of known values.

To avoid significant discontinuities in the cumulative deformation time series and ensure that the results conform to the physical mechanism of surface deformation, piecewise linearization algorithm is used for SBAS-InSAR. If the unknown phase is replaced by the average phase rate on the time baseline between continuous image pairs, the new unknown is defined as :

$$v^{\mathrm{T}}=\left[v_2,v_3,\dots ,v_M\right]=\left[{\displaystyle \frac{{\phi }_2}{t_2-t_1}},{\displaystyle \frac{{\phi }_3-{\phi }_2}{t_3-t_2}},\dots ,{\displaystyle \frac{{\phi }_M-{\phi }_{M-1}}{t_M-t_{M-1}}}\right]$$

(6)

Substituting into formula (3) to get:

$${\displaystyle \sum _{k=I{S}_i+1}^{I{E}_i}\left(t_k-t_{k-1}\right)}{v}_k,\forall i=1,2,\dots ,N$$

(7)

Therefore, a new matrix system will eventually be formed, and its expression is as shown in formula (8):

$$Bv=\delta \phi$$

(8)

In formula (8),$B$ is an $N\times (M-1)$ sized matrix [16]. Now the ($i,k$) elements are set, when $I{S}_i+1\le k<I{E}_i$, $\forall i=1,2,\dots ,N$,$B\left(i,k\right)=t_k+1-t_k$, or $B(i,k)=0$.

SBAS-InSAR constructs a small baseline interferogram set by processing a large number of SAR data sets to reconstruct the time series deformation corresponding to each interferogram subset. However, the lack of image connection ( no public image ) between different subsets often leads to the rank defect of the solution system. To overcome this problem, the singular value decomposition is employed to reorganize the observation equation, converting the original problem into a subset parameter estimation problem [17]. Subsequently, the minimum norm criterion can be applied to solve the surface deformation rate field to obtain deformation information.

2.3. Principle of PSO-Kalman Data Fusion Model

In this paper, PSO-Kalman filter model is proposed to assimilate InSAR and GNSS data. The principle of this model is that in the particle iterative calculation process of particle swarm optimization algorithm, the particle position is taken as the initial parameters $Q$ ( process noise matrix ) and $R$ ( observation noise matrix ) in Kalman filter algorithm, and the Kalman filter process is carried out. At the same time, the fitness function MSE is calculated and judged whether it meets the accuracy requirements, or whether the calculation process reaches the maximum number of iterations, so as to output the optimal initial parameters $Q$ and $R$ in Kalman filter, and then use the optimal $Q$ and $R$ as the initial parameters. The Kalman filter calculation is realized [18,19]. The main iterative calculation process of the model is shown in Figure 4.

Particle swarm algorithm is an evolutionary simulation computing algorithm proposed based on the foraging behavior of bird or fish. Its basic idea is to find the global optimal solution by leveraging the information sharing and collaborative search among individuals (particles) within the group [20]. In the calculation process, it is considered that each particle is a potential solution in the solution space. The movement of each particle will be affected by three factors : individual inertia ( current position ), individual post best position ( pbest ) and group historical best position ( gbest ). During the search iteration process, the particle will adjust its position and velocity based on pbest and gbest, and continuously optimize them using the following formula to achieve the global optimal solution.

$$v_i(t+1)=\omega *v_i(t)+c_1*r_1*(pbes{t}_i-x_i(t))+c_2*r_2*(gbes{t}_i-x_i(t))$$

(9)

$$x_i(t+1)=x_i(t)+v_i(t+1)$$

(10)

In the above formula, $v_i$ represents the speed of particle i, $t$represents the time, $\omega$ represents the weight, $c_1、c_2$ are learning factors, $r_1、r_2$ are random numbers within the interval [0, 1], $pbes{t}_i$ and $gbes{t}_i$ are pbest与gbest for particle i relatively, $x_i$ is the position of particle i.

The core feature of the Kalman filtering algorithm lies in its ability to effectively process noisy input data and observation signals based on a linear state space model, so as to accurately obtain the real state or signal of the system [21]. It possesses advantages such as recursiveness, adaptability and efficiency [22].

Equation of state:

$$X_t=A{X}_{t-1}+B{U}_t+W_t$$

(11)

Observational Equation:

$$Y_t=H{X}_t+V_t$$

(12)

In formula (11) and (12), the states of the system at time $t$ and its previous time point t - 1 are respectively $X_t$and $X_{t-1}$. the observed value of $t$ is $Y_t$, $U_t$represents the input signal of the system at time $t$, the process noise of the system at time $t$ is $W_t$, and the observation noise is $V_t$. $A、B、H$correspond to the system state transition matrix, control matrix and observation matrix respectively. In Kalman filtering, the value at time t can be predicted and calculated by the state equation and the optimal estimation value at time $t$-1, which is called a priori estimation :

$${\widehat{X}}_{\overline{t}}=A{\widehat{X}}_{\overline{t}-1}+B{U}_{t-1}$$

(13)

The corresponding prior estimation covariance is：

$$P_{\overline{t}}=A{P}_{t-1}{A}^T+Q$$

(14)

In this formula, $Q$ represents the corresponding variance of process noisy, the optimal result can be calculated through the PSO search. In order to obtain the optimal estimate Xt at the current time $t$, it is necessary to combine the formula ( 11 ) and ( 12 ) to calculate the optimal estimate at the time $t$, that is, the posterior estimate ${\widehat{X}}_t$. The formula is :

$${\widehat{X}}_t={\widehat{X}}_{\overline{t}}+K_t(Y_t-H{\widehat{X}}_{\overline{t}})$$

(15)

In formula ( 15 ), $K_t$ is called Kalman gain. After calculating the optimal estimation value at time t, it is necessary to use the following formula ( 16 ) to calculate and update the Kalman gain $K_t$, and use the following formula ( 17 ) to calculate the updated posteriori estimation covariance Pt to complete the iterative operation of the Kalman algorithm.

$$K_t=P_{\overline{t}}-H^T{(H{P}_{\overline{t}}{H}^T+R)}^{-1}$$

(16)

$$P_t=(I-K_{t}H)\ast P_{\overline{t}}$$

(17)

In the above formula ( 16 ),$R$ represents the observation noise matrix, which can be optimally calculated through the PSO search.

The advantage of particle swarm optimization algorithm to optimize the Kalman filter is that the parameter determination problem of the Kalman filter is transformed into a black box optimization problem. In the PSO calculation search process [23], a fitness function that can accurately reflect the filtering performance is used as a constraint, so that PSO can automatically and intelligently find the initial system noise parameters $Q$ and $R$ that make the Kalman filter achieve the best effect in a large number of iterations.

2.4. Data Prediction Theory

Based on the previous research, the long-term and short-term memory recurrent neural network is used to predict the future data of assimilation data and saturation line height data, which provides reliable data support for the stability analysis of key dams. Long Short-Term Memory Recurrent Neural Network LSTM [24] is proposed by adding memory units and gate modules to the Rerrent Neural Network ( RNN ), and its neuron structure is shown in Figure 5. When processing data, LSTM can avoid the gradient disappearance or explosion of RNN network in processing long sequence data, and can better process and learn to predict time series feature data with long interval and long delay. Figure 6 illustrates the working principle of the LSTM network.

By combining Figure 5 and Figure 6, it can be seen that the computation of a single neuron in the LSTM hidden layer consists of two parts: the update and output of the cell state. There are three gate functions in each neuron : forgetting gate, input gate and output gate. The calculation formulas are as follows [25].

$$f_t=\sigma (W_f{h}_{t-1}+W_f{x}_t+b_f)$$

(18)

$$i_t=\sigma (W_i{h}_{t-1}+W_i{x}_t+b_i)$$

(19)

$$o_t=\sigma (W_o{h}_{t-1}+W_o{x}_t+b_o)$$

(20)

In formula ( 18-20 ), $h_{t-1}$ represents the hidden state of the previous moment ; $\sigma$ is sigmoid function ; $f_t,i_t,o_t$represent the calculation results of the three gate modules respectively, $W_f,W_i,W_o$represent the weight matrices of the three door modules respectively, $b_f,b_i,b_o$represent the bias terms in the three gate modules. If the candidate value vector is selected for the calculation, the LSTM will update the entire cell state by using the product of the input values and the candidate value vector. The specific calculation process is as follows:

$$C_t=\tan h(W_c{h}_{t-1}+W_c{x}_t+b_c)$$

(21)

$$C_t=f_t{C}_{t-1}+i{C}_t$$

(22)

$$h_t=o_t\tanh (C_t)$$

(23)

$$f(x)={\displaystyle \frac{1}{1+e^{-x}}}$$

(24)

$$f(x)=\tanh (x)$$

(25)

In formula ( 21-25 ), $W_c$ is the weight matrix of the input unit, $b_c$is the bias term, $\tanh （\mathrm{x}）$ is called an activation function, $o_t$ is the output value of the neuron, $h_t$ represents the hidden state at the current moment.

The above represents the main implementation process of the LSTM neural network model. In practical application, the construction process of LSTM model can be divided into three key stages : the first is to build the basic framework of land subsidence prediction and complete the configuration of relevant model parameters ; then, based on the training data set and using the back error propagation algorithm for iterative optimization ; finally, the verification model is applied to the prediction of shape variables of different feature points [26].

2.5. Dam Stability Evaluation Model Based on Principal Component Analysis and Random Forest

For the stability of the dam body, the deformation data and water table data of the key dam body of the tailings reservoir were jointly calculated and analyzed to comprehensively evaluate its stability. To meet the evaluation requirements for the stability of tailings dams based on multi-dimensional monitoring data, the principal component analysis(PCA) is used for dimensionality reduction of the longitudinal deformation data and the depth data of the saturation line of each monitoring point, and then the dimensionality-reduction data is used to train the random forest model and test the training effect of the model, so as to determine the feasibility of principal component analysis combined with random forest model ( PCA-RF ) in the stability evaluation of tailings dam.

PCA is a multivariate statistical method that converts a high-dimensional data set with correlation into an unrelated low-dimensional data set through orthogonal transformation. It can achieve effective dimensionality reduction while maximizing the important information of the original data set [27]. Suppose that there are $m\times n$ data sets with the sample size of m and feature number of n, we first use the following formula to standardize it [28].This is mainly because there may be dimensional differences between multiple related variables contained in the sample data. When the dimensions of the variables are inconsistent, the influence of the variables with higher values in the analysis process will be amplified, while the influence of the variables with lower values will be reduced. In order to eliminate the deviations caused by dimensional differences and ensure the comparability of each feature variable, it is necessary to standardize the variables so that they all have the same scale [29].

$${\overline{x}}_j={\displaystyle \frac{{\displaystyle \sum _{i=1}^m{x}_{ij}}}{m}}$$

(26)

$$s_j=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{(x_{ij}-{\overline{x}}_j)}^2}}{m-1}}}$$

(27)

$${\tilde{x}}_{ij}={\displaystyle \frac{x_{ij}-\overline{x}}{s_j}}$$

(28)

In the formula, $i=1，2，3\dots m;j=1,2,3\dots n$. The mean value and standard deviation of the jth eigenvalue are ${\overline{x}}_j$ and $s_j$ relatively. After obtaining the standardized sample matrix $X$, the covariance matrix $R$ of the sample data is calculated by using the following formula ( 29 ).

$$R={\displaystyle \frac{X^{T}X}{n-1}}$$

(29)

In the formula, $X^{\mathrm{T}}$ is the transposed matrix of $X$.

Then, the eigenvalue decomposition of $R$is carried out using the following formula (30) to obtain its eigenvalues and eigenvectors.

$$R=V\Lambda V^T$$

(30)

In the formula, $\Lambda$ is the characteristic matrix which consists of eigenvalue ${\lambda }_i$, $V$is the orthogonal matrix.

Then the account of PCA is determined by eigenvalue contribution rate, the first k (dimension) data (principal components) are selected as the new orthogonal basis.

$${\theta }_k={\displaystyle \frac{{\lambda }_i}{{\displaystyle \sum _{k=i}^m{\lambda }_k}}}$$

(31)

In formula (31), $i=1,\mathrm{2...}m$. ${\theta }_k$represents the kth orthogonal basis(principal components), ${\lambda }_k$ represents the kth eigenvalue, the original data is projected onto the selected principal component axis in the end.

$$Z=X{V}_k$$

(32)

In formula ( 32 ), $Z$ is the data after dimensionality reduction. $v_k$ is a matrix composed of k principal components.

Random forest is a classification model with excellent generalization ability and robustness, consisting of multiple decision trees [30]. The model is shown in Figure 7. In practical applications, RF first conducts random sampling with replacement in the sample space through Bootstrap sampling to generate differentiated training subsets. Then, in the process of node splitting of each decision tree, N candidate features ( N ≤ P ) are randomly selected from the original feature space P to calculate the optimal splitting point, so as to reduce the correlation between features. Finally, take the mean of all decision tree prediction results or vote to obtain the final result [31].

3. Results

3.1. Results of Data Fusion

Data assimilation is to fuse the surface deformation data obtained by different monitoring methods into similar or identical [32,33], such as InSAR data and GNSS data. In this paper, considering that the Kalman filter algorithm is difficult to obtain the characteristic parameters of noise information in the calculation process [34], PSO is used to optimize the initial better noise parameters in the calculation of Kalman filter algorithm, so that the Kalman filter algorithm can achieve the best results, so as to realize the assimilation and fusion test between GNSS and InSAR time series monitoring data at monitoring points. The comparison between the ordinary Kalman filter data fusion assimilation results and the original data is shown in Figure 8, and the comparison between the PSO-based Kalman filter data assimilation results and the original data is shown in Figure 9, the vertical direction represents vertical deformation.

In order to verify the optimization of particle swarm optimization algorithm for Kalman filter data assimilation, this paper takes the GNSS monitoring data with high data accuracy in the original data as the true observation values, and selects the mean absolute error [35] ( Mean Absolute Error, MAE ) and root mean square error ( Root Mean Square Error, RMSE ) as the accuracy evaluation index to evaluate the data accuracy of the fusion results of the used fusion method. The accuracy index is calculated as formula ( 33 ), ( 34 ).

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n|y_i-{\widehat{y}}_i}|$$

(33)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}}\ast {{\displaystyle \sum _1^n\left(y_i-{\widehat{y}}_i\right)}}^2}$$

(34)

In formula ( 33 ) and ( 34 ), $y_i$ represents raw data, ${\widehat{y}}_i$ represents the predicted data, n represents the number of data.

To evaluate the effectiveness of data fusion, for the fused data obtained from the two fusion methods, this paper takes the GNSS monitoring data as the true value and each monitoring point as a unit, calculates the MAE and RMSE of the 30 sets of InSAR observation data and the fused data, as well as the overall MAE and overall RMSE of the multi-point and multi-period fused data, in order to test the overall data fusion effect. The comparison between the two fusion data of each monitoring point and the InSAR observation data MAE is shown in Figure 10. The comparison between the two fusion data of each monitoring point and the RMSE of InSAR observation data is shown in Figure 11. The overall MAE and RMSE of the data are shown in Table 1.

As can be seen from Figure 10, Figure 11 and Table 1, by using the PSO-Kalman filtering model to assimilate GNSS data and InSAR data, the accuracy of InSAR data has been significantly improved. Using the ordinary Kalman filter to assimilate InSAR and GNSS data, the overall MAE is reduced by 39.6 % and the overall RMSE is reduced by 41.5 % compared with the InSAR observation data. The PSO-Kalman filter is used to assimilate InSAR and GNSS data. Compared with the accuracy of InSAR observation data, the overall MAE decreases by 52.7 % and the overall RMSE decreases by 53.6 %. The overall accuracy of the assimilated data has improved significantly compared to the initial InSAR data. For the overall MAE and RMSE of the data obtained by the two fusion methods, the PSO-Kalman filter has a better data assimilation effect between GNSS and InSAR. Therefore, this paper will conduct subsequent data prediction experiments based on the fused data using PSO-Kalman filtering. Since the key dam body deformation obtained in this paper is in the longitudinal deformation of millimeter level, the MAE and RMSE of the data assimilation method used in this study are relatively small. However, the feasibility of the method can also be determined. This method is also feasible when the data volume is large or the error is significant. It also provides a reliable model for small deformation data assimilation analysis.

3.2. Results of Data Prediction

Based on the 30-period time series deformation data of each displacement monitoring point obtained by PSO-Kalman assimilation and the monitoring data of 8 saturation line monitoring points in the same period, the Adam optimization algorithm and Sigmoid activation function were selected to construct a one-dimensional single-step single-feature dimension architecture LSTM network. Through multiple experiments, it was found that training iterations for 300 will be the optimal number to avoid overfitting, and the data is normalized to accelerate the convergence speed of the network. The data from periods 1 to 26 in the basic dataset were used as the training data. The LSTM neural network was trained using an alternating input-output method for one time step. The trained network is used to output the deformation data and the saturation line depth data of each monitoring point in 27-30 periods, and the root mean square error (RMSE) is selected as the accuracy evaluation standard to test the prediction accuracy of the network after training. The comparison chart between the data of the 27 th-30 th period of the LSTM network prediction and the basic data is shown in Figure 12, and the RMSE value of the four-period prediction data of each monitoring point is shown in Table 2.

As can be seen from Figure 12(a), the LSTM network model can accurately predict the data change trends of the 27-30 periods of the tailings dam monitoring points. As can be learned from Table 2, the RMSE between the longitudinal deformation prediction data and the fused data of each monitoring point from 27 to 30 periods is all no more than 3mm, with the maximum being 2.2364mm and the minimum being only 0.5308mm. From Figure 12(b) of the saturation line monitoring point depth height prediction comparison chart and Table 2, it can be seen that the LSTM network can also accurately predict the change trend of the buried depth and height of the saturation line, and the RMSE of the predicted data and the original monitoring data is not more than 1.5 mm, the maximum is 1.3973 mm, and the minimum is only 0.5836 mm.

Based on the above analysis, the data errors predicted by the LSTM network are all at millimeter-level, and the RMSE of the predicted data of each saturation line monitoring point is not more than 1.5 mm, indicating a high prediction accuracy ; the maximum RMSE of the prediction data of each displacement monitoring point is not more than 2 mm. Considering that the main reason for the inability to accurately predict deformation may be that the dam body in the study area is a highway, and the displacement monitoring points are located on the outer side of the highway. The location is close to the highway, resulting in more complex deformation changes, and a relatively large prediction error is within the normal range. Therefore, it is considered that it is feasible to use the LSTM neural network to predict the GNSS and InSAR fusion data and the saturation line buried depth data of each monitoring point in the study area. In this paper, the trained LSTM network is used to predict the data for the next period of the monitoring points, that is, the dam deformation data and the buried depth data of the saturation line on December 15, 2024. The prediction results of each displacement monitoring point are shown in Table 3, and the prediction results of each saturation line buried depth monitoring point are shown in Table 4.

Based on the predicted data in Table 3 and Table 4, this paper will carry out the next evaluation and analysis of the stability of the dam in the future.

4. Discussion

4.1. Random Forest Dam Stability Evaluation

In this paper, the PSO-Kalman filter assimilation data of each displacement monitoring point and the depth data of the saturation line of each monitoring point in Section 3.1 are used as sample data, that is, a 30 × 14 sample data matrix is generated, where 30 represents the number of time series data dates and 14 represents the depth data of the saturation line monitored by 6 displacement monitoring points and 8 saturation line monitoring points. PCA is used to process the sample data. The absolute value of the load of each feature on the principal component is shown in Figure 13, and the variance of the principal component interpretation is shown in Table 5.

From Table 5, it can be seen that the cumulative explanatory variance of the first three principal components has reached 85.8463 %, and the cumulative explanatory variance of the first six principal components has reached 95.6995 %. In view of the fact that there are only 30 sample data in this paper, which is relatively small, so the first six principal components are selected as the training data for the subsequent random forest model. Then, a random forest classification model with a decision tree parameter of 50 trees was established. The result label was set as Y, and the data contained in "Y" was 30 periods of monitoring data. If there was an alarm, the label in Y was set to ' 1 ', and if there was no alarm, the label in Y was set to ' 2 '. Then, the 30 × 6 sample data after dimensionality reduction and 70 % ( 21 groups ) data in the result label Y are used for training, and the remaining 30 % ( 9 groups ) data are used for testing. The confusion matrix of the prediction results of the trained model in the test set is shown in Figure 14.

It can be seen from Figure 14. that among the 9 groups of random forest classification model accuracy test data used in this paper, 8 groups of test results are correct, and the accuracy rate is as high as 88.89 %, which confirms the feasibility of random forest classification model for dam stability evaluation. At the same time, the ROC curve is plotted and the AUC value of the data is calculated to evaluate the performance of the trained random forest classification model. Among them, ROC is a model performance evaluation tool that can reflect the relationship between the true positive(TP) rate and the false positive(FP) rate of the model at different classification thresholds. The AUC value is the area enclosed by the ROC curve and the coordinate axis, and the value ranged in [ 0,1 ]. The ROC curve is close to the upper left corner and the AUC value is close to 1, which proves that the model has good performance. The ROC curve drawn by the model is shown in Figure 15.

As can be seen from Figure 15, the AUC value of the trained random forest model obtained through calculation is 0.8999, which proves that the trained RF model has good performance and binary classification ability. For the trained random forest model, in this paper, the data for the next period predicted by the LSTM neural network is input as the input data into the model for the purpose of conducting an analysis on the stability assessment of the tailings dam in the future period. After inputting the prediction data, the model predicts that the result label is 2, which indicates that the tailings dam under study will remain stable in the future period.

4.2. Discussion on Stability of Tailings Dam

In this paper, the evaluation and prediction of the stability of the tailings dam is mainly based on the longitudinal deformation of the tailings dam and the buried depth of the saturation line at a specific time. In the training data of the random forest model, the determination of whether the dam is stable or not is based on whether an alarm occurs in the tailings reservoir safety monitoring system during the selected time period. It can be seen from the analysis that when one or more displacement monitoring points show significant displacement, the tailings dam safety monitoring system will issue an alarm. At this time, the staff inside the dam can take emergency reinforcement and other corresponding measures to reduce the risk of dam collapse. When the data of one or more saturation line buried depth monitoring points fall below the safety threshold due to continuous heavy rainfall drainage system failure, the tailings reservoir safety monitoring system will issue an alarm. At this time, the staff should take corresponding drainage measures to reduce the height of the saturation line and reduce the possibility of dam break.

5. Conclusions

In this paper, according to the requirements of tailings dam stability monitoring, the PSO-Kalman filter is used to assimilate the acquired time-series GNSS and InSAR data, and then the LSTM neural network model is used to predict the time series and verify its feasibility. Finally, the tailings dam deformation data and saturation line data are used as data support to train the safety assessment model and predict the safety of the tailings dam.

1) The PSO-Kalman filter algorithm is used to assimilate the multi-monitoring point time series GNSS and InSAR monitoring data. Compared with the original data, the overall MAE of the fused data is reduced by 52.7 % and the overall RMSE is reduced by 53.6 %.

2) The LSTM neural network prediction model used in this paper, when the number of training samples is limited, the RMSE of the test output set of the tailings dam deformation data is less than 2.5mm, and the RMSE of the test output set of the saturation line depth is less than 1.5mm. It has a good effect on the prediction of the time series deformation data of the dam body and the depth of the saturation line.

3) In this paper, a dam stability evaluation model is proposed by combining principal component analysis dimension reduction and random forest. The model is trained by using dam deformation data and saturation line buried depth data. Based on the model prediction accuracy of 88.89 %, and the model classification performance index AUC value of 0.8999, the safety of tailings dam in the future is predicted.

The research results of this paper form a multi-technology integrated system of ' monitoring-prediction-evaluation ', which is a relatively complete safety monitoring and stability evaluation method of tailings reservoir. It can be effectively applied to the stability analysis and evaluation of tailings dam in mining area, and provide reliable data and technical support.