# Genetic Programming-Based Prediction Model for Microseismic Data

Authors: Man Wang; Hongwei Zhou; Dongming Zhang; Yingwei Wang; Weihang Du; Beichen Yu
Journal: Geofluids
 (2022)

---

## Abstract

Microseismic monitoring is a rock breakdown monitoring technology, and it has become a major technical tool for underground disaster warning and prevention. However, the massive amount of data involved in multipoint monitoring raises the difficulty for microseismic research when the research object involves underground large-scale complex engineering on time and space scales. The risk predictions of rockburst based on microseismic parameters rely only on the experience of researchers and subjective factors by the human factor. The dependent variable cannot be related to the acoustic and energy signals in the form of a functional equation. Therefore, we collected a large amount of microseismic data obtained from the working faces of Shoushan Mine of Pingdingshan Coal Group in China and filtered the data, then we constructed a prediction model of microseismic data based on underground spatial three-dimensional coordinates using genetic programming (GP), which can realize real-time monitoring and disaster warning of microseismic signals. We established the magnitude and energy prediction formulas by GP and obtained the real-time evolution of each parameter of the optimal individual. The results showed that the actual seismic and energy data of the working faces exhibited good agreement with those predicted by the prediction formulas. The high energy and seismic distribution were caused by the mining stress at the edge of working faces due to the excavation and unloading effect of the surface, and the energy and seismic evolution predicted by GP could also show this phenomenon well.


---

## Body

## 1. Introduction

In the China’s 14th Five-Year Plan, the science and technology development will focus on frontier areas, i.e., deep earth, deep sea, aviation, and space technology [1]. Deep underground space is rich in coal, oil, and unconventional natural gas resources, but it is also under complex geological condition of high in situ stress, fluid pressure, and temperature with strong mining disturbance [2–4]. Therefore, the energy extraction is highly susceptible to rock burst and other dynamic hazards, which become a challenge to the efficient development of deep earth resources [5–9].Microseismic monitoring technology is a rock fracture monitoring technology and has become a major technical tool for early warning and prevention of subsurface hazards [10–12]. Based on the acoustic and energy signals transmitted by rock mass during breakdown, it analyzes the source characteristics and assesses the damage of rock mass, thus provides a theoretical basis for controlling rock dynamic hazards [13, 14]. Many scholars have evaluated the reservoir stability and production status based on microseismic monitoring [15–17]. Tang et al. [18] analyzed the distribution of rock burst in multipanel and multistope areas with simultaneous mining and designed sensor network spatial arrangement schemes with the optimized multiple microseismic monitoring systems. Maxwell et al. [19] monitored the deformation in Yibal oil field and other fields by identifying the microseismic signal characteristics and located the microseismic sources accordingly. Tang et al. [20] analyzed the rockburst mechanism in tunnel by microseismic monitoring and pointed out that the rockburst was induced by the interaction between rock and surrounding rock and deformation localization. Zhang et al. [21] used a microseismic monitoring system to obtain microseismic energy during mining and defined new rockburst prediction parameters. They established the relationship between microseismic energy and rock damage to effectively assess the stability of coal and rock mass.The acoustic and energy data obtained from microseismic monitoring provide an effective database for early warning and prevention and control of dynamic hazards. However, the large amount of data involved in multipoint monitoring raises the difficulty of microseismic research because the research objects generally involve the complex engineering on large temporal and spatial scales underground. Moreover, the risk prediction of rockburst based on microseismic parameters relies only on the experience of researchers and is greatly influenced by human subjective factors [22]. Therefore, using machine learning algorithm to process the microseismic data has become a research hotspot [23–25]. It can also ensure the accuracy and precision of the data. Wang et al. [26] proposed a self-optimization mechanism and introduced a machine learning algorithm called automatic detection network, which could monitor microseismic signals in real time. Tang et al. [27] proposed a new end-to-end training network architecture to automatically identify microseismic signals. The accuracy of intelligent microseismic monitoring was greatly improved by autonomously adjusting the weights of dependent variables without increasing the computational cost. Ma et al. [28] used a numerical method to simulate the progressive damage of coal seam floor during mining and constructed a microseismic monitoring system. They used virtual reality to realize 3D visualization of microseismic monitoring and analyzed the damage characteristics of the floor more accurately.Current scholars have used machine learning for effective processing of microseismic data. However, most of the studies often require much time to find the optimal or near-optimal network structure for microseismic signal monitoring and tuning, and they cannot relate the dependent variable to the acoustic and energy signals in the form of a functional equation [29, 30]. Genetic programming (GP), as an extension of genetic algorithm (GA), provides a new soft computing approach to solve this method [31–33]. It does not require users to adjust the weights and can automatically evolve the type, structure, and parameter of the required model. It explicitly derives predictive equation to describe the relationship between variables. Effective prediction of GP in geomechanics has been demonstrated [34–36]. Therefore, we collected a large amount of microseismic data obtained from the working faces of Shoushan Mine of Pingdingshan Coal Group in China and filtered the data, then we constructed a prediction model of microseismic data based on underground spatial three-dimensional coordinates using GP. The prediction formulas for microseismic data were obtained, which can realize real-time monitoring and disaster warning of microseismic data.


