Gate Road Support Deformation Forecasting Based on Multivariate Singular Spectrum Analysis and Fuzzy Time Series

: Underground mining engineers and planners in our country are faced with extremely difﬁcult working conditions and a continuous shortage of money. Production disruptions are frequent and can sometimes last more than a week. During this time, gate road support is additionally exposed to rock stress and the result is its progressive deformation and the loss of functionality of gate roads. In such an environment, it is necessary to develop a low-cost methodology to maintain a gate road support system. For this purpose, we have developed a model consisting of two main phases. The ﬁrst phase is related to support deformation monitoring, while the second phase is related to data analysis. To record support deformations over a deﬁned time horizon we use laser scanning technology together with multivariate singular spectrum analysis to conduct data processing and forecasting. Fuzzy time series is applied to classify the intensity of displacements into several independent groups (clusters).


Introduction
Due to hard working conditions and a chronic shortage of money, the underground mining in our country is very a difficult task. Most underground mines are characterized by low mechanized mining methods and production disruptions can happen often. Such mining methods are characterized by a long production cycle, and if we add disruptions then it is obvious that the gate road support is exposed to the rock stress much longer than is the case with the application of modern mechanized methods. In such a production environment, the maintenance of the gate road stability is almost impossible and is performed in an unsystematic way. In order to improve such a difficult situation, we use a laser scanner to obtain high accuracy deformation data, and multivariate singular spectrum analysis to process them and forecast the future state of deformations. Fuzzy time series is used to classify deformation data into several independent groups with respect to the magnitude of deformation intensity. In this way, we are enabling underground mining engineers to develop a plan to support maintenance. They can plan activities related to the identification of support states in the future, time of reconstruction, and type of support that can be used to decrease the progressive deformation. An effective deformation forecast can result in a decrease of costs and delays.
There are many different approaches to obtaining future states of support deformation. Ding-Ping Xu et al. [1] compared the predictions from the rock-soil composite material model, where an analytic method, along with existing data from physical model tests, was applied to acquire results. Good consistency, both in terms of strength and failure mode, was provided when a comparison of these three types of information was done [1]. Danial Jahed Armaghani et al. [2] applied three-nonlinear forecasting methods, namely a non-linear multiple regression, an artificial neural network, and an adaptive neuro-fuzzy inference system, in order to estimate the uniaxial compressive strength of rock.
Grey system prediction models were developed by Xiaobo Xiong [3] for forecasting rock tunnel displacement considering the non-linear parameters of the deformation of the tunnel. Ping Tang [4] forecasted the working face underground pressure by applying a grey system theory and genetic algorithm. Sun Yu et al. [5] proposed a time series prediction model based on a generalized regression neural network to predict the long-term potential deformation trends of surrounding rocks in a soft rock roadway tunnel. Haiming Chen et al. [6] built an intelligent displacement back analysis network of a deep mine roadway surrounding rock that is based on MATLAB NN toolbox. Francesca Bozzano et al. [7] used terrestrial SAR interferometry to support the management of each phase of tunneling. Qian Zhang et al. [8] made a numerical simulation for advanced displacement of tunnel with weak and broken surrounding rock in order to reveal the adjusting process of displacement in the main parts of a tunnel.
Grossauer et al. [9] developed an efficient way to predict displacements caused by tunnel excavation based on the use of analytical functions that describe the movements as a function of time and the face advance. A procedure based on analytical functions has been developed to support the prediction of the tunnel performance and surface movements. Merlini and Falanesca [10] illustrated the comparison between the prognosis and the crosscheck of the support methods and the replacement solutions. Several backanalyses were carried out in order to achieve the correct validation of the interventions and an optimal understanding of the deformation behavior of the rock mass. Bao-Zhen Yao et al. [11] has presented a multi-step-ahead prediction model, based on a support vector machine, for tunnel surrounding rock displacement forecasting. This paper is divided into five sections. In the section Materials and Methods, we describe the model of forecasting based on observed data and multivariate singular spectrum analysis. We introduce the fuzzy time series concept to make clusters describing the gate road support configuration deformation for every time point over observed and forecasted periods. A computational process is performed in the Numerical Example section to represent the possibilities of model and conclusions are provided in the Results and Conclusion section. Obtained results show a high level of correlation of original and reconstructed series of the support deformations. Based on these results, it can be concluded that the model is reliable and applicable for solving real-time problems in terms of predicting gate road support deformations.

