Abstract

In order to mitigate economic and safety risks during mine life, a microseismic monitoring system is installed in a number of underground mines. The basic step for successfully analyzing those microseismic data is the correct detection of various event types, especially the rock mass rupture events. The visual scanning process is a time-consuming task and requires experience. Therefore, here we present a new method for automatic classification of microseismic signals based on the Gaussian Mixture Model-Hidden Markov Model (GMM-HMM) by using only Mel-frequency cepstral coefficient (MFCC) features extracted from the waveform. The detailed implementation of our proposed method is described. The performance of this method is tested by its application to microseismic events selected from the Dongguashan Copper Mine (China). A dataset that contains a representative set of different microseismic events including rock mass rupture, blasting vibration, mechanical drilling, and electromagnetic noise is collected for training and testing. The results show that our proposed method obtains an accuracy of 92.46%, which demonstrates the effectiveness of the method for automatic classification of microseismic data in underground mines.

1. Introduction

Microseismic monitoring is becoming a popular technology with wide and successful applications in petroleum, mining, and geotechnical engineering [13]. In general, the microseismicity is recorded by accelerometers and/or geophones buried in the surrounding rock. Thus, by processing the recorded waveforms, microseismic monitoring can provide important insight into a rock mass and quantify where a certain magnitude of induced rock fracturing is occurring within the volume [4]. With decades of applications, microseismic monitoring plays a vital role in generating valuable information and mitigating economic and safety risks during mine life [57]. The system detects all types of vibration events (e.g., rock mass rupture, blasting vibration, mechanical drilling, and electromagnetic noise) in its monitored area. These records are then processed to evaluate the rock mass stability. However, identification of suspicious microseismic events is the first key step of processing microseismic data, which is usually done by experienced analysts through visual scanning of waveforms manually. Thus, microseismic records classification is a time-consuming and tedious task, needs experience, and may suffer from the subjective view of the observers [8]. For these reasons, an automatic technique of identification of event types is in great request.

Throughout the years, many automatic classification methods have been proposed to address this problem in seismic and microseismic fields [914]. For example, Scarpetta et al. use neural networks to distinguish volcano-tectonic earthquakes from local seismic signals [15]. Provost et al. propose a random forest supervised classier to identify the type of seismic sources based on 71 seismic attributes [16]. Ruano et al. build a classifier using support vector machines, aiming at distinguishing local and regional earthquakes and explosions from the other possibilities in earthquake early-warning system [17]. Shang et al. propose a hybrid technique based on principal component analysis and artificial neural networks (PCA-ANN) to discriminate between microseismic events and quarry blasts [18]. The PCA-ANN is trained on a dataset with 1600 events, and 22 source parameters are extracted from each event, such as corner frequency, seismic moment, energy, source radius, and static stress drop.

The results of these works are very encouraging because they demonstrate an alternative way to do the tedious scanning work. However, there is still room for further improvement. Most of these classifiers are trained on a large number of parameters [1821], which are acquired through experienced processing. In other words, these algorithms cannot classify an event unless basic processing procedures (e.g., P-wave arrival picking and epicenter location) are done. Thus, these algorithms are not suitable for real-time processing and early-warning system. In this article, we focus on presenting an automatic classification system that requires a minimum amount of training data while enabling to recognize highly variable event patterns. Mel-frequency cepstral coefficient (MFCC) is a widely used feature that has been successfully applied in speech recognition and volcano classification [22]. Hidden Markov Model (HMM) is a powerful tool in modeling any time-varying series [22, 23]. In particular, microseismic records can also be modeled as a time sequence of different microseismic events. Gaussian mixtures are capable of clustering data into different groups as a collection of multinomial Gaussian distributions. Each microseismic signal can be devised as a collection of multinomial distribution and HMM can model the intraslice dependencies between each time period. Therefore, in this research, we propose to utilize Gaussian Mixture Model-based HMM (GMM-HMM) for microseismic signal classification using only MFCC features extracted by waveforms.

The rest of the paper is organized as follows. Section 2 presents the implementation details of automatic classification of microseismic signals based on MFCC and GMM-HMM. In Section 3, we test our proposed method using the field data recorded in an underground copper mine. Finally, Section 4 summarizes the main conclusions.

