Microseismic Data-Driven Short-Term Rockburst Evaluation in Underground Engineering with Strategic Data Augmentation and Extremely Randomized Forest

Abstract

Rockburst is a common dynamic geological disaster in underground mining and tunneling engineering, characterized by randomness, abruptness, and impact. Short-term evaluation of rockburst potential plays an outsize role in ensuring the safety of workers, equipment, and projects. As is well known, microseismic monitoring serves as a reliable short-term early-warning technique for rockburst. However, the large amount of microseismic data brings many challenges to traditional manual analysis, such as the timeliness of data processing and the accuracy of rockburst prediction. To this end, this study integrates artificial intelligence with microseismic monitoring. On the basis of a comprehensive consideration of class imbalance and multicollinearity, an innovative modeling framework that combines local outlier factor-guided synthetic minority oversampling and an extremely randomized forest with C5.0 decision trees is proposed for the short-term evaluation of rockburst potential. To determine the optimal hyperparameters, the whale optimization algorithm is embedded. To prove the efficacy of the model, a total of 93 rockburst cases are collected from various engineering projects. The results show that the proposed approach achieves an accuracy of 90.91% and a macro F₁-score of 0.9141. Additionally, the local F₁-scores on low-intensity and high-intensity rockburst are 0.9600 and 0.9474, respectively. Finally, the advantages of the proposed approach are further validated through an extended comparative analysis. The insights derived from this research provide a reference for microseismic data-based short-term rockburst prediction when faced with class imbalance and multicollinearity.

1. Introduction

Rockburst is a typical dynamic geological disaster, frequently occurring in underground mining and tunneling engineering [1,2,3,4]. Rockburst is caused by the sharp release of elastic energy accumulated in rock mass due to the unloading effect and external disturbances, often accompanied by the violent ejection of rock fragments [5]. Rockburst brings serious threats to the safety of operators and equipment and also results in the deformation and even destruction of supporting structures and excavation systems. The world’s first recorded rockburst occurred in 1738 at a tin mine in England [6], and, after that, rockburst has been reported in many countries. Between 1936 and 1993, the United States experienced 172 rockburst incidents, followed by 78 fatalities [7]. From the early 1980s to the mid2000s, Germany went through 42 recorded rockburst incidents with death and injuries [8]. In addition, there were also many rockburst incidents in Australia, Canada, South Africa, and other mining countries [9,10,11]. The occurrence of rockburst is closely related to the characteristics of minerals and rocks. In general, the higher the content of siliceous or other hard minerals, the more prone to rockburst. In terms of rock types, highly brittle rocks such as quartzite, granite, marble, and sandstone are prone to rockburst. Due to the randomness in space and the suddenness in time, rockburst prediction is faced with great challenges.

The microseismic monitoring technique is developed based on acoustic emissions and seismology [12,13,14,15,16]. When underground rocks fracture due to human or natural factors, seismic waves will be generated and spread around. By arranging multiple detectors in the three-dimensional space around the fracturing zone, the microseismic monitoring technique can realize the real-time in situ observation of the fracturing process. Starting from the initial deformation stage, the microseismic monitoring technique can capture the full failure process of rock mass, including fracture initiation, fracture propagation, and structure instability. This makes it so that the microseismic monitoring technique can detect the precursor to rockburst and provide an important basis for its early warning [17,18,19,20].

However, as a real-time monitoring technique, massive microseismic data bring great challenges to the timeliness and accuracy of rockburst early warning. To this end, many researchers have combined artificial intelligence with microseismic monitoring, that is, using machine learning to mine microseismic data and establish the mapping relationship with rockburst. Jin et al. [21] extracted six parameters from microseismic data, including cumulative number, cumulative energy, cumulative apparent volume, changing rate of cumulative number, changing rate of cumulative energy, and changing rate of cumulative apparent volume, and used a support vector machine to build the short-term rockburst prediction model. Similarly, based on these six microseismic parameters, Liang et al. [22] employed a gradient boosting decision tree and a random forest to build the short-term rockburst prediction model, and Feng et al. [23] adopted a probabilistic neural network to construct the short-term rockburst prediction model. Subsequently, Zhou et al. [24] and Liang et al. [25] added a new factor, namely incubation time, to the above parameter system and developed the ensemble learning-based short-term rockburst prediction model. In practical applications, these models have made indelible contributions. Class unbalance and multicollinearity are two common data defects, and their effects on model performance cannot be ignored [26,27,28,29,30,31]. Consequently, this study, on the basis of the comprehensive consideration of class imbalance and multicollinearity, proposes a novel hybrid intelligent modeling approach for short-term assessment of rockburst potential.

