Prediction Method for Mine Earthquake in Time Sequence Based on Clustering Analysis

Abstract

Under the background of the intelligent construction of a coal mine, how to efficiently extract effective information from the massive monitoring data of mine earthquakes, and improve prediction accuracy, is a research hotspot in the field of coal mine safety production. In view of this problem, more and more machine learning methods are being applied to the prediction on mine earthquakes. Considering that clustering analysis can enhance the correlation between microseism data, we propose a method whose main idea is to cluster microseism data before establishing the prediction model, and then train the model, so as to improve prediction accuracy. Specifically, microseism events on a working face are divided into clusters in advance by the Spatial Temporal-DBSCAN(ST-DBSCAN) algorithm, then a prediction model is established with Support Vector Regression (SVR) to predict the occurrence location and daily frequency of high-energy mine earthquake events. A set of engineering experiments were conducted in H Coal Mine, and the results show that the spatial-temporal clustering analysis of microseism events can indeed improve the prediction accuracy of machine learning methods on mine earthquakes.

1. Introduction

A mining induced earthquake refers to the seismic activity caused by surface or underground mining, which is referred to as a mine earthquake [1]. The mine earthquake is an abnormal state of instability in the stress field of the surrounding rock in the mining process. It is manifested as the vibration of coal and rock mass, caused by the sudden release of local elastic energy, which can be caused by the destruction of local ore bodies, such as spalling rib, rock burst, etc., or the sliding activity of faults. In recent years, with the increase of mining depth and the mining intensity of coal resources in China, the phenomenon of mine earthquake has begun to show up with unprecedented frequency, intensity, and complexity. Strong mine earthquakes have occurred in many mining areas across the country (the maximum magnitude is 3.1). In the Ordos mining area alone, in 2021 there were six 2.0 magnitude or above mine earthquakes [2], leading to mine earthquake becoming a sensitive topic, and even causing social panic. Strong mine earthquake activities are not only easy to cause dynamic disasters such as underground rock bursts [3], but also cause consequences such as ground shaking, collapse, and building damage [4]. Therefore, it is necessary to carry out monitoring and early warning research on the mine earthquake phenomena.

Since the 1990s, microseism monitoring technology has been widely used in the monitoring and early warning of mine earthquakes by coal mining enterprises, due to its advantages of high sensitivity and strong practicability [5]. In mining and production, the dynamic and static loads cause the deformation and instability of the coal and rock mass, or else the geological fault, fold, and other defective structures, are activated to trigger geological activities, which will release vibration waves with different frequencies and energy levels. The principle of the microseism monitoring system is to capture the relevant information of seismic waves, and determine the spatial-temporal location of the source, the strength and frequency of microseism activities, and other high-energy release activity information through calculation and inversion [6]. At present, microseism monitoring technology has made great progress in station network optimization, seismic wave pickup, inversion positioning, and waveform real-time monitoring, and has made a series of important research achievements [7]. The microseism monitoring system can identify the micro vibration wave and its dynamic change process in the whole coal mine, obtain the specific information of the time, space, and intensity of the microseism event, and directly reflect the real-time mechanical characteristics evolution of the coal mine.

With the increasingly mature and extensive application of microseism monitoring technology, it has become an effective means for the early warning of coal mine disasters, such as rock burst. It works by using a large amount of microseism monitoring data to analyze the breaking and migration laws of coal seam mining surrounding rocks, and to summarize the spatial scale, combined structure, the mechanical parameters of key stratum, and the correlation between the breaking process and high-energy microseism events. He et al. [8] used microseism monitoring technology to analyze the temporal and spatial evolution law of microseism events before and after the breaking of key stratum, laying a foundation for the prediction on rock burst using the combination of strata movement theory and microseism monitoring technology. Yuan et al. [9] analyzed the time sequence characteristics of microseism signals during the period of rock burst, and then obtained the spectrum characteristics and distribution change laws of microseism signals. Wang et al. [10] analyzed the waveform signals of microseism events and obtained the power spectrum evolution characteristics before and after the occurrence of rock burst. Xiao et al. [11] revised the maximum effective amplitude based on the attenuation characteristics of microseism waves in deep buried tunnels, taking the relative effective amplitude and the maximum effective frequency as spectrum analysis parameters. Wei et al. [12] analyzed the time-domain waveform and frequency spectrum characteristics of microseism signals when rock burst occurred in the coal mine by using the fast Fourier transform method. Peng [13] classified the microseism signal categories of coal mine working faces and analyzed the distribution law of microseism events in time and space. However, in the face of microseism monitoring signals containing a large amount of information, it is difficult to manually define and extract parameters to reflect all the features of microseism events, which makes it easy to cause a large amount of effective information to be ignored. As a result, current microseism monitoring can only monitor and reflect the microseism events that have occurred, and are occurring, and it is difficult to accurately judge the potential microseism events that may occur.