2. Methods

2.1. MFCC

The first step of the classification process is feature extraction, which converts the microseismic waveform to a parametric representation with less redundant information. Feature extraction consists of choosing those features which are most effective for preserving class separability. Davis and Mermelstein [24] first proposed the MFCC, which is a representation of the short-time energy spectrum of the signal waveform; it is obtained by projecting the logarithmic power spectrum of the microseismic signal onto the nonlinear Mel scale by linear cosine transformation. The transformation relationship between the Mel scale and the frequency is as follows:where f is the frequency of the signal. The MFCC is well suited to compensate for signal distortion.

The steps for extracting the MFCC are as follows:(1)High-pass filtering: Passing the waveform of microseismic signal through a high-pass filter can effectively enhance the high-frequency portion, reducing the spectral fluctuation of the signal waveform and allowing any band of the spectrum, regardless of high or low frequencies, to be obtained based on a similar signal-to-noise ratio.(2)Waveform framing: The microseismic waveform is segmented every N sampling points to form a new waveform unit, called a frame. According to the requirements of feature extraction and the length of the signal, the value of N is generally selected as 256 or 512. In addition, the front and back frames of each frame are intersected with a small portion thereof to avoid an excessive difference between consecutive frames.(3)Adding the hamming window: To increase the continuity between each frame and its adjacent frames after the microseismic waveform is framed, each frame of the waveform is multiplied by the Hamming window. Assuming that the microseismic signal is , multiplying the signal by the Hamming window giveswhere N is the number of frames of the framed microseismic signal.(4)Fast Fourier transform: The difference in energy distribution can represent the features of different signals, so the microseismic signals are converted into an energy spectrum in the frequency domain for comparison. After the microseismic signals are continuously overlapped (50% overlap between successive windows) and framed, the fast Fourier transform is performed on each of the decomposed signal frames to calculate an energy spectrum in the frequency domain.(5)Triangular band-pass filtering: Using a series of triangular band-pass filters Hm(k), the energy spectrum obtained by fast Fourier transformation is converted to the Mel scale to obtain a set of coefficients m1, m2, …, mP. The series of filters are a series of triangular windows that are spaced evenly with overlapping on the Mel-frequency axis.(6)Calculation of the logarithmic energy spectrum: Calculate each filter bank, and take the logarithm of the result. The obtained value is the logarithmic energy, and the logarithmic power spectrum of the corresponding frequency band can be obtained. The calculation formula is as follows:where X(k) is the energy spectrum of the microseismic signal and Hm(k) is the filter bank, where m = 1, 2, …, P, P is the number of filters and s(m) is the logarithmic energy.(7)Discrete cosine transform: The discrete cosine transform is used to transform the spectrum from the frequency domain to the time domain. The result is the standard MFCC. The mathematical formula for calculating the cepstral coefficient is as follows:where n is the number of frames calculated, with , and m is the number of Mel-frequency cepstral coefficients, with .

As the lower order coefficients contain most of the information about the overall spectral shape of the source-filter transfer function, it has become customary in many signal applications to select the first 12 MFCCs because they are considered to carry enough discriminative information in the context of various classification tasks. Consequently, we use 12 features calculated by equation (4) and 12 difference cepstral parameters to form a 24-dimensional feature parameter vector in order to improve the classification performance. The difference cepstral parameters obtained by the difference operation show the variation of the original 12 Mel-frequency cepstral coefficients in the time domain. The method for calculating the difference cepstral parameter Dt(n) is as follows:where Dt(n) represents the t-th first-order differential cepstral parameter of the Mel-frequency cepstral coefficient C(n) calculated in the n-th frame of the signal. Θ denotes the time difference of the first derivative in the expression; generally, Θ = 2, and .

Following the above steps, the 24-dimensional feature parameter vector can be successfully extracted from the original complex microseismic signal.

2.2. GMM-HMM

To realize the automatic classification of microseismic signals, the basic procedure is to extract the features of the waveforms by using algorithms and then use these features in combination with machine learning. The classification system presented in this paper is based on GMM-HMM.