In this approach, the local outlier factor-guided synthetic minority oversampling technique and extremely randomized forest with C5.0 decision trees are proposed from the data and algorithm aspects, respectively. Additionally, the whale optimization algorithm is embedded to synergistically optimize the critical hyperparameters of the data strategy and algorithm strategy. The rest of the paper is organized as follows: Section 2 describes the proposed modeling framework in detail; Section 3 depicts the collected database; Section 4 gives the modeling results and conducts some profound discussion; Section 5 summarizes the main conclusions.

2. Proposed Modeling Framework

Class imbalance and multicollinearity are two common types of data defects in geotechnical and geological engineering. In order to remove their adverse effects on the prediction performance of machine learning models, improvements are made at the data level and algorithm level. At the data level, we propose local outlier factor-guided synthetic minority oversampling (LOF-SMO). At the algorithm level, we propose an extremely randomized forest with C5.0 decision trees (C5.0DT-ERF). To determine the optimal hyperparameters in the LOF-SMO and C5.0DT-ERF, the whale optimization algorithm (WOA) is merged. As a result, the WOA-LOF-SMO-C5.0DT-ERF modeling framework is devised (Figure 1).

2.1. Local Outlier Factor-Guided Synthetic Minority Oversampling

2.1.1. Local Outlier Factor

Local outlier factor (LOF) is a critical indicator to distinguish outliers and non-outliers in the database [32]. An object is determined to be an outlier if its LOF exceeds the threshold of 1.5. In order to capture the dominant topology of the data in subsequent oversampling, synthetic minority oversampling is implemented in the database without outliers. The calculation principle of the LOF is described as follows:

(1)

Determine the k-distance neighborhood $N_k\left(O\right)$ of the object O, consisting of the k nearest neighbors of O in the database D.

(2)

Calculate the k-distance $dis{t}_k\left(O\right)$ of O by:

$$dis{t}_k\left(O\right)=\max \left\{dist\left(O,P\right)\left|P\in N_k\left(O\right)\right.\right\}$$

(1)

where P is a neighbor from $N_k\left(O\right)$; $dist\left(O,P\right)$ is the Euclidean distance between O and P.

(3)

Calculate the reachable distance $dis{t}_{reach}\left(O,P\right)$ between O and P by:

$$dis{t}_{reach}\left(O,P\right)=\max \left\{dis{t}_k\left(O\right),dist\left(O,P\right)\right\}$$

(2)

(4)

Calculate the local reachable density $lrd\left(O\right)$ of O by:

$$lrd\left(O\right)={\displaystyle \frac{k}{{\displaystyle {\sum }_{P\in N_k\left(O\right)}dis{t}_{reach}\left(O,P\right)}}}$$

(3)

(5)

Calculate the LOF $lof\left(O\right)$ of O by:

$$lof\left(O\right)={\displaystyle \frac{{\displaystyle {\sum }_{P\in N_k\left(O\right)}\frac{lrd\left(P\right)}{lrd\left(O\right)}}}{k}}$$

(4)

where $lrd\left(P\right)$ is the local reachable density of P.

2.1.2. Synthetic Minority Oversampling

Synthetic minority oversampling (SMO) is an oversampling method used to deal with class imbalance and increased sample size, especially when the number of minority samples is small [33]. Its basic principle is to generate new minority samples by linear interpolation, as shown in Figure 2. The detailed flow of the SMO is described as follows:

(1)

Select minority samples: each minority sample in turn is selected as the root sample for synthesizing new samples.

(2)

Find nearest neighbors: for each root sample, using Euclidean distance as the standard, its distance to all other minority samples is calculated to obtain its K nearest neighbors.

(3)

Select auxiliary samples: one sample is randomly selected from the K nearest neighbors of each root sample as the auxiliary sample for synthesizing new samples.