## 2. Geological Information

In this study, the microseismic data were taken from the 12070 and 12090 working faces of Shoushan Mine of Pingdingshan Coal Group in China. Shoushan Mine is located in the northeast of Pingdingshan City, with a shaft field of 6.1 km long in east-west direction and 4.4 km wide in north-south inclination. The average depth of the coal seam in the well field is 750 m, with an inclination angle of 8°~12°. The coal seams H15 and H16-17 are both coal seams with outburst risk. The working face of 12070 is located in the east of Shoushan Mine’s hexagonal mining area, with a strike length of 1 580 m and a design mining length of 240 m. The coal seam of H15-17 was mined with a mining height of 4.0 m. The comprehensive histogram of roof and floor slate of coal seam in 12070 working face is shown in Figure 1. The average burial depth of 12090 working face is 455 m, and the average coal thickness is 4.5 m. The surrounding rock is generally monoclinic structure with a slightly westward inclination to the north. The dip angle is around 10°. The original gas content is 4.32 m3/t, and the Protodyakonov coefficient of the coal seam is 0.17. The gas pressure was measured as 0.21 MPa. According to the available geological data, there is no obvious structure in this working face, but the top and bottom plates are undulating. It has a certain impact on the excavation. The lithology of roof and floor slate of coal seam in 12090 working face is shown in Table 1. In this study, 7851 microseismic data near 12070 and 12090 working faces are selected for this study, and their magnitude distribution is from 0.41 to 6.86. The depth distribution of microseismic data is from 0 to -1400 m. The rule of sensor recording data is that when any sensor received a signal exceeding the set threshold, the data from 500 ms before to 1000 ms after the reception of this signal was intercepted by all sensors as a microseismic data. The next step is to filter the energy and magnitude data obtained from the field microseismic data monitoring to exclude the interference of noise.Figure 1
Comprehensive histogram of roof and floor slate of coal seam in 12070 working face.Table 1
Lithology of roof and floor slate of coal seam in 12090 working face.
TitleRockThickness/mLithologyRoofMain roofDazhan sandstone3~12.1Deep gray, gray black, quartz mainly argillaceous cementation thick layer; local microsandy mudstone.Immediate roofSandy mudstone4~10.9; average 7.4Deep gray, gray black, medium-thick stratiform mica flakes, and carbonaceous local sandy mudstones containing fossil plants.False roofMud or sandy mudstone0.5~3; average 1.7Soft, broken, and easy to fall.FloorImmediate floorMud or sandy mudstone0~2.6Black ash, black, thin-layer carbonaceous, and plant fossils, undulating. It may have an impact on recovery.Main floorSandy mudstone or fine sandstone3.1~8.3Gray black, gray, pyrite-bearing nodules and argillaceous bands visible in large quantities of pyrite nodules and siderite nodules.


## 3. Filtering of Sample Data