The HMM is a probabilistic model of time series. An HMM can typically be represented by five parameters: λ = (N, M, π, A, B), where N is the size of the Markov state chain in the HMM and is a fixed value in actual use. Let N states be θ1, θ2, …, θN; the state at time n is then qn, where qn  ∈  (θ1, θ2, …, θN). M is the number of observations that may correspond to each state in the Markov state chain. Let M observations be V1, V2, …, VM; then, the observed value at time n is on, where on ∈ (V1, V2, …, VM). π is the initial state probability distribution, and π ∈ (π1, π2, …, πN), where πi = P(q1 = θi) and 1 ≤ i ≤ N. The parameter A is the state transition probability matrix: A = [aij]N×N, where aij = P(qn+k = θj | qn = θi) and 1 ≤ i, j ≤ N, indicating that at any time n, if the state is θi, then the probability of the state at the next time instant is θj. B is the observed value probability matrix: B = [bij]N×M, where bij = P(on = Vk | qn = θi) with 1 ≤ I ≤ N and 1 ≤ j ≤ M, indicating the probability that the observed value Vk is acquired at any time n if the state is θi.

To better identify complex microseismic signals, the GMM-HMM is constructed with the probability density function of the observed values by using the Gaussian Mixture Model based on the original HMM technique, and bjk is modified to the Gaussian distribution probability density function between the current state and observation, that is,where μjm is the mean, Ujm is the variance, and cjm is the Gaussian distribution weight; thus, a GMM-HMM is constructed.

2.3. Classification Strategy: Implementation Details

Our proposed automatic classification method is divided into three major steps:(1)Feature extraction of microseismic signals: In this paper, according to the extraction process in Section 2.1, the microseismic signals are framed, and the Mel-frequency cepstral coefficients of each frame are calculated to obtain 12 Mel-frequency cepstral coefficient values. Then, the differences between the 12 Mel-frequency cepstral coefficients are calculated. As a result, a 24-dimensional feature vector is extracted from a microseismic signal to provide a database for automatic classifier.(2)Training of an automatic classification model: The microseismic signals marked in a certain period of time are selected to extract the feature vectors by the MFCC, and the feature vectors are used as training data set for the corresponding event type. After the model parameters N, M, π, A, and B are obtained by iterations, the GMM-HMM classifier is constructed. The flowchart for microseismic signal feature extraction and classification is shown in Figure 1.

(3)Application: According to the classification model for the four microseismic events (rock mass rupture, blasting vibration, rig drilling, and electromagnetic interference) obtained from training, the Mel-frequency cepstral coefficients are calculated for the microseismic events generated in the mine production activities. The obtained feature vectors are substituted into the classification model to calculate the probability values of the corresponding GMM-HMM. The type of event corresponding to the model with the largest probability value is taken as the type of the microseismic event to be evaluated, thus realizing the automatic classification of the microseismic signals in the production activities. The flowchart for application of the automatic classifier is shown in Figure 2.

3. Field Application

3.1. Microseismic Monitoring System in Dongguashan Copper Mine

Dongguashan Copper Mine is located in the east of Tongling, Anhui Province, China, and the deposit is in the deep part of the Shizishan copper mining area. This mine is the first hard rock metal mine in China, with a mining depth of up to several kilometers and prominent susceptibility to rock burst. The mine has the characteristics of large deposits and high burial depth. The design mining mode is high-yield enhanced recovery that is characterized by the large number of mining panels and stopes, also leading to a wide distribution of the mining stopes. In addition, the multistope parallel mining mode leads to an excessively fast fitting speed and an excessive number of rockburst events with a wide distribution. Therefore, the microseismic monitoring for the prediction and early warning of rockburst activities during the mining process was introduced. In October 2017, the microseismic monitoring system was officially operated at the mine. The system monitored the production activity of Dongguashan in real time, and extensive research on microseismic activities based on monitored data was carried out.

3.2. Training of the Classification Model Based on Microseismic Data of Dongguashan Copper Mine

In this study, the microseismic data of Dongguashan Copper Mine before December 13, 2017, are selected as the training data sets. Four types of events, rock mass rupture events, blasting vibration events, rig drilling events, and electromagnetic interference events, are manually labelled by experienced experts. The total number of the events is 1400, with each type of 350 events.