(4)

Synthesize new samples: between the root sample and the auxiliary sample, a new sample is generated by linear interpolation. The interpolation formula is expressed as:

$$x_i^{ne}=x_i^{ro}+{\lambda }_i\cdot \left(x_i^{au}-x_i^{ro}\right)$$

(5)

where $x_i^{ne}$, $x_i^{ro}$ and $x_i^{au}$ are the i-th feature of the new sample, the root sample, and the auxiliary sample, respectively; ${\lambda }_i$ is a random number in [0, 1].

(5)

Repeat generation: for each root sample, the above steps are repeated until the number of new samples meets the requirements.

(6)

Add new samples: all generated new samples are added to the original dataset, thus increasing the number of minority samples and making the dataset more balanced.

2.2. Extremely Randomized Forest with C5.0 Decision Trees

2.2.1. C5.0 Decision Tree

A decision tree is a group of procedures designed to classify the input data into more homogenous subsets using generated rules [34]. In the training process, the decision tree aims to maximize the information gain and minimize the entropy in the generated subsets. Figure 3 demonstrates the basic components of a decision tree. Initially, all data are gathered in a root node and then split into relatively more homogenous subsets by internal nodes based on feature thresholds. The splitting process continues until the decision tree reaches the leaf nodes. At this final stage, labels are assigned to the leaf nodes.

A C5.0 decision tree (C5.0DT) [35] is a representative decision tree algorithm that uses the information gain ratio to select split attributes and thresholds in nodes and finish tree growth [36]. In this process, the decrease in the entropy is maximized, and consequently the purity of leaf nodes reaches the highest. Referring to Equation (6), the calculation principle of the information gain ratio is presented. In order to reduce the error when processing new data after tree growth, a decision tree must be pruned. In the C5.0DT, pruning automatically starts in reverse from the leaf node and extends upwards to the entire tree based on the information gain ratio [37].

$$Gain\_ratio\left(S,A\right)={\displaystyle \frac{Gain\left(S,A\right)}{-{\displaystyle {\sum }_{i=1}^v\frac{\left|S_i\right|}{\left|S\right|}{\log }_2\frac{\left|S_i\right|}{\left|S\right|}}}}$$

(6)

where $Gain\_ratio\left(S,A\right)$ is the information gain ratio obtained by splitting the node S using the attribute A; $Gain\left(S,A\right)$ is the information gain obtained by splitting the node S using the attribute A, calculated as [38]; v is the number of child nodes generated by splitting the node S; $S_i$ is the i-th child node generated by splitting the node S; ∣ ∣ is the number of samples contained in the node.

2.2.2. Extremely Randomized Forest

An extremely randomized forest (ERF) is made up of a series of standalone trees, as shown in Figure 4 [39]. Each tree generates the prediction results independently according to the input data, and then on this basis, the prediction results of the ERF are determined by voting on all trees, that is, majority rules. In this study, a C5.0 decision tree is used as the basic tree unit of the ERF for the subsequent analysis. The modeling process of the ERF is described as follows:

(1)

Conduct bootstrap sampling: bootstrap sampling is performed on the balanced dataset of Section 2.1 so as to generate L subdatasets, where the number of subdatasets is the same as the number of tree models in the ERF.

(2)

Build tree models: L tree models are built individually based on L subdatasets. To determine the appropriate split attribute and split value during tree growth, the information gain ratio is used as the metric.

(3)

Increase the diversity: To increase the diversity of the tree models, the extremely randomized strategy is implemented during modeling. Specifically, for each tree model, the optimal split attribute is produced by winning from the attribute subset randomly selected in all candidate attributes, and the split value is from the value subset randomly generated in the candidate range.

(4)

Integrate tree models: to determine the output of the ERF, voting is carried out on the basis of the prediction results of all tree models.

2.3. Hyperparameter Optimization

The whale optimization algorithm (WOA) is a metaheuristic optimization algorithm that simulates the hunting behavior of humpback whales [40]. In the exploration phase, the mathematical model is described as:

$$X_i\left(t+1\right)=X_{\mathrm{rand}}\left(t\right)-A\left|C\cdot X_{\mathrm{rand}}\left(t\right)-X_i\left(t\right)\right|$$