Dynamic of Support Deformations
The dynamic behavior of support is induced by the rock stress surrounding the gate road. Deformation dynamics is the study of time varying response of support under dynamic loads. These loads are primarily considered a change of rock stress intensity induced by mining activities. As a stope approaches the gate road, the rock stress intensity increases.
In order to describe deformations of a gate road support, we consider the time dependent displacements only in the vertical cross-section (planar problem). Let S(t) be a surface of the support bounded by B(t) =[B l (t),B u (t)] at the current time t, where B l (t) represents the lower edge and B u (t) represents the upper edge of the support. Without loss of generality, we equal the upper edge with the lower and transform the surface of the support into one line, i.e., we reduce a real two-dimensional support to an artificial one-dimensional support; S(t = 0) ∈ R 2 → B l (t = 0) ∈ R , (Figure 1).
The initial configuration of B l (t), t = 0 is undeformed and known (recorded by laser scanning). The motion of each point on B l (t), t = 0, from the initial to the current configuration, is completely defined by a time dependent mapping function. Since the equation of this function is unknown and it is necessary to make monitoring of support deformations at equal or nonequal time intervals. Monitoring at equal intervals is much more applicable for forecasting deformations of a gate road. It is very important to use the same coordinate system to describe deformations over the time of monitoring.
where → v (t)-the position vector of the marker m in the current support configuration.  In the xz-coordinate plane, the position of the marker m(x, z, t) is defined as follows: The gate road support domain B l (t) is divided into a finite number of segments "connected" at the marker's position. Support configuration at current time t is defined by the set B l (t) ={m t,i },i = 1, 2, . . . , M, where M is the total number of markers. Let us define a section of the initial support configuration by the set D l (t = 0) ⊂ B l (t = 0); D l (t = 0) = {m 0,1 , m 0,2 , m 0.3 }. At time t = 1, the section of initial support configuration will be transformed and defined by the set D l (t = 1) ={m 1,1 , m 1,2 , m 1,3 }. At time t = 2, the topology of the section is defined by the set D l (t = 2) ={m 2,1 , m 2,2 , m 2,3 }. The timedependent path of the section topology transformation for t = 0, 1, 2 is defined as follows: Figure 3). In conventional matrix form, the data of position of all markers obtained by scanning for 1 ≤ t ≤ N can be expressed by the time-dependent position vectors as follows: The position of each marker is represented by each row of the previous matrix, while each column represents the configuration of the gate road support over time. According to Equation (2), x and z coordinates of observed markers are defined as follows: As can be seen, the obtained data can be treated as a multichannel time series. Multivariate singular spectrum analysis is a very useful tool to process such datasets.