The original microseismic waveform obtained from the sensor contained large background noise, so the original waveform needed to be filtered. Filtering is an operation that filters out specific wave frequencies (high or low-frequency) from the signal to extract the desired signal [37]. Generally, high-frequency signals in images are the points in the image that have significant grayscale changes and present contours or noise, while low-frequency signals show flat and insignificant grayscale changes in the image. According to the high and low frequencies of the image, corresponding high and low pass filters can be designed. High pass filtering detects sharp and obvious changes in the image, while low pass filtering makes the image smooth and filters out the noise in the image. The filter is often designed as a discrete-time system. So the system input is a time series of xt, and the output is yt. yt versus xt can be expressed as follows [38]:
(1)yt=ht⊗xt+nt,where ht denotes the response between the input and output, ⊗ denotes the convolution sign, and nt denotes the additive noise caused by the inverse filtering to achieve the best recovery of smoothing metric. nt is not correlated with xt. Fourier transform of the above equation can be obtained:
(2)Yjω=Hjω⊗Xjω+Njω,where Hjω is the transfer function of the filter, which reflects the frequency characteristics of the filter. Njω is the power spectrum of the additive noise.Low pass filtering involves the linear mean filter, Gaussian filter, nonlinear bilateral filter, and median filter; high pass filtering has various filters based on Canny and Sobel operators. Low and high pass filtering algorithms often contradict each other, so it is often necessary to reduce noise by low pass filtering before edge detection and adjust the parameters to deal with more image noise without losing high-frequency edges [39].The ideal low pass filter model can be expressed as follows:(3)Hx,y=1,Dx,y≤D00,Dx,y≥D0,where D0 is the passband radius, and Dx,y is the distance to the center of the spectrum, also known as the Euclidean distance; it can be calculated as follows: Dx,y=x−X/22+y−Y/22, where X and Y denote the size of the spectrum image, and (X/2,Y/2) denote the center of the spectrum.The Gaussian low pass filter can be expressed as follows:(4)Hx,y=e−D2x,y2D02.After obtaining the low pass filter expression, the high pass filter expression can be construct by subtracting the low pass filter model by 1.For better noise removal, the response functionHjω needs to be continuously adjusted. Therefore, the convolution function Gjω needs to be found to further improve the input signal Xjω according to the output:
(5)X^jω=Gjω∗Yjω.The filter function in the frequency field can be expressed as follows:(6)Gf=H∗fRfHf2Rf+Nf=H∗fHf2+Nf/Rf,where Gf and Hf denote the Fourier transform of ht and gt in the frequency field, and Rf and Nf denote the power spectrum of the input function xt and additive noise nt in the frequency field, respectively. ∗ denotes the conjugate complex number. Equation (6) can be further rewritten as follows:
(7)Gf=1HfHf2Hf2+Nf/Rf.where Nf/Rf denotes the signal-to-noise ratio. When the signal-to-noise ratio tends to infinity, i.e., the noise closes to zero, the above equation is equal to 1, so the filtering is simplified as an inverse filtering process. However, when the noise increases, the signal-to-noise ratio decreases, and the values of the above parameters decrease as well. This suggested that the bandpass frequency of the filtering depends on the signal-to-noise ratio. Substituting Equation (7) into Equation (5), the input signal X^f in the frequency field can be obtained:
(8)X^f=1HfHf2Hf2+Nf/Rf∗Yf.After obtainingX^f, the deconvolution result x^f can be obtained by the inverse Fourier variation.The original microseismic waveform obtained from the sensor was filtered, and the microseismic energy and magnitude data obtained were shown in Figures2 and 3. The filtered data will be used to build a prediction model by GP for the energy and magnitude data in next section.Figure 2
Energy and magnitude of 12070 working face coal seam.
(a)
power(b)
magnitudeFigure 3
Energy and magnitude of 12090 working face coal seam.
(a)
power(b)
magnitude


## 4. Energy and Magnitude Prediction Modeling