(7)

where $X_i\left(t\right)$ and $X_i\left(t+1\right)$ are the position of the i-th whale at the t-th and (t + 1)-th iterations, respectively; $X_{\mathrm{rand}}\left(t\right)$ is the position of a whale randomly selected from the population; and A and C are changeable coefficients, defined as:

$$A=2a\cdot r-a$$

(8)

$$C=2r$$

(9)

where a decreases linearly from 2 to 0 with iteration; r is a random number in [0, 1].

In the exploitation phase, the mathematical model is described as:

$$X_i\left(t+1\right)=\begin{array}{cc}{X}_{\mathrm{best}}\left(t\right)-A\left|C\cdot X_{\mathrm{best}}\left(t\right)-X\left(t\right)\right|& \mathrm{when}\hspace{0.33em}p<0.5\end{array}$$

(10)

$$X_i\left(t+1\right)=\begin{array}{cc}{e}^{bl}\left|X_{\mathrm{best}}\left(t\right)-X\left(t\right)\right|\cos \left(2\pi l\right)+X_{\mathrm{best}}\left(t\right)& \mathrm{when}\hspace{0.33em}p\ge 0.5\end{array}$$

(11)

where $X_{\mathrm{best}}\left(t\right)$ is the optimal position searched by the population; b is the constant that defines the logarithmic spiral shape; l is the random number in [−1, 1]; p is a random number in [0, 1].

In order to balance the global and local search, the exploration and exploitation phases alternate, as demonstrated in Figure 5. In the proposed short-term evaluation framework for rockburst potential, there are five critical hyperparameters that need to be optimized, including the size of the k-distance neighborhood in the LOF (referred to as x₁), the number of used nearest neighbors in the SMO (referred to as x₂), the minimum number of split samples in the C5.0DT (referred to as x₃), the maximum number of leaf nodes in the C5.0DT (referred to as x₄), and the number of C5.0 decision trees in the ERF (referred to as x₅). Therefore, the position of the whale is a five-dimensional vector, denoted as $\mathit{X}=\left[x_1,\hspace{0.33em}{x}_2,x_3,x_4,x_5\right]$.

3. Database Description

3.1. Case Collection

To demonstrate the validity of the proposed approach, 93 rockburst cases in total are collected from Feng et al. [41] (see Supplementary Material). According to

Table 1. Classification criteria of rockburst intensity [22,41].

Rockburst Intensity	Failure Characteristics
None	No rockburst occurs. No abnormality in surrounding rock. Normal construction.
Slight	The depth of rockburst crater is <0.5 m. Slight spalling or slabbing. The size of ejected rock fragment is 10~30 cm. Slight cracking sound.
Moderate	The depth of rockburst crater is 0.5~1.0 m. Severe spalling and slabbing. The size of ejected rock fragment is 30~80cm. Detonator blasting-like sound.
Strong	The depth of rockburst crater is >1.0 m. Extensive spalling and slabbing. The size of ejected rock fragment is >80 cm. Explosion-like sound with an impact wave.

, rockburst intensity is split into four grades: none, slight, moderate, and strong. Among these 93 rockburst cases, 34 cases (accounting for 36.56%) belong to none rockburst, 21 cases (accounting for 22.58%) belong to slight rockburst, 25 cases (accounting for 26.88%) belong to moderate rockburst, and 13 cases (accounting for 13.98%) belong to strong rockburst. As observed in Figure 6, the class imbalance is prominent.

In order to predict rockburst intensity, six microseismic source parameters are used as the input parameters of the model, including cumulative number (unit), cumulative energy (J), cumulative apparent volume (m³), changing rate of cumulative number (unit/day), changing rate of cumulative energy (J/day), and changing rate of cumulative apparent volume (m³/day). During the rock fracturing process, the microseismic monitoring technique captures the stress wave emitted by the source in real time and records it as the microseismic event [42,43]. The cumulative number of microseismic events indicates the microseismic activity, the cumulative energy represents the microseismic strength, and the cumulative apparent volume reflects the damage degree of rock mass [24]. On the whole, these three microseismic source parameters jointly characterize the frequency, strength, and scale of rock fracture and are able to comprehensively describe the fracture situation inside the rock mass. Considering that rockburst incubation is a dynamic process, three microseismic source parameters with the time effect are also taken into account in the input system, namely the changing rate of cumulative number, changing rate of cumulative energy, and changing rate of cumulative apparent volume [21].