Traditionally, mine earthquake prediction usually uses geophysical methods to monitor some precursory signals, and the synthetic index method with artificially defined and extracted parameters is used to evaluate the possibility of the occurrence of high-energy mine earthquakes. The determination of the index and corresponding weight involves subjectivity and inconsistency. Machine learning can well overcome the problems caused by the synthetic index method. Its modeling process does not involve too much subjective decision-making, and is a data-driven strategy. Using the machine learning model, researchers do not need to pay attention to the weight of each index and the corresponding classification standard, they just need to know the specific value of each index, which is objective and measurable. At the same time, the use of intelligent devices such as microseism monitoring systems will generate massive real-time data carrying effective information. Empirically driven and mechanistically driven mine earthquake prediction methods are not enough to use these data, resulting in the loss of effective information. Using a data-driven method to solve the problem of mine earthquake prediction will become a breakthrough point for data-driven systems to enter the traditional engineering field, and is also a key step for mining to enter the smart mine and data mine era.

After the great development of computer hardware and the upsurge of machine learning, researchers continue to try to use data-driven methods, such as machine learning, to predict and warn of mine earthquake, achieved good results. Using the machine learning method to analyze microseism monitoring signals can get the most effective information, most of which cannot be obtained by an explicit algorithm. The machine learning method also has unique advantages in establishing the relationship between monitoring parameters and the time, location, and intensity of a mine earthquake. It can analyze the automatic monitoring signal into a high-dimensional matrix, without the need to manually determine the type of extraction parameters, and maximize the retention of signal features. Vallejos and Dong [14,15] established a logistic regression model to identify microseism events and mining blasting activities, respectively, with a higher accuracy than other analogy models. Del [16] established two neural network models to identify microseism events, one of which is used to extract signal features to construct training samples, and the other is used for event classification. In addition, machine learning methods such as random forest [17], Bayesian network [18], and support vector machine [19,20], are widely used for monitoring microseism events. However, these studies do not consider the spatial-temporal aggregation and activity correlation of microseism events, and the prediction accuracy of mine earthquake is not high; there is a lot of room for improvement. Microseisms are not isolated events. Mine earthquakes caused by a stress release usually have a similar spatial location, short time interval, and have the law of time series. If we first cluster the microseism data and then use the machine learning method to establish the prediction model, we will get better results.

Therefore, we propose to combine clustering analysis and machine learning methods to predict the high-energy mine earthquake in time sequence, including the occurrence location prediction and energy frequency prediction. In the selection of clustering algorithm, considering that microseism events include space and time information, the ST-DBSCAN algorithm [21] is selected to cluster the microseism events occurring in the working face. After clustering, there are no large number of microseism events in a single cluster, so SVR [22], which is suitable for a small sample data set, is used to predict mine earthquakes in sequence [23,24]. As a comparison experiment, the number of microseism events included in the working face is large before clustering, so Long Short-Term Memory (LSTM) [25] with a good prediction effect of time series on a large sample data set is selected to predict mine earthquake [26].

The main contributions of this paper are summarized as follows.

(1) A time sequence prediction method for mine earthquake based on ST-DBSCAN using SVR is proposed. Clustering analysis can source microseism clusters and enhance the correlation of microseism data. This data preprocessing method can improve the accuracy of subsequent machine learning prediction models.

(2) Contrast engineering experiments in the H Coal Mine are conducted to evaluate the performance of the proposed method. The results show that, compared with the classical time sequence prediction method (LSTM) on unclustered microseism data, the SVR model has a better prediction effect of mine earthquake on clustered microseism data.

The remainder of this paper is organized as follows. The related machine learning algorithms are introduced in Section 2. The proposed approach is described in Section 3. A real case study is illustrated in Section 4. Conclusions and future work are summarized in Section 5.

2. Methodology Background

2.1. Spatial Temporal-DBSCAN

The ST-DBSCAN algorithm adds time dimension on the basis of the DBSCAN algorithm to form a spatial-temporal neighborhood with spatial distance as radius and time interval as height, which is used for the clustering analysis of spatial-temporal data. As shown in Figure 1, the algorithm identifies the density of samples through three parameters, namely space radius $Eps$, time window $\mathsf{\Delta }T$ and density threshold $Minpts$. Based on this, the data meeting the conditions are divided to form the final cluster set.

Where $Eps$ is used to measure the distance between two spatial-temporal data in the spatial dimension and $\mathsf{\Delta }T$ is used to measure the distance in the time dimension. For example, $A\left(x_1,y_1,z_1,t_1\right)$ and $B\left(x_2,y_2,z_2,t_2\right)$ represent two spatial-temporal data, where $\left(x_1,y_1,z_1\right)$ and $\left(x_2,y_2,z_2\right)$ are spatial attributes and $t_1$ and $t_2$ are temporal attributes, then the calculation formulas of $Eps$ and $\mathsf{\Delta }T$ are as follows:

$$Eps=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2}$$

(1)

$$\mathsf{\Delta }T=\left|t_1-t_2\right|$$

(2)

The ST-DBSCAN algorithm uses the threshold parameter $Minpts$ to represent the minimum number of sample points required for clustering in the spatial-temporal neighborhood, which can be obtained from the total number $N$ of microseism events. The calculation formula is as follows:

$$Minpts=lnN$$

(3)

According to the relevant definitions of the ST-DBSCAN algorithm, the algorithm starts loop judgment after inputting the sample data set. If the number of adjacent points within the space radius $Eps$ and time window $\mathsf{\Delta }T$ of the core point is not less than $Minpts$, a cluster will be formed. Otherwise, the next point will be judged until all points have completed traversal and judgment. The specific construction steps are described as follows:

(1)

Algorithm inputs: sample data set $D$, spatial radius parameter $Eps$, time window parameter $\mathsf{\Delta }T$, density threshold $Minpts$.

(2)

Establish a tag list $L$, which is used to mark whether the sample point has completed the traversal state, and initially set the traversal state $L=0$ of all points in $D$.

(3)

Start traversal, randomly select point $d$ in $D$, and set its traversal state to $L=1$.

(4)

Start loop 1, if the number of points in the $Eps$ and $\mathsf{\Delta }T$ spatial-temporal neighborhood of $d$ is not less than $Minpts$, create a new cluster $C$, add $d$ to $C$, otherwise mark $d$ as a noise point.

(5)

$D^{*}$ represents the set of all points in the $Eps$ and $\mathsf{\Delta }T$ spatial-temporal neighborhood of $d$.

(6)

Start loop 2, traversing each point $d^{*}$ in $D^{*}$.

(7)

If $L=0$ of $d^{*}$, select this point and set $L=1$ of $d^{*}$.

(8)

If the number of points in the $Eps$ and $\mathsf{\Delta }T$ spatial-temporal neighborhood of $d^{*}$ is not less than $Minpts$, add these points to the set $D^{*}$.

(9)

If $d^{*}$ does not belong to any cluster, add $d^{*}$ to cluster $C$.

(10)

End loop 2, and the algorithm outputs the set of cluster $C$.

(11)

End loop 1, the termination condition of the algorithm is that the traversal states of all data points are $L=1$, that is, all data points have been processed, either classified into a cluster or marked as noise points.

2.2. Support Vector Regression

In recent years, the support vector machine (SVM) has been a hot research topic in the field of machine learning applications. It is a machine learning method based on statistical learning theory and the structural risk minimization principle. Because of its good learning ability, the SVM is widely used in many fields of life and production [27,28], mainly to solve the problems of classification and regression analysis. The SVM has a solid mathematical theoretical foundation and excellent performance, so can well solve many key problems of modeling based on neural network, especially the problem of neural network training. The SVM structure is shown in Figure 2:

The traditional SVM is a generalized linear classifier, which is used to solve data classification problem. However, the theory of the SVM can be adjusted so that it can also be applied to regression analysis problem. This improved method is called Support Vector Regression (SVR). It is flexibly applied to the modeling of nonlinear regression problem and can show unique advantage when processing multidimensional or high-dimensional data, because the nonlinear kernel function in the algorithm can map the sample data to a higher dimensional feature space and can quickly find the linear regression function. The SVR algorithm has many advantages, such as low generalization error rate, low computational complexity, not falling into local optimization, and high accuracy in dealing with nonlinear prediction and regression problem. As shown in Figure 3, the objective of the algorithm is to find the function $f\left(x\right)$, which wants to maximize the “soft edge $\mathsf{\epsilon }$” to include as many target values as possible while maintaining a certain flatness.

The regression prediction of data can be realized by finding the decision boundary that maximizes the edge. The function can be expressed as follows:

$$f\left(x\right)={\omega }^{T}x+b$$

(4)

where $\omega$ is the weighting coefficient, $b$ is the threshold, also known as the bias term, flatness means that we need to find a $\omega$ as small as possible (in the sense of Euclidean norm, that is ${\left|\left|\omega \right|\right|}^2$). This problem can be transformed into a convex optimization problem:

$$min\frac{1}{2}{\parallel \omega \parallel }^2$$

(5)