Gateroad Support Deformation Forecasting Algorithm
Our forecasting algorithm is based on the methodology of the multivariate singular spectrum analysis (MSSA). In the paper by Harris and Yuan [12], the basic SSA (singular spectrum analysis) algorithm was presented. Hassani and Zhigljavsky [13] showed that SSA is a powerful method for time-series analysis and forecasting. Hassani and Mahmoudvand [14] pointed out that this method can be applied both to single series and jointly for several series (MSSA). When it is a case of a two or more series, we are talking about MSSA. The application of the SSA technique was found to be very useful for time series analysis in many different fields, such as medicine, biology, genetics, finance, engineering and many others [15]. The basic concept of the singular spectrum analysis is well presented in the paper by Hassani et al. [16].
According to Equation (2), x and z coordinates of the marker m over time correspond to one channel time series. If we take into consideration the total number of markers is M, then we are faced with an M-channel time series.
Consider an M-channel time series with a series length of N i ; V (i) where N is the time of monitoring. Here, we provide the procedure for the x coordinates of the markers. The procedure for the z coordinates is completely the same.
At the first stage of the algorithm, the decomposition, includes the following two steps: embedding and singular value decomposition (SVD) [14,16].
Embedding is a mapping that transfers a one dimensional time series of coordinates is known as a Hankel matrix. Thus, the above procedure for each series separately provides M different L i × K i trajectory matrices Y (i) (i = 1, . . . , M). To form a new block Hankel matrix in a vertical form, we need to have K 1 = . . . = K M = K. The result of this step is the following block Hankel trajectory matrix: where Y V indicates that the output of the first step is a block Hankel trajectory matrix formed in a vertical form; index V means vertical.
In the second step, we perform the SVD of Note also that the structure of the matrix Y V Y T V , is as follows: The structure of the matrix Y V Y T V is similar to the variance-covariance matrix in the classical multivariate statistical analysis literature. The matrix Y (i) Y (i)T , which is used in the univariate SSA, for the series X . U V i are called factor empirical orthogonal functions and V V i are the left and right eigenvectors of the trajectory matrix, often called principal components.
The second stage of the algorithm, called reconstruction, includes the following two steps: grouping and diagonal averaging or Hankelization [14,16].
The grouping step corresponds to splitting the matrices Y V 1 , . . . , Y V Lsum into several disjointed groups and summing the matrices within each group. The split of the set of indices {1, . . . , L sum } into disjointed subsets I 1 , . . . , I k corresponds to the representation of Y V = Y I 1 + . . . + Y I k . The procedure of choosing the sets I 1 , . . . , I k is called grouping. For a given group I, the contribution of the component Y V I is measured by the share of the corresponding eigenvalues: In a simple case where we have only signal and noise components, we use two groups of indices, I 1 = {1, . . . , r} and I 2 = {r + 1, . . . , L sum } and associate the group I 1 with the signal component and the group I 2 with the noise.
The purpose of diagonal averaging is to transform the reconstructed matrixŶ V i into a Hankel matrix, which can subsequently be converted into a time series. LetỸ (h) be the approximation of Y (i) obtained from the diagonal averaging step. Ifỹ The third stage of the algorithm concerns the future positions of the markers and is based on the vertical multivariate singular spectrum analysis recurrent procedure (VMSSA-R).
Let us have M-channel series X (i) and corresponding window length Optimal values of the window length are discussed in chapter 4 of the paper by Hassani and Mahmoudvand [14].
The VMSSA-R forecasting algorithm for the h-step ahead forecast is as follows: For a fixed value of K, construct the trajectory matrix for each single series separately; construct the block trajectory matrix Y V as follows: . . .
define matrix W as follows: if the matrix I M×M −WW T −1 exists and r ≤ L sum −M, then the h-step ahead VMSSA forecasts exist and is achieved by the following formula: where M). It should be noted that Equation (12) indicates that the h-step ahead forecasts of the refined series are obtained by a multi-dimensional linear recurrent formula (LRF).