Specifically, the microseismic energy and apparent volume are calculated by Equations (12) and (13), respectively [21,44,45]. The time interval for the calculation of the cumulative value is the incubation time of rockburst [24]. Assuming that the incubation time is $N$ days, the cumulative value $C_v$ is the sum of the value $c_n\hspace{0.33em}\left(n=1,2,\cdots ,N\right)$ of each day in these $N$ days, calculated as Equation (14). Correspondingly, the changing rate $C_r$ is the ratio of the cumulative value $C_v$ to the incubation time $N$, calculated as Equation (15) [21]. Taking the strong rockburst that occurred at milestone SK8+709 of the Jinping Ⅱ hydropower station on 11 January 2011 as an example, Figure 7 shows the evolution of the cumulative number of microseismic events from 6–10 January, where $C_v=49$ and $C_r=12.25$ according to Equations (14) and (15), respectively [41].

$$E=8\pi \rho v{\displaystyle {\int }_0^{\infty }{s}^2\left(f\right)df}$$

(12)

$$V_A=\mu P^2/E$$

(13)

where $E$ is the microseismic energy; $V_A$ is the apparent volume; $\rho$ is the density of rock mass; $v$ is the wave velocity; $s^2\left(f\right)$ is the velocity power spectrum; $f$ is the frequency; $\mu$ is the shear modulus of rock mass; $P$ is the seismic potency.

$$C_v=c_1+c_2+\cdots +c_N$$

(14)

$$C_r=C_v/N$$

(15)

where $C_v$ is the cumulative value; $C_r$ is the changing rate; $c_n\hspace{0.33em}\left(n=1,2,\cdots ,N\right)$ is the individual value of each day during rockburst incubation; $N$ is the incubation time of rockburst.

Table 2. Statistical description of input parameters.

Rockburst Intensity	Input Parameter	Statistical Index
		Number	Minimum	Maximum	Mean	Median	Skewness	Kurtosis
None	CN	34	1	17	3.94	3	2.04	5.22
	CE	34	0.78	5.82	3.16	3.55	−0.14	−1.30
	CAV	34	2.51	4.86	3.62	3.58	0.27	−0.87
	CNR	34	0.11	2.50	0.85	0.76	1.21	0.94
	CER	34	0.18	4.78	2.51	2.81	−0.14	−1.20
	CAVR	34	1.67	4.31	2.97	2.96	−0.19	−0.40
Slight	CN	21	3	29	10.14	8	1.48	2.02
	CE	21	3.54	5.56	4.54	4.53	0.05	−0.13
	CAV	21	3.50	4.94	4.18	4.13	0.12	−0.57
	CNR	21	0.54	4.00	1.48	1.11	1.29	1.23
	CER	21	2.84	4.80	3.71	3.67	0.46	−0.36
	CAVR	21	2.39	3.99	3.35	3.50	−0.64	−0.54
Moderate	CN	25	3	36	15.12	14	0.80	1.61
	CE	25	3.54	5.98	5.13	5.10	−0.93	0.96
	CAV	25	3.52	4.87	4.48	4.57	−1.41	3.32
	CNR	25	0.43	4.00	1.70	1.71	0.88	1.92
	CER	25	2.29	5.08	4.12	4.25	−1.01	0.60
	CAVR	25	2.67	4.02	3.51	3.55	−0.58	0.61
Strong	CN	13	10	70	37.31	42	−0.09	−0.80
	CE	13	4.11	7.09	5.94	6.15	−1.02	1.02
	CAV	13	3.62	5.17	4.87	4.98	−2.76	8.72
	CNR	13	1.25	12.25	4.53	3.73	1.48	2.88
	CER	13	3.41	5.89	5.01	5.15	−0.92	−0.03
	CAVR	13	2.93	4.39	3.94	4.08	−1.52	3.07