The training process is performed according to the flow shown in Figure 1. First, the MFCC is extracted for each event waveform. Figures 36 show the waveforms randomly selected from each event type and its corresponding MFCC values. Then, the 24-dimensional feature vector composed of the manually labelled event types and the MFCC values is used as input data to determine the GMM-HMM parameters N, M, π, A, and B. Finally, the optimal parameters of the model are obtained by an iterative algorithm, thus obtaining a trained automatic classification model for the microseismic signals.

3.3. Results and Discussions

In this study, the microseismic monitoring system of Dongguashan Copper Mine is selected to test the GMM-HMM classification model using the microseismic events recorded between December 13, 2017, and January 17, 2018 (due to system maintenance, no data were recorded by the monitoring system on December 23–24, 2017). All events have been manually confirmed and type-marked as a basis for evaluating the classification accuracy of our method. Figure 7 shows the number of various types of microseismic events monitored daily during the abovementioned period at the mine.

In this paper, according to the flow of application of the automatic classifier in Figure 2, the types of events shown in Figure 7 are classified using our proposed method, and the accuracy of the model is tested in comparison with the types of manually labelled events.

For the constructed automatic classification model, we need to determine the relevant model parameters according to the actual situation of the mine: (1) the number of Gaussian equations in the hybrid Gaussian method Q; (2) the number of Markov chain states N; and (3) the number of observations M in the Markov chain. Among them, the observed value M is the dimension of the MFCC in this paper, that is, M = 24. We choose the number of Gaussian equations Q by testing the identification accuracy under different numbers. As shown in Figure 8, when Q ≥ 3, a high identification accuracy can be achieved, but the use of more Gaussian equations increases the calculation cost of the model, so Q = 4 is used in this paper. Similarly, we choose the final number of states N by comparing the identification accuracies with different trial N values. As shown in Figure 9, the identification accuracy increases as N becomes larger, but the subsequent increase rate is small. Considering the computational efficiency, N = 6 is used in the example.

The actual microseismic signals are classified based on the selected model parameters. The identification results of each event are shown in Figure 10.

Figure 7 shows that there are a total of 981 events between December 13, 2017, and January 17, 2018. Among them, there are 467 rock mass rupture events, 138 blasting vibration events, 108 rig drilling events, and 268 electromagnetic interference events. As shown in Figure 10, the numbers of incorrectly identified events are 26, 32, 2, and 14 for rock failure events, blasting vibration events, rig drilling events, and electromagnetic interference events, respectively. Therefore, the model for automatic identification and classification of microseismic events based on MFCC and GMM-HMM has a high identification accuracy of 92.46% as shown in the test. The results show that the automatic classification method realizes the high timeliness of the microseismic monitoring system in mines and provides a technical support for real-time analysis by the microseismic monitoring system. It should be noted that the identification accuracy for blasting vibration events is lower than the others. Among all incorrectly identified blasts, 22 events are identified to rock mass rupture events, which is consistent with fact that there are many similarities between blast and high-energy rock mass rupture. Additionally, due to the relatively small data used for our training, the accuracy will be improved by a large dataset collected for a long time.

4. Conclusions

In this paper, we present a new method for automatic classification of microseismic signals based on GMM-HMM by using only MFCC features extracted from the waveform. The performance of this algorithm has been tested by its application to microseismic events selected from Dongguashan Copper Mine (China). Using optimal parameters, the results show that our method obtains an accuracy of 92.46%, which outperforms many other algorithms. The input of the method is only the waveform, which means it is suitable for real-time processing. Requiring a minimum amount of preparation time and workload, the method has several advantages over classical techniques. Until now, the use of microseismic warning systems for mining disasters is mainly limited by the incorrect identification of the signal source. By using the suggested approach, this problem can be overcome, leading to an automatic detection of rock mass rupture and other sources.

Data Availability

The data used to support the findings of this study were provided by Dongguashan Copper Mine under license and so cannot be made freely available. Access to these data will be considered by the author upon request, with permission of Dongguashan Copper Mine.

Conflicts of Interest

The authors state that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

This work was supported by the National Key R&D Program of China (No. 2017YFC0602905).