Displacement Time Series Clustering
The time displacement intensity vector of marker m is defined as follows: The universe of displacement (clustering) is defined as the union of the two following sets: where ∆V is related to the observed displacement series and ∆V (i) N i +h to the the h-step ahead forecasts of the refined displacement series. Note that clustering is performed for each set separately.
To create the clusters within a defined universe of clustering, we apply the concept of fuzzy time series. Such series, which were first proposed by Song and Chissom, ref. [17] are based on fuzzy set theory proposed by Zadeh [18,19]. The most important advantage of fuzzy time series approaches is to be able to work with a very small set of data and not to require the linearity assumption. In the process of measuring the cross-section of the underground roadways, certain measurement errors can occur (conditioned by the error of the measurement method, instrumental error, personal error by the operator of the instrument, or the error due to the working environment itself), the nature of fuzzy time series is a good choice for a credible presentation of displacements of the markers and its clustering.
Let ∆V M×N be the universe of observed displacement. Data obtained by monitoring can be clustered in different ways, for example, a cluster composed of marker displacement over time separately (∆ v 12 , ∆v 32 , ∆v 24 , . . .), a cluster composed of each row or column data of the matrix ∆V. We focus on each column data in order to see how displacement intensity of configuration of the gate road support changes over time. For that purpose, we apply the concept of fuzzy time series with multiple observations [20][21][22].
A fuzzy set A of ∆V M×N is defined as follows: where f A is the membership function of A, and f A : ∆V → [0, 1] . The symbol "+" denotes the union operator. Membership function represents the degree of membership Let ∆V (t) be the universe of displacement on which fuzzy sets f i (t), (i = 1, 2, . . .) are defined and F(t) is a collection of f 1 (t), f 2 (t), . . . , F(t) that is referred to as a fuzzy time series of ∆V (t). Here, F(t) is viewed as a linguistic variable and f i (t) represents possible linguistic values of F(t). If F(t) is caused by F(t -1) only, the relationship can be expressed as [23]. If the maximum degree of membership of F(t) belongs to A i then F(t) is considered to be A i . Hence, The algorithm of a fuzzy time series model for multiple observations is composed of the following steps [23]: Step 1: define the fuzzy sets for the fuzzy time series. We set the beginning and the end of the universe of displacement as l =min∆V M×N and u =max∆V M×N , respectively. The observed data are sorted in ascending order. The distance of any two consecutive values of displacement is calculated as follows: Now it is necessary to compute the corresponding average value d and standard deviation σ d . The intention is to eliminate the outliers from the sorted data and obtain an average distance value free of distortion. Outliers are values that are either abnormally high or abnormally low in the sorted dataset. The process of elimination of outliers is performed as follows: Since the process of elimination is completed, a revised average distance value d R is calculated for the remaining values in the sorted data set. Accordingly, the universe of displacement is also revised and defined as follows: The number of equal intervals n is given by the user and each interval is characterized by an adequate linguistic variable. The variable whose values are words or sentences in a natural or artificial language is called a linguistic variable. We use five linguistic variables n = 5 to describe displacement intensity of configuration of the gate road support; A 1 =(very small), A 2 =(small), A 3 =(medium), A 4 =(high), and A 5 =(very high). The length of interval can be calculated by the equation: Each interval is obtained as: ∆v 1 = [lR, lR + L], ∆v 2 = [lR + L, lR + 2L], . . . ,∆v n = [lR + (n − 1)L, lR + nL]. Linguistic variable A i (i = 1, 2, . . . , 5) can now be defined according to defined intervals as follows: The triangular fuzzy number is used to quantify the linguistic variable. It can be defined as a triplet (a, b, c). The corresponding membership function is defined as [24]: Step 2: determine a fuzzy observation of displacement. Triangular fuzzy number O(t) =(a(t), b(t), c(t)) is also used to represent displacement of all markers at time point In other words, for every observation (column of matrix), it is necessary to separate minimum, maximum, and average values of displacement for all markers.
Step 3: calculate the fuzzy relationship, Y(t), between each fuzzified observation and defined fuzzy time series as in Step 1.
Fuzzified observation of displacement at time point t belongs to cluster A i if and only if their intersection has the highest value of membership function; µ i (t) → max .
Each µ i (t) is calculated according to Figure 4, representing the intersection between a fuzzy observation of displacement O(t) and fuzzy time series consisting of A 1 , A 2 , . . . , A i . The calculation of the membership function value of intersection is based on the following approach: From Equations (24) and (25), we can calculate the value of µ i (t) as follows:

Numerical Example
Artificial data related to displacement of markers are generated in order to verify the validity of the proposed forecast model. A cross section of the gate road support and positions of the markers are represented by Figure 5. The coordinate system origin is located at the left lower corner of the gate road support (marker 1). The time resolution of observation is five days.
Coordinate data obtained after six observations are represented in Tables 1 and 2 and Figure 6.   Measurements for t = 0, 1, . . . , 5 are used to create a forecast model, while the measurement for t = 6 is used as a control series. Detail calculation is provided for the x coordinates. Figure 7 represents x coordinates of the marker 10 obtained by observation.
The inputs needed to perform multivariate singular spectrum analysis are: • the number of markers (series) M = 19, • the time series length needs to be equal to each marker (series) N = 6, • the window length for each marker (series) L = 3, • parameter K also needs to be equal to each marker (series) K = 4.
The trajectory matrix of marker 1 is obtained as follows: from the original series X (1) = (0.000 0.011 0.025 0.039 0.055 0.068), we create the following trajectory matrix of dimension L × K = 3 × 4.
The above procedure is performed for each series separately and provides nineteen different L × K = 3 × 4 trajectory matrices Y (i) (i = 1, 2, . . . , 19). The result of this step is the following vertical formed block Hankel trajectory matrix: . . .
The eigenvalues of Y V Y T V , arranged in decreasing order are λ V 1 = 1066.7105; λ V 2 = 0.2704;λ V 3 = 0.0135;λ V 4 = 0.0080. The values of the corresponding eigenvectors U V 1 , U V 2 , U V 3 , U V 4 are represented in Table 3.  The rank of the matrix Y V is d = 4. The contribution of the component Y V I is measured by the share of the corresponding eigenvalues (Table 4). Table 4. Share of the eigenvalues. To perform the reconstruction stage, it is necessary to calculate the principal compo-