Note: The CN, CE, CAV, CNR, CER, and CAVR denote cumulative number, cumulative energy, cumulative apparent volume, changing rate of cumulative number, changing rate of cumulative energy, and changing rate of cumulative apparent volume, respectively. Particularly, the values of the CE, CAV, CER, and CAVR are expressed in logarithmic form with base 10.

shows the statistical description of input parameters for various rockburst intensities, and the visual distribution of input parameters is presented in Figure 8.

3.2. Correlation Analysis

Generally, input parameters should be independent of each other, and the stronger the parameter correlation, the higher the information redundancy. The Pearson correlation coefficient is adopted to analyze the correlation of the six input parameters selected in Section 3.1, calculated by Equation (16) [46]. Based on the magnitude of the Pearson correlation coefficient, the correlation strength is divided into five levels, among which 0.0~0.2, 0.2~0.4, 0.4~0.6, 0.6~0.8, and 0.8~1.0 indicate extremely weak, weak, moderate, strong, and extremely strong correlation, respectively. In particular, a positive Pearson correlation coefficient represents positive correlation, and a negative Pearson correlation coefficient stands for negative correlation.

$${\rho }_{X,Y}={\displaystyle \frac{{\displaystyle \sum \left(X-\overline{X}\right)\left(Y-\overline{Y}\right)}}{\sqrt{{\displaystyle \sum {\left(X-\overline{X}\right)}^2{\displaystyle \sum {\left(Y-\overline{Y}\right)}^2}}}}}$$

(16)

where ${\rho }_{X,Y}$ is the Pearson correlation coefficient between input parameters X and Y; $\overline{X}$ and $\overline{Y}$ are the mean of input parameters X and Y, respectively.

Referring to Figure 9, there exists a strong correlation between cumulative number and its changing rate, with the Pearson correlation coefficient reaching 0.7866. Similarly, a strong correlation can also be found between cumulative energy and its changing rate, as well as cumulative apparent volume and its changing rate, with the Pearson correlation coefficient reaching 0.9747 and 0.8891, respectively. Regarding other parameter pairs, no obvious correlation appears. In essence, cumulative number, cumulative energy, and cumulative apparent volume are static indicators, while their changing rate, as a dynamic indicator, reflects the evolution trend of static indicators during rockburst development. Therefore, static and dynamic indicators have an inevitable physical correlation. Since both static and dynamic metrics have clear physical meaning, no additional measures, such as dimensionality reduction, are taken to eliminate correlations herein. Moreover, the dimension of the current input combination is 6, which does not cause a dimensional disaster in modeling.

3.3. Multicollinearity Analysis

When there is multicollinearity between input parameters, small changes in the data can lead to significant changes in the parameter estimates, resulting in a decrease in the stability of the model. Variance inflation factor (VIF) is a common measure of multicollinearity severity, and VIF greater than 10 is generally considered to indicate strong multicollinearity [47]. VIF is calculated by:

$$VI{F}_i={\displaystyle \frac{1}{1-{R_i}^2}}$$

(17)

where $R_i$ is the determination coefficient of the linear fit of the input parameter $X_i$ with respect to other input parameters.

Referring to Figure 10, cumulative energy, cumulative apparent volume, changing rate of cumulative energy, and changing rate of cumulative apparent volume behave with strong multicollinearity, with VIF reaching 153, 35, 130, and 29, respectively. And, cumulative number and changing rate of cumulative number behave with weak multicollinearity, with VIF reaching 6 and 5, respectively.

4. Results and Discussion

4.1. Evaluation Metrics

In order to quantitatively evaluate the generalization performance of the proposed model, two evaluation metrics are selected, including the accuracy and macro F₁-score. The accuracy is calculated by:

$$\mathrm{accuracy}={\displaystyle \frac{{\displaystyle {\sum }_{k=1}^K{u}_{k,k}}}{{\displaystyle {\sum }_{i=1}^K{\displaystyle {\sum }_{j=1}^K{u}_{i,j}}}}}$$

(18)

where K is the number of sample classes; $u_{i,j}$ is the number of samples belonging to the i-th class but predicted as the j-th class.

The macro F₁-score is defined as the mean of the F₁-score in each class, calculated by Equation (19). Unlike the accuracy, the macro F₁-score comprehensively considers the local generalization performance of the model in each class.