$$subject to \{\begin{array}{c}{y}_i-\omega x_i-b\le \epsilon \\ \omega x_i+b-y_i\le \epsilon \end{array}$$

(6)

where $\epsilon$ is the allowable error value of the regression function, and the smaller $\epsilon$ is, the smaller the error of the regression function will be, indicating that the sample points are more concentrated. However, not every sample point in the training set can meet the constraint conditions. Therefore, relaxation variables ${\xi }_i^{+}$ and ${\xi }_i^{-}$ can be introduced to alleviate the optimization constraint problem. The above formula can be rewritten as follows:

$$min L\left(\omega ,{\xi }_i^{+},{\xi }_i^{-}\right)=min\left(\frac{1}{2}\parallel \omega {\parallel }^2+C{{\displaystyle \sum }}_{i=1}^n\left({\xi }_i^{+}+{\xi }_i^{-}\right)\right)$$

(7)

$$subject to\{\begin{array}{c}{y}_i-\omega x_i-b\le \epsilon +{\xi }_i^{+}\\ \omega x_i+b-y_i\le \epsilon +{\xi }_i^{-}\\ {\xi }_i^{+},{\xi }_i^{-}\ge 0\end{array}$$

(8)

where $C$ is the penalty factor, which is a positive value, indicating the penalty degree of the samples exceeding the allowable error, and is used to adjust the trade-off between the relaxation variable and the boundary. Lagrange multiplier method and duality theory are applied to the above formula, and Lagrange multipliers ${\alpha }^{+}$ and ${\alpha }^{-}$ are introduced. The objective function can be changed as:

$$Max{{\displaystyle \sum }}_{i=1}^l\left[y_i\left({\alpha }_i^{+}-{\alpha }_i^{-}\right)-\epsilon \left({\alpha }_i^{+}+{\alpha }_i^{-}\right)\right]-\frac{1}{2}{{\displaystyle \sum }}_{i,j=1}^l\left({\alpha }_i^{+}-{\alpha }_i^{-}\right)\left({\alpha }_j^{+}-{\alpha }_j^{-}\right)k\left(x_i,x_j\right)$$

(9)

$$subject to\{\begin{array}{c}{\displaystyle \sum }_{i=1}^l\left({\alpha }_i^{+}-{\alpha }_i^{-}\right)=0\\ 0\le {\alpha }_i^{+},{\alpha }_i^{-}\le c\end{array}$$

(10)

where $k\left(x_i,x_j\right)$ is the kernel function, so the final fitting function can be expressed as:

$$f\left(x\right)={{\displaystyle \sum }}_{i,j=1}^l\left({\alpha }_i^{+}-{\alpha }_i^{-}\right)k\left(x_i,x_j\right)+b$$

(11)

Common kernel functions include linear kernel, polynomial kernel, sigmoid kernel, and Radial Basis Function (RBF) kernel function. Among them, RBF is the most commonly used kernel function in the SVR algorithm, because it can map sample data to high-dimensional space nonlinearly, and can handle the situation that there is a nonlinear relationship between class labels and sample attributes.

2.3. Long Short-Term Memory

In the Recurrent Neural Network (RNN), information is allowed to be transmitted in different time steps, which is conducive to learning order dependency in input data, making it more suitable for long-time series. Gradient descent is often used in the training of neural networks. However, with the increase of time step, the RNN suffers more and more from gradient disappearance and explosion, which aggravates the difficulty of model training. In order to solve this problem, the LSTM is developed based on the RNN model. The LSTM model structure is shown in Figure 4, including a core unit (memory part) and three gating structures (forget, input and output gates), which guide the information flow inside the LSTM unit.

The calculation steps of LSTM are as follows:

(1)

Determine the information to discard in the cell state. Input the output $h_{t-1}$ at the last moment and the input $x_t$ at the current moment, through the forget gate, add the Sigmoid function, and output a vector whose element values are between 0 and 1, which represents the stored proportion of the cell state $C_{t-1}$ at the last moment. The calculation formula is:

$$f_t=\sigma \left(W_f\left[h_{t-1},x_t\right]+b_f\right)$$

(12)

(2)

Determine the new information stored in the cell state. There are two steps: the input gate determines the updated value, and the candidate value vector $i_t$ is then added to the cell state. The calculation formula of the input gate value and the new candidate value vector ${\tilde{C}}_t$ are as follows:

$$i_t=\sigma \left(W_i\left[h_{t-1},x_t\right]+b_i\right)$$

(13)

$${\tilde{C}}_t=tanh\left(W_C\left[h_{t-1},x_t\right]+b_C\right)$$

(14)

(3)

Update the states of old cells. $C_{t-1}$ is updated to $C_t$ by discarding some of the information and adding new candidate values. The calculation formula of cell state value at the current moment is:

$$C_t=f_t\odot C_{t-1}+i_t\odot {\tilde{C}}_t$$

(15)

(4)

Determine the value of the output. It includes two steps: determine the value $O_t$ of the output gate, and then determine the final output result value $h_t$ at the current moment. The calculation formula is as follows:

$$O_t=\sigma \left(W_o\left[h_{t-1},x_t\right]+b_o\right)$$

(16)

$$h_t=O_t\odot \tan \mathrm{h}\left(C_t\right)$$

(17)

where $W_f$, $W_i$, $W_C$ and $W_o$ are the weight terms to be learned, $b_f$, $b_i$, $b_C$ and $b_o$ are the offset terms to be learned, $\sigma$ is the Sigmoid activation function, $\odot$ is the Hadamard product, where the components of vectors are multiplied one by one. The key advantage of using the LSTM unit is that the cell state will be updated dynamically with the passage of time, which increases the processable sequence length of time-series neural network and has better effect on prediction tasks with long-term regularity.

3. Model Design

With regard to the prediction and early warning of high-energy mine earthquake, it is difficult for the traditional experience-driven and mechanism-driven strategies to accurately answer the three key questions of when, where, and energy level [29]. This gives machine learning a chance. In fact, the data-driven method has a better solution to those questions. At present, most research follows a fixed model, that is, select a machine learning model and find the best results through parameter optimization. However, if the microseism data are directly input into the training model, the final prediction accuracy is often not high, because the microseisms are not isolated events, but have obvious spatial-temporal aggregation. Therefore, if the spatial-temporal clustering analysis of microseism events is carried out first, the closely related microseism events will be divided into clusters, so the correlation degree between microseism data will be improved. On this basis, the machine learning method is used to fit this rule and establish a prediction model, which can significantly improve the accuracy of mine earthquake prediction. Accordingly, we propose a time sequence prediction method for mine earthquake combined with clustering analysis. The specific process is shown in Figure 5, which mainly includes three parts.

(1)

The clustering of microseism events. Similar to earthquake, coal mine microseisms are in essence a kind of dynamic phenomenon caused by stress release. Generally, they do not appear as isolated events, but many microseism events with different energy levels will occur successively in a certain region and time period, to form swarm sequences. The microseism events included in a swarm sequence are caused by a stress release, so the occurrence time, location, and energy levels of these microseism events have a certain regularity. In the ST-DBSCAN algorithm, the density is obtained by calculating the number of adjacent points in the designated spatial-temporal neighborhood around the sample point. The points whose density are higher than the designated threshold will be constructed as clusters, which are very suitable for processing the spatial-temporal data recorded in the microseism monitoring system. Clustering analysis can identify the active areas of microseism events in the working face and divide the microseism events into closely related clusters. The time, location, and energy information of microseism events contained in each cluster can be regarded as time series data [30] and have a certain regularity.

(2)

The establishment of a prediction model. After the completion of cluster analysis, microseism events on the working face will be divided into clusters. Although the correlation degree between microseism events in the cluster increases, the number of microseism events in a cluster is not large. This is equivalent to dividing a large sample data set containing microseism events of the whole working face into many small sample data sets. Therefore, the SVR algorithm suitable for small samples is selected to predict the time series of microseism events within the cluster. Mine earthquake prediction includes two aspects: one is the location of high-energy events; the second is the daily frequency of high-energy events with different energy levels, that is, the information of microseism events that have occurred in a specific period of time is used to predict the possible occurrence location (X, Y and Z coordinates) and the daily frequency of subsequent high-energy events.

(3)

The comparative analysis of the model. In order to verify the effectiveness of the proposed method which uses SVR to predict mine earthquake on the basis of clustering analysis, the LSTM, a very classical method in time sequence prediction, is selected for comparison experiment. Similarly, the occurrence location (X, Y and Z coordinates) and the daily frequency of high-energy events are predicted. At the same time, the effects of the two prediction models are compared and evaluated by using the Mean Square Error (MSE) and R-square (R²). The data used by the LSTM here is the microseism data of the whole working face, without clustering analysis. On the one hand, it is to compare the effect of clustering analysis on improving the accuracy of machine learning prediction model, and on the other hand, it also meets the requirement of method for large sample size.

The proposed framework consists of the ST-DBSCAN and SVR algorithms, so its time complexity depends on the time complexity of the two algorithms. When the ST-DBSCAN runs, it needs to traverse all the points in the dataset and calculate the number of density-reachable points of each point, so its time complexity is $O\left(N^2\right)$, where N is the number of the training samples. Meanwhile, the time complexity of SVR is $O\left(N^3\right)$ [31]. Only the item with the fastest growth rate is considered, so the time complexity of proposed framework is $O\left(N^3\right)$.

4. Engineering Experiments

4.1. Microseism Data Acquisition

The data in this paper comes from the real-time microseism data collected by the SOS microseism monitoring system from Poland, specifically the microseism data of 103# working face of the H Coal Mine, with a total of 8497 records. As shown in

Table 1. Records in microseism monitoring system.

Index	Date	Time	X (m)	Y (m)	Z (m)	Energy (J)
0	7 November	19:22:47	495,987	891,476	359	57,397
1	7 November	21:09:05	495,918	891,335	366	9585
2	7 November	21:45:21	496,058	891,415	366	10,067
3	7 November	22:00:10	495,912	891,456	354	57,543
…	…	…	…	…	…	…

, each record contains the time, location, and energy information of microseism.

These microseism events are presented in the mining engineering plan, as shown in Figure 6. It can be seen that microseism events cover the whole target working face.

The distribution of microseisms in three-dimensional space is shown in Figure 7.

4.2. Clustering of Microseism Events

The ST-DBSCAN algorithm is used to cluster microseism events in the target working face. According to the empirical formula (3), $Minpts=lnN\approx 9$ can be calculated. The spatial radius $Eps$ and time window $\mathsf{\Delta }T$ of the cluster are obtained from the K-Dist graph [32], where K takes the value of $Minpts$, which is 9 at this time. The specific steps are to calculate the spatial distance and temporal distance between each sample point and the ninth closest point, and draw the corresponding 9-Dist graph in descending order, the distance value corresponding to the inflection point is the effective radius of dividing the noise point and the non-noise point. The 9-Dist graph is shown in the following Figure 8.

By observing Figure 8, the vertical axis coordinate values corresponding to the two inflection points are 80 and 20, respectively. However, considering the influence of human factors, there may be subjective errors in the vertical axis coordinate values. In order to find the parameters with the best clustering results, different distance values corresponding to inflection points were selected to carry out the parameter combination clustering tests. $Eps$ was taken as [70 80 90], and $\mathsf{\Delta }T$ was taken as [15 20 25]. In

Table 2. Parameter combinations of $Eps$ and $\mathsf{\Delta }T$.

$$\mathit{E}\mathit{p}\mathit{s}\left(\mathbf{m}\right)$$	$$\mathbf{\Delta }\mathit{T}\left(\mathbf{h}\right)$$
	15	20	25
70	ID1	ID4	ID7
80	ID2	ID5	ID8
90	ID3	ID6	ID9

, nine parameter combinations can be obtained based on the different values of $Eps$ and $\mathsf{\Delta }T$, numbered ID1 to ID9. respectively.

In order to measure the effect of clustering, the Silhouette Coefficient (SC) is selected as the clustering effectiveness assessment index [33]. Generally speaking, the higher the average SC of the samples in the clusters, the better the clustering quality. For any point i in the cluster, the SC is calculated as follows:

$$SC\left(i\right)=\frac{b\left(i\right)-a\left(i\right)}{\max \left\{a\left(i\right),b\left(i\right)\right\}}$$

(18)

where $a\left(i\right)$ is the degree of dissimilarity within the cluster, representing the average value of the degree of dissimilarity from point $i$ to other points in the same cluster, reflecting the degree of cohesion; $b\left(i\right)$ is the degree of dissimilarity between clusters, representing the minimum value of the average degree of dissimilarity from point $i$ to other clusters, reflecting the degree of separation. It can be seen that the value of SC is between [−1, 1], and the closer it is to 1, the better cohesion and separation degree are. The average of the SC of all points is the comprehensive SC of the clustering result.

The clustering algorithm program is run on the nine parameter combinations in Table 2, and the output results are the specific cluster sets and noise points that are not classified into any clusters. Different parameter combinations of $Eps$ and $\mathsf{\Delta }T$ will get different clustering results. Their SC are calculated respectively and the noise rate is counted. The results are shown in

Table 3. SC and noise rate of parameter combinations.

Combinations	ID1	ID2	ID3	ID4	ID5	ID6	ID7	ID8	ID9
SC	0.284	0.212	0.271	0.296	0.32	0.229	0.274	0.213	0.246
Noise Rate	29.50%	23.10%	18.50%	24.20%	12.20%	14.70%	20.50%	15.50%	19.00%

and Figure 9.

It can be easily seen from Figure 9 that the SC of ID5 is high and the noise rate is low. Therefore, when $Eps$ is set to 80 and $\mathsf{\Delta }T$ is set to 20, the clustering result is the best. Different colors are used to represent different clusters, and the clustering result of this parameter combination is presented in three-dimensional space, as shown in Figure 10.

At this time, the noise rate is 12.2%, and 47 microseism clusters are obtained, and the largest cluster contains 1077 events. Due to space limitation, some statistical results are shown in

Table 4. Clustering results.

Cluster Number	Amount	Cluster Number	Amount
1	60	11	1077
2	527	12	95
3	128	13	71
4	749	14	948
5	149	15	151
6	807	16	914
7	74	17	193
8	135	18	286
9	98	19	59
10	104	20	52
…	…	…	…

.

Several clusters containing more microseism events are projected into the mining engineering plan. The results are shown in Figure 11. Compared with Figure 6, it can be seen that microseism events have obvious aggregation. Clustering analysis of microseism events can be regarded as a preprocessing process for subsequent prediction data. The purpose is to remove interference factors, enhance the correlation between sample data and improve the accuracy of prediction models.

4.3. Prediction of Mine Earthquake

The microseism clusters can be regarded as mine earthquake sequences. The location and frequency corresponding to different energy levels of events in the cluster have a time-series law, that is, the location and frequency of subsequent microseism events are related to historical microseism information. The largest microseism cluster containing 1077 events was selected, and the SVR was used to establish a model to conduct prediction experiment in sequence for the microseism events in the cluster. The prediction includes two aspects: occurrence location, and daily frequency of high-energy event. Generally, when the energy released by mine earthquake reaches above ${10}^4 \mathrm{J}$, the rock burst may be induced, causing serious injury to underground roadway and personnel [34]. Therefore, the high-energy event refers to the mine earthquake event that releases energy greater than ${10}^4 \mathrm{J}$. In data preprocessing, considering the strong correlation of microseism events within the cluster and the limited number of events, the sampling step of position prediction is 10 times, the sampling step of frequency prediction is 3 days, and the moving step is 1. In order to evaluate the prediction effect of model, MSE and R² are selected as assessment indexes [35], and their calculation formulas are as follows:

$$MSE=\frac{1}{n}\underset{i=1}{\stackrel{n}{{{\displaystyle \sum }}^{ }}}{\left(y_i-{\widehat{y}}_i\right)}^2$$

(19)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(20)

where $n$ is the total number of samples, ${\widehat{y}}_i$ is the predicted value, $y_i$ is the true value of the sample, $\overline{y}$ is the mean of the true value. The smaller the value of MSE is, the more accurate the model prediction is. R² is between (−∞, 1], and the larger the value is, the better the effect of model fitting ability is.

4.3.1. Prediction for Location of High-Energy Events

It is of great significance to predict the location of high-energy mine earthquake for the prevention and control of underground dynamic disasters in coal mines. The specific method is extracting the high-energy events with release energy greater than ${10}^4 \mathrm{J}$ from the cluster, using the related information of the previous 10 mine earthquakes to predict the location of the next event. The precursory pattern sequence features of the model include the X, Y, Z coordinate values and the distance between two events immediately adjacent to each other of the previous 10 events. The corresponding label is the location of the next event, that is, the X, Y, Z coordinate values. As shown in Figure 12 below, $P=\left(p_1,p_2,p_3,\dots ,p_n\right)$ is the time series set of microseism location information, where $p_i$ is the location information of one event. When the sampling step is 10 and the moving step is 1, the precursory pattern sequence set $S=\left(s_1,s_2,s_3,\dots ,s_u\right)$ can be obtained, where $s_i$ is a precursory pattern sequence, which contains the location information of the previous ten microseism events. The label set corresponding to $S$ is $L=\left(l_1,l_2,l_3,\dots ,l_u\right)$,where $l_i$ is the label corresponding to $s_i$ and contains the location coordinates of the next microseism event.

The statistical result was $n=$ 97. After $P$ was processed as shown in the Figure 12, $u=$ 87 was obtained, that is, the prediction model would have 87 samples and corresponding labels. The samples were divided into training set and test set according to the ratio of 7:3. The training set was used for the simulation training of the model, and the test set was used to test the effect of the training. Four kernel functions of SVR, Linear, Sigmoid, Polynomial and RBF, were selected for the experiment, and the grid search method was used to optimize the hyperparameters of model. Overfitting can be suppressed by two parameters. One is the penalty factor C, the higher its value is, the more intolerant the model is to errors and easy to overfit. The other is the parameter gamma in the kernel function. With the increase of gamma, the prediction effect of the training set becomes better, and the prediction effect of the test set becomes worse. At the same time, the complexity of the model increases, and the generalization ability becomes worse, which leads to the overfitting. The overfitting issue is controlled by adjusting the above two parameters combined with the final prediction effect of the model.

The predicted values outputted from model and the true values were compared together to get Figure 13 below. In view of the equal position of X, Y and Z coordinates in location prediction, this case only presents the prediction of X coordinates due to space limitation. The MSE and R² of the four prediction models with different kernel functions in the figure were calculated respectively, and the statistical results are shown in

Table 5. Model performance with different kernel functions.

Assessment Indexes	Kernel Functions
	Linear	Sigmoid	Polynomial	RBF
MSE	0.0355	0.0403	0.0081	0.0076
R2	0.3642	0.2772	0.8517	0.8636

.

Combined with Figure 13 and Table 5, it can be seen that the SVR model with RBF kernel function has the best effect on the location prediction of high-energy mine earthquakes.

4.3.2. Prediction for Frequency of High-Energy Events

For the microseism events in the cluster, the occurrence frequencies corresponding to different energy levels are counted every day. The energy levels are divided according to ${10}^2 \mathrm{J}$, ${10}^3 \mathrm{J}$, ${10}^4 \mathrm{J}$, ${10}^5 \mathrm{J}$... As with the location prediction, only high-energy events with energy above ${10}^{4 }\mathrm{J}$ are selected as the prediction objects. The specific method is using the related frequency information of the previous three days, to predict the occurrence frequency of high-energy events in the next day. Therefore, the precursory pattern sequence features of the model are the daily occurrence frequency of microseism events with different energy levels in the previous three days, and the corresponding label is the occurrence frequency of high-energy mine earthquakes with different energy level in the next day. As shown in Figure 14, $D=\left(d_1,d_2,d_3,\dots ,d_m\right)$ is the time series set of microseism energy information, where $d_i$ is the energy frequency data within a day. When the sampling step is 3 and the moving step is 1, the precursory pattern sequence set $S=\left(s_1,s_2,s_3,\dots ,s_v\right)$ can be obtained, where $s_i$ is a precursory pattern sequence, which contains the energy frequency information of the previous three days. The corresponding label set is $L=\left(l_1,l_2,l_3,\dots ,l_v\right)$, where $l_i$ is the label corresponding to $s_i$ and contains the occurrence frequency of mine earthquakes with energy level above ${10}^4\mathrm{J}$ in the following day.

The statistical result was m = 31. After $D$ was processed as shown in the Figure 14, v = 28 was obtained, that is, the prediction model would have 28 samples and corresponding labels. The samples were divided into training set and test set according to the ratio of 7:3. The training set was used for the simulation training of the model, and the test set was used to test the effect of the training. Four kernel functions of SVR, Linear, Sigmoid, Polynomial and RBF, were selected for the experiment, and the grid search method was used to optimize the hyperparameters of model. The predicted values outputted from model and the true values were compared together to get Figure 15. In view of the equal position of different energy levels in energy frequency prediction, this case only presents the prediction of daily frequency of ${10}^4 \mathrm{J}$ level mine earthquake due to space limitation. The MSE and R² of the four prediction models with different kernel functions in the figure were calculated respectively, and the statistical results are shown in

Table 6. Model performance with different kernel functions.

Assessment Indexes	Kernel Functions
	Linear	Sigmoid	Polynomial	RBF
MSE	1.8984	1.8063	0.6344	0.2251
R2	−0.3962	−0.3301	0.5328	0.8342

.

Combined with Figure 15 and Table 6, it can be seen that the SVR model with RBF kernel function has the best prediction effect on the daily frequency of high-energy mine earthquakes with different energy levels.

4.4. Comparison of Prediction Models

As a comparison experiment, the LSTM was used to build a prediction model, and the model was trained on the microseism dataset (including 8497 event records) of 103# working face without clustering process, and used to predict the location of high-energy mine earthquakes and the daily occurrence frequency corresponding to different energy levels. In the SVR model above, the microseism data used was distributed from April 23 to May 23, so the event samples in the test set of the LSTM prediction model were also set within this time period. The time series sets of microseism location information and energy frequency information were also processed respectively as shown in Figure 12 and Figure 14 to obtain $n=$ 1320, $u=$ 1310, $m=$ 562, $v=$ 559. The Early Stopping method was used to prevent the overfitting issue. The specific method is to calculate the accuracy of validation data at the end of each epoch training. When the accuracy does not improve any more, the model training will be stopped to avoid overfitting. Similarly, the grid search method was used to optimize the hyperparameters of model, and the prediction of X coordinates of high-energy mine earthquakes and daily occurrence frequency corresponding to ${10}^4 \mathrm{J}$ energy level events are shown in Figure 16:

In the SVR model, the RBF kernel function has the best prediction effect. Therefore, the MSE and R² of the LSTM model in predicting the X coordinates of high-energy events and the daily frequency corresponding to ${10}^4 \mathrm{J}$ energy level events are calculated and compared with the former. The results are shown in

Table 7. Comparison of model performance.

Assessment Indexes	X Coordinate		$\mathbf{Daily} \mathbf{Frequency} \mathbf{of} {10}^4 \mathbf{J} \mathbf{Level} \mathbf{Events}$
	LSTM	SVR_RBF	LSTM	SVR_RBF
MSE	0.0314	0.0076	0.5387	0.2251
R2	0.4602	0.8636	0.8098	0.8342

.

Through comparative experiments, it can be concluded that the SVR model with RBF kernel function on the basis of clustering analysis is better than the LSTM model without clustering analysis, both for the prediction of the location of high-energy mine earthquakes and the prediction of daily frequency of high-energy mine earthquakes with different energy levels.

4.5. Out-of-Sample Prediction Experiment

Save the trained model and feed it the microseism data of a new cluster, which was not involved in the previous model training process. The final performance is shown in

Table 8. Model performance for out-of-sample prediction.

Assessment Indexes	X Coordinate		$\mathbf{Daily} \mathbf{Frequency} \mathbf{of} {10}^4 \mathbf{J} \mathbf{Level} \mathbf{Events}$
	LSTM	SVR_RBF	LSTM	SVR_RBF
MSE	0.0298	0.0471	0.6127	1.2062
R2	0.4954	0.4383	0.7351	0.5281

. It can be seen that the prediction effect of the trained SVR model with RBF kernel on the new cluster is not as good as that of the LSTM model using global data. This indicates that the time series laws of data in different clusters may not be the same, due to the influence of the changing geological conditions and other factors, and the trained SVR model cannot be simply and directly applied to the new cluster. Its parameters need to be adjusted to get a good performance on different clusters.

5. Conclusions

It is easy for mine earthquake activity to cause coal mine dynamic disasters, resulting in serious consequences. The microseism monitoring system can monitor the microseism events in real-time and record the time, space, and intensity information when the mine earthquake occurs. In this paper, the microseism data of 103# working face of the H Coal Mine were selected. Firstly, the ST-DBSCAN algorithm was used to find out the microseism clusters, and then SVR was used to predict the occurrence location and daily frequency of high-energy events on this basis. MSE and R² were selected as the assessment indexes of the model prediction effect. The results (Table 5 and Table 6) showed that among the four different kernel functions, the SVR prediction model with RBF kernel function has the best performance. Finally, as a comparison experiment, the LSTM model was used to predict high-energy mine earthquakes in the same period using the 103# working face microseism data without clustering process. Through comparative analysis, its performance (Table 7) was inferior to the SVR prediction model with RBF kernel function. This engineering experiment showed that the spatial-temporal clustering analysis of microseism events in advance can improve the prediction accuracy of the machine learning method on high-energy mine earthquakes.

However, results of out-of-sample prediction experiment (Table 8) showed that the trained SVR model does not perform well on the new cluster. The reason could be that the time series laws of data in different clusters are not exactly the same, due to the influence of the changing geological conditions and other factors. In future work, we will conduct research on the model’s generalization ability on different microseism events clusters and mines.