Eigenvalue
The results of this calculation are represented in Table 5. Table 5. Values of V V i . The reconstructed trajectory matrix of all markers achieved from four eigentriples is presented in Table 6.
The reconstructed series of marker M 1 expressed in the form of principal components achieved by arithmetic averaging ofŶ V 4 are represented in Table 7 and Figure 8 separately.    Table 7. Reconstruction stage of principal components of marker M 1 . Within the grouping step, we only have signal and noise components and we use two groups of indices, I 1 = λ V 1 , λ V 2 and I 2 = λ V 3 , λ V 4 and associate the group I 1 with signal component and the group I 2 with noise. The reconstructed series of x coordinates of the marker M 1 composed of the signal components are represented in Table 8. Table 8. The reconstructed series of the marker M 1 composed of the signal components.
0.011073 0.025299 0.039269 0.054865 0.069102 Figure 9 represents the original and reconstructed series of x coordinates of the marker M 1 using only signal group (I 1 ). The reconstructed seriesX (i) N i , i = 1, 2, . . . , 19 achieved by arithmetic averaging of Y V (I 1 ) are represented in Table 9. Table 9. Reconstructed series of x coordinates of markers based on the signal components.
Marker t = 0 t = 1 t = 2 t = 3 t = 4 t = 5 The error of the reconstruction stage is calculated according to the following equation of mean absolute percentage error (MAPE): The aggregated mean absolute percentage error (AMAPE) of the reconstruction stage is calculated as follows: MAPEi and AMAPE are represented in Table 10. j , j = 1, 2, 3, 4, is represented in Table 11.  Forecasted x coordinates of the markers for h = 1 step ahead are obtained according to Equation (12) and represented in Table 13. The time displacement intensity vector of marker m for t = 0, 1, . . . , 5 and t = 6, 7, 8 is represented in Table 16. The universe of displacement is defined as follows: D =[l, u], where l =min∆V 19×19 = 0 and u =max∆V 19×19 = 0.2020. The observed and forecasted data are sorted in the fol-lowing ascending order: 0, 0, . . . , 0.0063, 0.0090, . . . ,0.1958, and 0.2020. The distance of any two consecutive values is: 0, 0, . . . , 0.0027, . . . ,0.0062. The corresponding average value d and standard deviation σ d of obtained distances are: d= 0.0011, σ d = 0.0038. The process of elimination of outliers is performed with respect to the following condition: −0.00265 ≤ d ≤ 0.00503. A revised average distance is d R = 0.000695. A revised universe of displacement is ∆V R = [−0.0007, 0.2027 ]. According to Equation (19) the length of interval is 0.0338.

Results
Finally, the obtained sequences of displacement intensity of configuration of the gate road support over observed and forecasted time is as follows: A2→A2→A2→A2→A2→A1 →A1→A1, see Figure 12.

Conclusions
Efficiency of underground mine production directly depends on functionality of the main gate roads. Due to hard geological conditions and mining activities, deformations of gate roads are very intensive and happen frequently. Having the ability to monitor and forecast deformations is recognized as a key role to maintaining the stability of gate roads, i.e., to maintain their functionality. The signal and noise of deformations are isolated by multivariate singular spectrum analysis and it is also used to forecast future deformations with respect to different combinations of obtained principal components. Clustering based on fuzzy time series enables us to see how the configuration of the gate road support has generally changed over time and defines the intervals of deformation by adequate triangular fuzzy numbers. Making decisions and plans based on the intervals is much more reliable and suitable than performing such activities based on crisp values.
The developed model enables mine planners to create an effective plan to support repairs, including duration time and costs of reparation with respect to forecasted deformations. Future research will be focused on carrying an in situ measuring of the deformations in an underground mine to quantify the quality of the method. Furthermore, the focus will be on the making of a 3-D model of gate road support deformation forecasting.

Conflicts of Interest:
The authors declare no conflict of interest.