$$\mathrm{macro}\hspace{0.33em}{F}_1-\mathrm{score}={\displaystyle \frac{1}{K}}{\displaystyle {\sum }_{k=1}^K{F}_1-{\mathrm{score}}_k}$$

(19)

where $F_1-{\mathrm{score}}_k$ is the F₁-score in the k-th class, calculated by:

$$F_1-{\mathrm{score}}_k={\displaystyle \frac{2\cdot Pr{e}_k\cdot Re{c}_k}{Pr{e}_k+Re{c}_k}}$$

(20)

where $Pr{e}_k$ and $Re{c}_k$ are the precision and recall in the k-th class, respectively, calculated by:

$$Pr{e}_k={\displaystyle \frac{u_{k,k}}{{\displaystyle {\sum }_{i=1}^K{u}_{i,k}}}}$$

(21)

$$Re{c}_k={\displaystyle \frac{u_{k,k}}{{\displaystyle {\sum }_{j=1}^K{u}_{k,j}}}}$$

(22)

4.2. Evaluation Results

Out of 93 collected rockburst cases, 71 cases (accounting for 75%) were used as the training set to construct the model, and the rest 22 cases (accounting for 25%) were used as the test set to evaluate the model. In order to ensure the representativeness of the training set and the test set, stratified sampling was used to split the database. After that, the LOF-SMO was performed on the training set, as described in Section 2.1. On the one hand, class imbalance was eliminated. On the other hand, the number of training samples was increased from 71 to 104. Based on the balanced and extended training set, the C5.0DT-ERF was established for subsequent rockburst prediction, as described in Section 2.2. Particularly, the data were standardized by Min-Max standardization before modeling, described in Equation (23). In order to determine the optimal hyperparameters in the LOF-SMO and C5.0DT-ERF, including the size of the k-distance neighborhood in the LOF, the number of used nearest neighbors in the SMO, the minimum number of split samples in the C5.0DT, the maximum number of leaf nodes in the C5.0DT, and the number of C5.0 decision trees in the ERF, the WOA was implemented combined with a 10-fold cross-validation on the augmented training set, more technical information of which can be found in Section 2.3.

$$x^{\ast }={\displaystyle \frac{x-x_{\min }}{x_{\max }-x_{\min }}}$$

(23)

where $x^{\ast }$ is the standardized value; $x$ is the original value; $x_{\max }$ and $x_{\min }$ is the maximum and minimum, respectively.

Figure 11a shows the confusion matrix on the test set, which describes the predictive behavior of the model in detail. For six moderate rockburst cases and three strong rockburst cases, there were no wrong predictions. For eight none rockburst cases and five slight rockburst cases, one out of them was wrongly predicted in each. Based on the confusion matrix, the accuracy and macro F₁-score of the model were calculated by Equations (18) and (19), reaching 90.91% and 0.9141, respectively (Figure 11b).

Further, the local generalization performance of the model on low-intensity (none + slight) and high-intensity (moderate + strong) rockburst was dissected, as illustrated in Figure 12. For the former, the model achieved the F₁-score of 0.9600. For the latter, the model achieved the F₁-score of 0.9474. The results show that the model has a high prediction reliability for both low-intensity and high-intensity rockburst, with an F₁-score of more than 0.9.

In addition, referring to Figure 11a and Figure 12a, misjudged samples are generally from the low intensity category and predicted as high intensity. Compared to low-intensity rockburst, the model has superior sensitivity to high-intensity rockburst. From the perspective of engineering safety and risk management, this misjudgment is conducive to rockburst prevention, but at the same time, it may also cause an unnecessary increase in prevention costs.

4.3. Comparative Analysis

For ease of expression, the model proposed in this study is denoted as the LOF-SMO-C5.0DT-ERF. As the data-level improvement strategy, the LOF-SMO aims to eliminate class imbalance and overcome multicollinearity by increasing the data size. To probe the effectiveness of the LOF-SMO, a comparative analysis was conducted with the single C5.0DT-ERF, in which the LOF-SMO was not integrated. Referring to Figure 13, it can be found that the accuracy of the LOF-SMO-C5.0DT-ERF increased by 18.18% compared with the single C5.0DT-ERF, and, meanwhile, the macro F₁-score improved by 0.2105.