GP, as an improvement of GA [40], has been widely used in geotechnical engineering for reservoir mechanics prediction and simulation [41–43]. The principle is based on the survival of the fittest of the genetic law, which simulates the reproduction and evolution of populations (e.g., replication, crossover, and mutation) to generate optimal individuals for given condition. In contrast to GA, GP utilizes tree structure codes to represent individuals in a population instead of binary code lines of GA [35, 44]. The tree structure individuals consist of predictor variables, mathematical symbols, and numbers, and the expressions of output variables can be constructed under fixed output rules. Genetic operators (replication, crossover and mutation, etc.) are used on tree-structured individuals in population to create new offspring programs. It can generate the best individual, also known as the individual with the highest fitness, i.e., the equation with the best prediction, under the rule of the survival of the fittest [45]. The sum of the absolute error between the measured and predicted values of an individual (Sumer) is used as a measure of fitness in GP. This means that a lower Sumer corresponds to a higher fitness, and the Sumer equation can be expressed as follows:
(9)Sumer=∑i=1Nmsmi−mspi,where msmi is the microseismic data obtained through the field; mspi is the microseismic data predicted by GP. N is the overall number of samples.Table2 shows the parameters used in the GP modeling. In this study, the three-dimensional coordinates (x, y, and z) of the subsurface space were used as predictor variables, and the energy and magnitude data were used as output variables to be predicted using the coordinates. In the initialization of the microseismic population, the ramped half-and-half method was used to generate microseismic individuals with different breadth and dimensions to enhance diversity [36]. In the selection of microseismic individuals, the tournament method was used to randomly select a certain number of highly adapted individuals from the microseismic population as parents using high fitness as the criterion, and then create offspring of microseismic individuals using genetic operators including replication, crossover, and mutation [36]. The highest fitness individuals in each generation are directly entered into the offspring by the replication operator, and the number of replications can be specified by the user. For the remaining positions of next generation individuals, the microseismic individuals are genetically manipulated according to the programmed crossover and mutation probabilities.Table 2
Parameters of the GP model.
GP parameterValuesTerminal setx, y, zFunction set+, -, ×,/, sin, cos, tan, exp, sqrt, ln, powerPrediction setPower and magnitudeFitness functionSumerMaximum generation100Population size100Tree initializationRamped half-and-halfSelection methodTournamentGenetic operatorsMutation, crossover, replicationMutation rateDynamic (initial value 0.5)Crossover rateDynamic (initial value 0.5)Reproduction rateDynamic (initial value 0)Termination criterionMaximum generationSubtree crossover is the selection of random crossover points (nodes) of microseismic parents and replacement to create new microseismic individuals, and the crossover is shown in Figure4(a). Mutation is the selection of a random mutation node of a microseismic individual and replacing the subtree of that node using the system randomly generated microseismic individuals, and the mutation is shown in Figure 4(b). The initial values of crossover and mutation probabilities are the default values in the system settings. The initial values of mutation and crossover probabilities are set to 0.5, and they are adjusted at any time during the population iteration to reflect the operator performance. An increase in the operator probability indicates that the produced microseismic individuals are more adaptive than the previous generation; at the same time, it decreases the probability of producing individuals that are poorly adapted to the environment. The microseismic individuals with the highest fitness obtained during the final run of the algorithm are the predicted results of GP.Figure 4
Schematic diagram of individual (a) crossover and (b) mutation.
(a)(b)The energy and magnitude data are substituted into the GP separately, and the following prediction curves can be obtained. The functions of the algorithm running of population iteration can be presented separately, taking the magnitude prediction of working face 12090 as an example. Figure5 shows the dynamic evolution of the crossover and mutation probabilities during the algorithm run, and the cumulative occurrence of crossover, mutation, and replication operators. When the individuals with high fitness of population appeared, the replication operator would be introduced to replicate the high fitness individuals directly to the next generation to increase the average fitness of the population, and the occurrence of the replication operator in the last generation is 87. It can be found from the figure that the crossover operator has a higher ability than the mutation operator to produce highly adapted individuals in the early iterations of the population. However, as the population evolves, the fitness stabilizes, and the mutation ability of the magnitude individuals needs to be improved to create better individuals. The ability of mutation operator to create the best individuals was higher than that of crossover in the late evolutionary stage, with 474 and 525 occurrences of crossover and mutation operators in the last generation, respectively.Figure 5
Evolution of operator probabilities and cumulative occurrences.Figure6 shows the maximum, optimum, and average sizes and depths of the population individuals in the 12090-magnitude modeling. As the population multiplies, the size and depth of the population individuals increase. The magnitude prediction expressions, represented by the magnitude individuals, also gradually increase in complexity but also gradually stabilize as the generation increases. The purple curve in the figure represents the average fill rate of population individuals, which is closely related to population diversity. The average filling rate represents the filling degree of each individual, and the average filling rate of unevenly developed individuals is low. At the beginning of the program, the average fill rate presents a sudden decrease, but the decrease gradually tends to be balanced. The lower average fill rate implies higher population diversity, which shows uneven development of individuals. This phenomenon is related to the ramped half-and-half method invoked to generate microseismic individuals with different breadth and dimensionality. Figure 7 shows the fitness, number of nodes, and depth evolution of the optimal individuals in each generation. Since we used the sum of absolute error (Sumer) between the measured and predicted values of individuals as a measurement of fitness, it can be found that the fitness of the optimal individuals gradually increases, showing a continuous decrease in Sumer. The best individual appeared at 90th generation of population evolution, and the fitness of the best individual was 212.20, the depth of the tree was 14, and the number of nodes was 24.Figure 6
Evolution of tree size and depth in population.Figure 7
Fitness, depth, and nodes evolution of the optimal individual.The tree expression of the magnitude individual with the highest fitness generated at the end of the population iteration for a given 3D coordinate environment is shown in Figure8. To obtain the prediction formula for the magnitude, we traversed each node starting from the bottom and left node, based on the rule of bottom to top and left to right, until the top of the tree. The final generated prediction formulas are shown in Table 3. Table 3 shows the depth, number of nodes, iterative generations, fitness (Sumer), R2, and the final prediction formula for the energy and magnitude optimal individuals of 12070 and 12090 working faces. From the equations, the coordinates y and z have a greater influence on energy and magnitude and show a significant role in the prediction. In contrast, the x-direction does not show a significant difference between energy and magnitude, so the x-coordinate does not appear in the prediction formula for magnitude and energy. Substituting the coordinates into the prediction formula, the obtained energy and magnitude data for the two working surfaces are shown in Figures 9 and 10. Due to the excavation and unloading effect on the working face, the mining stress causes the high energy and magnitude at the edge of the working face, and the energy and magnitude evolution predicted by GP can also describe this phenomenon well. Comparing Figure 3, Figure 4, and Figure 9 with Figure 10, it can be found that the actual magnitude and energy data of the working face fit well with the predicted data, and R2 does not show an accurate fit due to the large amount of data obtained, but it is within the acceptable range. This implied that GP shows good performance in predicting the energy and magnitude signals in the field.Figure 8
Optimal tree expression generated by GP.Table 3
Prediction parameters of rock strength obtained by GP.
Modeling objectDepthNodesGenerationFitnessR2Equation12070 power12611001.97e80.56power=10.97−Y+30Z+21062.42+sin204Y−137059+lncos6+lnY·15Y−35.1012070 magnitude819801665.830.74magnitude=coscosln43+Z+0.80394Y+0.8912090 power26531001.42e70.78power=27×30×ln24×ln116−24×ε−37−Z−12−Z+30ε=24×ln52−24×ln−27−Z−37−37−Z−2712090 magnitude1424100212.320.88magnitude=coscoscoscoscoslncos1+coslncos19Y−6+coscoslnZ−6Figure 9
Prediction of power and magnitude of 12070 working face coal seam.
(a)
power(b)
magnitudeFigure 10
Prediction of power and magnitude of 12090 working face coal seam.
(a)
power(b)
magnitude


## 5. Conclusion

In this study, we established the magnitude and energy prediction formulas by GP and obtained the real-time evolution of each parameter of the optimal individual based on a large amount of energy and magnitude data from the 12070 and 12090 working faces of Shoushan Mine of Pingdingshan Coal Group in China. The data were filtered, and a prediction model of microseismic data was constructed based on underground spatial three-dimensional coordinates. The main conclusions of this study are as follows:(i)
Since the original microseismic waveform contains large background noise, the original waveform was filtered. The filter function was designed to convolve the functionGjω in the frequency field, and the input signal Xjω can be continuously improved according to the output, and x^f after deconvolution can be obtained by the inverse Fourier variation, and then the noise elimination was achieved(ii)
The magnitude and energy prediction formulas are established by GP. The real-time evolution of each parameter of the optimal individual is obtained. The actual magnitude and energy data of the working face show good agreement with those predicted by the GP prediction formulas. Due to the excavation and unloading effect of working face, the mining stress causes the high energy and magnitude at the edge of the surface, and the energy and magnitude evolution predicted by GP can also show this phenomenon well



---

*Source: 2525923-2022-10-12.xml*