As the algorithm-level improvement strategy, the C5.0DT-ERF aggregates multiple standalone C5.0 decision trees through the extremely randomized strategy to accomplish the complicated prediction task. To verify the strength of the C5.0DT-ERF, a comparative analysis between the LOF-SMO-C5.0DT-ERF, LOF-SMO-MLP (multi-layer perceptron), and LOF-SMO--SVM (support vector machine) was conducted, as presented in Figure 14. Particularly for the MLP, a hidden layer was set, and the optimized hyperparameters included the number of neurons in the hidden layer and learning rate. And, for the SVM, the radial basis function was used, and the optimized hyperparameters included the penalty factor and kernel function width. Compared to the LOF-SMO-MLP and LOF-SMO-SVM, the LOF-SMO-C5.0DT-ERF showed an increase of 22.73% and 18.18% in accuracy, respectively. Regarding the macro F₁-score, the LOF-SMO-C5.0DT-ERF achieved an increase of 0.2442 and 0.2082, respectively.

4.4. Sensitivity Analysis

In this section, we analyzed the sensitivity of the model performance for each input parameter. This was achieved by removing each input parameter individually from the input system and then observing the change in model performance, as described in Equation (24). The more significant the performance change, the more sensitive the input parameter. The results showed that the LOF-SMO-C5.0DT-ERF is the most sensitive to cumulative energy, and the accuracy of the model decreased by 13.64% when cumulative energy was removed, followed by cumulative number, cumulative apparent volume, and changing rate of cumulative energy, which caused a 9.09% decrease in accuracy when they were removed. In contrast, the least sensitive input parameters are changing rate of cumulative number and changing rate of cumulative apparent volume, which each reduced the model’s accuracy by 4.55% when removed. Figure 15 summarizes the results of the sensitivity analysis.

$$\Delta Ac{c}_u=Ac{c}_S-Ac{c}_{u\notin S}$$

(24)

where $\Delta Ac{c}_u$ is the change in model accuracy when the input parameter u is removed from the input system; $Ac{c}_S$ is the accuracy of the model using the entire input system; $Ac{c}_{u\notin S}$ is the accuracy of the model using the input system without the input parameter u.

5. Conclusions

Rockburst is a common dynamic geological disaster in underground engineering that seriously threatens the safety of workers and equipment and leads to damage to the excavation. As an in situ real-time monitoring technology, microseismic monitoring is widely used to analyze the fracturing behavior of rock mass. This study aims to utilize machine learning to establish the quantitative mapping relationship between microseismic parameters and rockburst intensity so as to achieve the short-term evaluation of rockburst potential.

Class imbalance and multicollinearity are two common types of data defects in geotechnical and geological engineering. To eliminate their adverse effects on rockburst prediction, this study makes some improvements from two levels of data and algorithm. Correspondingly, a WOA-LOF-SMO-C5.0DT-ERF modeling framework is proposed. The results indicate that this model achieves an accuracy of 90.91% and a macro F₁-score of 0.9141. Additionally, the local F₁-scores on low-intensity and high-intensity rockburst are 0.9600 and 0.9474, respectively. Finally, through a comparative analysis with the single C5.0DT-ERF, LOF-SMO-MLP, and LOF-SMO-SVM models, the advantages of the proposed model were validated.

However, this study still has some limitations. First, the size of the used database is limited. It is well known that the generalization performance of machine learning models is closely related to database size. In general, the larger the database size, the better the generalization performance. Second, there is a reverse data structure in the used database. Normally, non-rockburst cases are more than rockburst cases, but here non-rockburst cases are fewer than rockburst cases. This may result in the application of the model producing more false predictions of non-rockburst cases. Third, this study only utilizes microseismic data to predict rockburst, neglecting the geological data (e.g., rock types, stress level, and the orientation of the fault and structural plane) and construction data (e.g., support conditions and the distance to the working face). Fourth, the rockburst type discussed here is dominated by stress rockburst, and there is a lack of in-depth analysis of fault rockburst and stroke rockburst. In future research, we will pay more attention to the establishment of the comprehensive database from the above four aspects, so as to further strengthen the performance of the model by improving the database quality.