Quantitative prediction method of three spatial and temporal elements of microseismic events in rock burst mines

Abstract

A quantitative prediction method for the three elements of spatiotemporal intensity of microseismic events in rock burst mines, comprising the following steps: 1. Importing the spatiotemporal intensity historical data of microseismic events. 2. Preprocessing the data of step 1. 3. Using the data of step 2 to create a training set and a test set. 4. Constructing a microseismic prediction model. 5. Inputting the training set data of step 3 to the prediction model training of step 4. 6. Inputting the test set data of step 3 to the prediction model of step 4 for prediction, and evaluating the prediction effect of the model. 7. Inputting the spatiotemporal intensity historical data of microseismic events and the current data to the qualified prediction model in step 6 to realize the spatiotemporal intensity prediction of microseismic events. The present invention can realize the prediction of the three elements of time, spatial position and intensity of microseismic events through the above steps, and has the advantage of quantitative prediction.

Description

Technical Field

[0001] The technical field to which the present invention belongs is the field of microseismic prediction for coal mine safety, and specifically provides a quantitative prediction method for the temporal and spatial intensity three elements of rock burst mine microseismic events. Background Art

[0002] Prediction is the prerequisite and foundation for prevention and control of destructive mine tremors and rock burst disasters, and plays a crucial role. Accurate quantitative prediction is a measure to mitigate the severity of these disasters. Currently, regional microseismic monitoring systems have been established in all mines with rock burst in my country. These systems continuously monitor the occurrence of microseismic events caused by rock failures underground 24 hours a day. The data recorded by these systems is limited to the calculation of the three elements of time, spatial location, and energy of microseismic events that have already occurred. Traditional prediction methods are mostly qualitative, such as assessing rock burst hazards based on coal mine geological conditions and mining conditions, dividing hazardous areas and hazard levels, or predicting the likelihood of large-energy microseismic events based on the frequency and energy trends of monitored microseismic events. Quantitative prediction is difficult to achieve. With the deepening of coal mining in my country, the occurrence of destructive mine tremors and rock burst disasters is becoming increasingly serious. Therefore, quantitative prediction of the three elements of time, spatial location, and energy of microseismic events in rock burst mines is of great significance for the prediction of destructive mine tremors and rock burst. The present invention proposes a method for realizing prediction driven by microseismic big data. By deeply learning the time, space and energy evolution laws of a large amount of information on microseismic events that have occurred, the three elements of time, space and energy of future microseismic events are predicted, especially the induction of destructive mining earthquakes and rock burst disasters, so as to achieve the purpose of reporting earthquakes with earthquakes and rock bursts with earthquakes, and taking emergency measures in advance to reduce disaster losses. Summary of the Invention

[0003] The present invention extracts the characteristics of microseismic historical data based on a deep learning algorithm and makes quantitative predictions on the occurrence time, spatial position coordinates and released energy of microseismic events.

[0004] The present invention is achieved through the following technical solutions:

[0005] The quantitative prediction method of the three elements of temporal and spatial intensity of rock burst mine microseismic events is as follows:

[0006] 1) Import historical microseismic data, including earthquake time, source coordinates, and microseismic energy;

[0007] 2) Preprocessing of microseismic data;

[0008] 3) Using the processed microseismic data to create microseismic training datasets and microseismic test datasets;

[0009] 4) Construct a BWF-BiGRU-BO microseismic prediction model based on Butterworth filter (BWF), bidirectional gated recurrent unit (BiGRU), and Bayesian optimization (BO);

[0010] 5) Input the microseismic training set data into the prediction model in step 4 to train it;

[0011] 6) Input the microseismic test data set into the model trained in step 5) and use the evaluation index to verify the effect of the model in step 5. If the effect meets the expectations, proceed to step 7). Otherwise, increase the number of iterations and retrain the model through step 5).

[0012] 7) Input the actual microseismic sequence data into the prediction model trained in step 5 for quantitative prediction.

[0013] In the step 1), the historical microseismic data are specifically as follows: when a microseismic event occurs, the rock mass fractures, energy is released, and vibration waves are generated. After multiple sensors collect the vibration waves, the program calculates three elements: time series data of the earthquake occurrence time, the source coordinates, and the microseismic energy; the earthquake occurrence time is composed of year / month / day hour: minute: second, and each moment corresponds to a microseismic event; the source coordinates are composed of three components of X, Y, and Z, which are the relative offset of the coordinate origin defined by the mining area, and the unit is meter; the microseismic energy represents the amount of energy released by this microseismic event, and the unit is joule.

[0014] The pretreatment method in step 2) is specifically:

[0015] 2.1) Using the difference algorithm to process the earthquake time series data into time interval series data I;

[0016] I=t next -t now (1)

[0017] where t next is the time when the next microseismic event occurs, t now is the time when the current microseismic event occurs;

[0018] 2.2) Analyze the autocorrelation of microseismic data. Calculate the autocorrelation of the earthquake occurrence time, source coordinates, and microseismic energy data. Specifically:

[0019] Input the time series data of the time interval obtained by formula 1 into the autocorrelation function for calculation to obtain the autocorrelation analysis graph of the time interval interval; input the time series data of coordinates X, coordinates Y, coordinates Z and energy into the autocorrelation function to obtain the autocorrelation analysis graphs corresponding to coordinates X, coordinates Y, coordinates Z and energy respectively;

[0020] The autocorrelation function is independent of any specific time t and is a function of the time interval h, expressed as ρh express:

[0021]

[0022] Among them, h is the time interval, cov(X t ,X t+h ) represents the variable sequence X t ,X t+h The covariance of t σ t+h Represents standard deviation.

[0023] The horizontal axis of the autocorrelation analysis graph represents the time interval h, and the vertical axis represents the autocorrelation function value. When a point in the graph is outside the shaded area, it is considered that the time series data has autocorrelation of the time interval h; when all the points in the graph are within the shaded area, it is considered that the time series data has weak autocorrelation data and the autocorrelation needs to be enhanced.

[0024] For weak autocorrelation data, the logarithm method is used to enhance the autocorrelation. The formula is:

[0025] data=lg(time_series) (3)

[0026] Where time_series is the time series of microseismic events.

[0027] If any point in the autocorrelation plot is outside the shaded area, this method can be used and step 3 is performed. If all points in the plot are within the shaded area, the data autocorrelation needs to be enhanced using formula (3). If no point in the enhanced autocorrelation plot is still outside the shaded area, the data is not suitable for prediction using this method and the process ends.

[0028] When analyzing the autocorrelation of microseismic data:

[0029] When calculating the autocorrelation of earthquake occurrence time, first use formula (1) to obtain the earthquake occurrence time interval series data I t , will I t+h As the time interval sequence I t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24;

[0030] 2.2.2) When calculating the autocorrelation of the earthquake source coordinates, the coordinates X, Y, and Z are used as the time series data X t 、Y t 、Z t , X t+h 、Y t+h 、Z t+h As coordinates X, Y, Z time series X t 、Yt , Z t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24;

[0031] 2.2.3) When calculating microseismic energy, the microseismic energy is used as time series data E t , E t+h As the microseismic energy time series E t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis diagram when the value of h is 0 to 24.

[0032] The specific method of step 3) is:

[0033] 3.1) Each column of the microseismic data obtained after preprocessing in step 2) is normalized using the maximum and minimum normalization method. The formula is:

[0034]

[0035] where x max is the maximum value of the sample data, x min is the minimum value of the sample data, and x is the value to be normalized.

[0036] 3.2) For the time interval I, coordinates X, Y, Z, and energy E of the microseismic data, the corresponding relationship between the microseismic event precursor sequence and the future value is constructed to obtain a supervised dataset. The construction rules are as follows:

[0037] I xn ={i n ,i n+1 ,i n+2 ,...,i time_step+n-1} (5)

[0038] I yn =i time_step+n (6)

[0039] Among them I xn Indicates the nth model input data, I yn Indicates the label corresponding to the nth model input data, i n Represents the value of the time interval series data of the microseismic event on the nth day, time_step is a fixed constant;

[0040] X xn ={x n ,x n+1 ,x n+2 ,...,x time_step+n-1} (7)

[0041] Xyn =x time_step+n (8)

[0042] where X xn Represents the nth model input data, X yn Indicates the label corresponding to the nth model input data, x n Represents the value of the time series data of the microseismic event coordinate X on day n. time_step is a fixed constant;

[0043] Y xn ={y n ,y n+1 ,y n+2 ,...,y time_step+n-1} (9)

[0044] Y yn =y time_step+n (10)

[0045] where Y xn Represents the nth model input data, Y yn Indicates the label corresponding to the nth model input data, y n Represents the value of the time series data of the microseismic event coordinate Y on day n. time_step is a fixed constant;

[0046] Z xn ={z n ,z n+1 ,z n+2 ,...,z time_step+n-1} (11)

[0047] Z yn =z time_step+n (12)

[0048] where Z xn Represents the nth model input data, Z yn Indicates the label corresponding to the nth model input data, z n Represents the value of the time series data of the microseismic event coordinate Z on day n, time_step is a fixed constant;

[0049] E xn ={e n ,e n+1 ,e n+2 ,...,e time_step+n-1} (13)

[0050] E yn =e time_step+n (14)

[0051] Among them E xnIndicates the nth model input data, E yn Indicates the label corresponding to the nth model input data, e n Represents the value of the time series data of the microseismic event energy E on the nth day, and time_step is a fixed constant;

[0052] 3.3) The first 70% of the constructed dataset is used as the training dataset, and the last 30% is used as the test dataset.

[0053] The microseismic prediction model in step 4) is specifically:

[0054] Model structure: Butterworth filter BWF is used to decompose the data, and the input is I obtained in step 3 xn 、X xn 、Y xn 、Z xn 、E xn ; Output is I xn 、X xn 、Y xn 、Z xn 、E xn The decomposed low-frequency, medium-frequency and high-frequency data; the bidirectional gated recurrent unit BiGRU is used to extract data features, and there are three BiGRUs of low frequency, medium frequency and high frequency, which are used to extract low-frequency data, medium frequency data and high frequency data features respectively; the low-frequency BiGRU input is low-frequency data, and the output is the low-frequency prediction value; the medium-frequency BiGRU input is medium-frequency data, and the output is the medium-frequency prediction value; the high-frequency BiGRU input is high-frequency data, and the output is the high-frequency prediction value; the three-band prediction values ​​are added together to form the final prediction result, and the final prediction result needs to be combined with the I in step 3 yn 、X yn 、Y yn 、Z yn 、E yn For comparison, the Adam optimization algorithm is used for backpropagation to update the BiGRU weight matrix; Bayesian optimization is used to optimize the hyperparameters of BWF and BiGRU, and the trained model is used for prediction in subsequent steps.

[0055] 4.1) Decompose the data using Butterworth filter (BWF);

[0056] The low-pass, band-pass, and high-pass filters of the Butterworth filter BWF are used to decompose the microseismic event time series data into three different frequency bands: low frequency, medium frequency, and high frequency. The system function of the Butterworth filter is:

[0057]

[0058] Where N is the filter order and Ω is the normalized cutoff frequency.

[0059] The filtering process is the process of solving a linear difference equation with constant coefficients, specifically:

[0060]

[0061]

[0062] Among them, the x(n) sequence is the signal sequence before filtering, a k ,b m is the system array of the denominator and numerator of the system function, and y(n) is the filtered signal sequence.

[0063] 4.2) Use BiGRU to extract the features of the three frequency bands respectively;

[0064] BiGRU is implemented by two layers of GRU structure through propagation in different directions. Finally, the output results of the two layers of GRU hidden layers are merged into one result as the BiGRU hidden layer output. The two layers of GRU are described as forward GRU and reverse GRU respectively.

[0065] Forward GRU calculation process:

[0066] The gated recurrent unit GRU contains two gate structures, namely the update gate z t and reset gate r t The update gate is used to control the extent to which the historical data of the microseismic event at the previous moment is brought into the current state. The larger the value of the update gate, the more historical data of the microseismic event at the previous moment is brought into the current state. The formula is:

[0067]

[0068] The reset gate is used to control the degree of ignoring the historical data information of microseismic events at the previous moment. The smaller the reset gate value, the more historical data information of microseismic events at the previous moment is ignored. The formula is:

[0069]

[0070] in is the input at the current moment, is the hidden state at the previous moment, To update the gate parameter matrix to be trained, is the parameter matrix to be trained for the reset gate, and σ is the sigmoid activation function.

[0071] When calculating the hidden state at the current moment, a candidate state is first calculated And consider resetting the gate value. Candidate state The formula is:

[0072]

[0073] in is the parameter matrix to be trained, is the reset gate output, is the input at the current moment, It is the output at the previous moment.

[0074] Finally, the current output formula of the model is:

[0075]

[0076] in is the update gate output, is the output of the previous moment, is the update gate output, Candidate status.

[0077] Reverse GRU calculation process:

[0078]

[0079]

[0080]

[0081]

[0082] Output the forward GRU hidden layer With the reverse GRU hidden layer output Superposition is performed to obtain the BiGRU hidden layer output h t for:

[0083]

[0084] 4.3) Use Bayesian optimization algorithm to optimize the hyperparameters of BWF and BiGRU;

[0085] The Bayesian optimization algorithm process is:

[0086] 4.3.1) According to the existing parameter adjustment history H = (x 1:t ,f(x 1:k )) , establish a probability distribution model M;

[0087] 4.3.2) Select the next hyperparameter x based on the acquisition function t+1 ;

[0088]

[0089]

[0090] in It represents the best y value in the parameter adjustment history library after t iterations, EI represents the conditional expectation, M t Represents the probability distribution model generated based on the existing t observation value;

[0091] 4.3.3) The new observation (x t+1 ,f(x t+1 )) added to H;

[0092] 4.3.4) Repeat steps 4.3.1) to 4.3.3) until the maximum number of iterations is reached.

[0093] The specific method of step 5) is:

[0094] 5.1) Optimize the hyperparameters of the prediction model in step 4:

[0095] Specify the Bayesian optimization algorithm objective function and continuously reduce the objective function to find the optimal hyperparameter value of the model. The objective function is expressed as:

[0096]

[0097] Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, high, middle, and low represent high-frequency, medium-frequency, and low-frequency prediction models respectively, and n represents the number of samples.

[0098] 5.2) Load the BiGRU weight matrix of the prediction model trained in step 4) with the optimal hyperparameters and construct a mapping relationship between the input microseismic sequence data and the output value;

[0099] The Adam optimization algorithm is used to train the BiGRU weight matrix, and the loss function is:

[0100]

[0101] Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, and n represents the number of samples.

[0102] The specific method of step 6) is:

[0103] Select evaluation indicators suitable for regression prediction problems: absolute percentage error, root mean square error;

[0104] The absolute percentage error (APE) is the difference between the true value and the predicted value of the microseismic event time series data, and then divided by the true value. It is mainly used to measure the relative error between the true value and the predicted value of each sample point in the microseismic event time series data. The smaller the value, the higher the accuracy of the model. The absolute percentage error (APE) calculation formula is:

[0105]

[0106] The root mean square error (RMSE) is mainly used to measure the average deviation between the true value and the predicted value of multiple sample points in the microseismic event time series. The smaller the value, the higher the accuracy of the model. The root mean square error (RMSE) calculation formula is:

[0107]

[0108] The mean absolute percentage error is the sum of the absolute percentage errors of each sample point and the average value. The smaller the value, the higher the accuracy of the model. The calculation formula of MAPE is as follows:

[0109]

[0110] where y i is the true value of the microseismic event time series data, is the predicted value of microseismic event time series data, and N is the number of samples.

[0111] When the MAPE value is less than 10%, the model training is considered complete; when the MAPE value is greater than 10%, it is necessary to return to step 5), increase the number of iterations of the Bayesian optimization algorithm, re-optimize the hyperparameters, and retrain the model.

[0112] The specific method of step 7) is:

[0113] 7.1) Load the optimal hyperparameters for the Butterworth filter and input the microseismic event data into the Butterworth low-pass, band-pass, and high-pass filters, respectively. The Butterworth filter decomposes the raw data into low-frequency, medium-frequency, and high-frequency data.

[0114] 7.2) Normalize the microseismic event data for each frequency band and input them into the trained low-frequency, medium-frequency, and high-frequency bidirectional gated recurrent units (BiGRUs) to obtain the predicted values ​​for each frequency band.

[0115] 7.3) Denormalize the prediction results and add the prediction values ​​of the three frequency bands to obtain the final prediction results of the microseismic event data.

[0116] The beneficial effects of the present invention are:

[0117] 1. A deep learning regression method was used to quantitatively predict the occurrence time, X, Y, and Z coordinates of microseismic events, and their energy. In particular, two high-energy microseismic events were predicted in this embodiment. This invention addresses the current widespread problem of a lack of quantifiable mapping between field monitoring parameters and the time, space, and intensity of disaster events.

[0118] 2. The parameters in the prediction model are obtained using an algorithmic optimization method without manual specification. This quantitative method solves the problem of the difficulty of manual parameter adjustment in conventional prediction models.

[0119] 3. The prediction method can be trained based on massive microseismic data to obtain prediction models with different knowledge for quantitative prediction, which solves the problems of poor generalization and poor dynamic prediction of traditional mathematical and physical models. BRIEF DESCRIPTION OF THE DRAWINGS

[0120] Figure 1 This is the original data form diagram;

[0121] Figure 2 Flowchart of the present invention;

[0122] Figure 3 This is the structure diagram of the BWF-BiGRU-BO model;

[0123] Figure 4 It is a quantitative prediction flow chart;

[0124] Figure 5 This is the data autocorrelation analysis diagram;

[0125] Figure 6 Comparison between the actual and predicted values ​​of the time interval between two microseismic events in a certain mine working face;

[0126] Figure 7 Comparison between the actual value and predicted value of the microseismic source coordinate X of a certain mine working face;

[0127] Figure 8 Comparison between the true value and predicted value of the microseismic source coordinate Y of a certain mine working face;

[0128] Figure 9 Comparison between the actual value and predicted value of the microseismic source coordinate Z of a certain mine working face;

[0129] Figure 10 Comparison between the actual and predicted values ​​of microseismic energy released at a certain mine working face;

[0130] Figure 11 Hyperparameter optimization results for the BWF-BiGRU-BO model. DETAILED DESCRIPTION

[0131] The present invention provides a quantitative prediction method for the temporal and spatial intensity three elements of rock burst mine microseismic events, which specifically includes the following steps:

[0132] 1) Import historical microseismic data, including the time of occurrence (year, month, day, hour, minute, second), focal coordinates, and microseismic energy;

[0133] When a microseismic event occurs, the rock mass breaks, energy is released, and vibration waves are generated. After multiple sensors collect the vibration waves, the program calculates three elements: the time of the earthquake, the coordinates of the earthquake source, and the time series data of the microseismic energy. The data format is as follows: Figure 1 The earthquake time is composed of year / month / day hour:minute:second, and each time corresponds to a microseismic event. The source coordinates are composed of three components: X, Y, and Z, which are the relative offset of the coordinate origin defined by the mining area, in meters. The microseismic energy represents the amount of energy released by the microseismic event, in joules.

[0134] 2) Preprocessing of microseismic data;

[0135] 2.1) Using the difference algorithm to process the earthquake time series data into time interval series data I;

[0136] I=t next -t now (1)

[0137] where t next is the time when the next microseismic event occurs, t now is the time when the current microearthquake occurs.

[0138] 2.2) Analyze the autocorrelation of microseismic data. Calculate the autocorrelation of the earthquake occurrence time (year-month-day hour: minute: second), the focal coordinates, and the microseismic energy. Specifically,

[0139] Input the time series data of the time interval obtained from Formula 1 into the autocorrelation function for calculation to obtain the autocorrelation analysis graph of the time interval interval. Input the time series data of coordinates X, Y, Z, and energy into the autocorrelation function to obtain the autocorrelation analysis graphs corresponding to coordinates X, Y, Z, and energy, respectively.

[0140] The autocorrelation function is independent of any specific time t and is a function of the time interval h, expressed as ρ h express:

[0141]

[0142] Where h is the time interval, cov(X t ,X t+h ) represents the variable sequence X t ,X t+h The covariance of t σt+h Represents standard deviation.

[0143] The form of the autocorrelation analysis diagram is as follows Figure 3 As shown in the figure, the horizontal axis represents the time interval h, and the vertical axis represents the autocorrelation function value. When a point in the figure is outside the shaded area, the time series data is considered to have autocorrelation for the time interval h. When all the points in the figure are within the shaded area, the time series data is considered to have weak autocorrelation and needs to be enhanced. The autocorrelation value for the time interval h of 0 is always outside the shaded area and is not discussed here.

[0144] For weak autocorrelation data, the logarithm method is used to enhance the autocorrelation. The formula is:

[0145] data=lg(time_series) (3)

[0146] Where time_series is the time series of microseismic events.

[0147] When calculating the autocorrelation of earthquake occurrence time, first use formula (1) to obtain the earthquake occurrence time interval series data I t , will I t+h As the time interval sequence I t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24;

[0148] 2.2.2) When calculating the autocorrelation of the earthquake source coordinates, the coordinates X, Y, and Z are used as the time series data X t 、Y t 、Z t , X t+h 、Y t+h 、Z t+h As coordinates X, Y, Z time series X t 、Y t 、Z t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24;

[0149] 2.2.3) When calculating microseismic energy, the microseismic energy is used as time series data E t , E t+h As the microseismic energy time series E t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis diagram when the value of h is 0 to 24.

[0150] If any point in the autocorrelation plot is outside the shaded area, this method can be used and step 3 is performed. If all points in the plot are within the shaded area, the data autocorrelation needs to be enhanced using formula (3). If no point in the enhanced autocorrelation plot is still outside the shaded area, the data is not suitable for prediction using this method and the process ends.

[0151] 3) Using the processed microseismic data to create microseismic training datasets and microseismic test datasets;

[0152] 3.1) Each column of the microseismic data obtained after preprocessing in step 2 is normalized using the maximum and minimum normalization method. The formula is:

[0153]

[0154] where x max is the maximum value of the sample data, x min is the minimum value of the sample data, and x is the value to be normalized.

[0155] 3.2) For the time interval I, coordinates X, Y, Z, and energy E of the microseismic data, the corresponding relationship between the microseismic event precursor sequence and the future value is constructed to obtain a supervised dataset. The construction rules are as follows:

[0156] I xn ={i n ,i n+1 ,i n+2 ,...,i time_step+n-1} (5)

[0157] I yn =i time_step+n (6)

[0158] Among them I xn Indicates the nth model input data, I yn Indicates the label corresponding to the nth model input data, i n Represents the value of the time interval series data of the microseismic event on the nth day, time_step is a fixed constant;

[0159] X xn ={x n ,x n+1 ,x n+2 ,...,x time_step+n-1} (7)

[0160] X yn =x time_step+n (8)

[0161] where X xn Represents the nth model input data, Xyn Indicates the label corresponding to the nth model input data, x n Represents the value of the time series data of the microseismic event coordinate X on day n. time_step is a fixed constant;

[0162] Y xn ={y n ,y n+1 ,y n+2 ,...,y time_step+n-1} (9)

[0163] Y yn =y time_step+n (10)

[0164] where Y xn Represents the nth model input data, Y yn Indicates the label corresponding to the nth model input data, y n Represents the value of the time series data of the microseismic event coordinate Y on day n. time_step is a fixed constant;

[0165] Z xn ={z n ,z n+1 ,z n+2 ,...,z time_step+n-1} (11)

[0166] Z yn =z time_step+n (12)

[0167] where Z xn Represents the nth model input data, Z yn Indicates the label corresponding to the nth model input data, z n Represents the value of the time series data of the microseismic event coordinate Z on day n, time_step is a fixed constant;

[0168] E xn ={e n ,e n+1 ,e n+2 ,...,e time_step+n-1} (13)

[0169] E yn =e time_step+n (14)

[0170] Among them E xn Indicates the nth model input data, E yn Indicates the label corresponding to the nth model input data, e nRepresents the value of the time series data of the microseismic event energy E on the nth day, and time_step is a fixed constant;

[0171] 3.3) The first 70% of the constructed dataset is used as the training dataset, and the last 30% is used as the test dataset;

[0172] 4) Construct a BWF-BiGRU-BO microseismic prediction model based on Butterworth filter (BWF), bidirectional gated recurrent unit (BiGRU), and Bayesian optimization (BO);

[0173] The model structure is as follows Figure 3 As shown: Butterworth filter BWF is used to decompose the data, and the input is I obtained in step 3 xn 、X xn 、Y xn , Z xn 、E xn ; Output is I xn 、X xn 、Y xn , Z xn 、E xn The decomposed low-frequency, medium-frequency and high-frequency data; the bidirectional gated recurrent unit BiGRU is used to extract data features, and there are three BiGRUs of low frequency, medium frequency and high frequency, which are used to extract low-frequency data, medium frequency data and high frequency data features respectively; the low-frequency BiGRU input is low-frequency data, and the output is the low-frequency prediction value; the medium-frequency BiGRU input is medium-frequency data, and the output is the medium-frequency prediction value; the high-frequency BiGRU input is high-frequency data, and the output is the high-frequency prediction value; the three-band prediction values ​​are added together to form the final prediction result, and the final prediction result needs to be combined with the I in step 3 yn 、X yn 、Y yn , Z yn 、E yn For comparison, the Adam optimization algorithm was used to backpropagate and update the BiGRU weight matrix. Bayesian optimization was used to optimize the hyperparameters of the BWF and BiGRU, and the trained models were used for prediction in subsequent steps.

[0174] 4.1) Decompose the data using Butterworth filter (BWF);

[0175] The low-pass, band-pass, and high-pass filters of the Butterworth filter (BWF) are used to decompose the microseismic event time series data into three different frequency bands: low frequency, medium frequency, and high frequency. The system function of the Butterworth filter is:

[0176]

[0177] Where N is the filter order and Ω is the normalized cutoff frequency.

[0178] The filtering process is the process of solving a linear difference equation with constant coefficients, specifically:

[0179]

[0180]

[0181] Among them, the x(n) sequence is the signal sequence before filtering, a k ,b m is the system array of the denominator and numerator of the system function, and y(n) is the filtered signal sequence.

[0182] 4.2) Use BiGRU to extract the features of the three frequency bands respectively;

[0183] BiGRU is implemented by two layers of GRU structure through propagation in different directions. Finally, the output results of the two layers of GRU hidden layers are merged into one result as the BiGRU hidden layer output. The two layers of GRU are described as forward GRU and reverse GRU respectively.

[0184] Forward GRU calculation process:

[0185] The gated recurrent unit GRU contains two gate structures, namely the update gate z t and reset gate r t The update gate is used to control the extent to which the historical data of the microseismic event at the previous moment is brought into the current state. The larger the value of the update gate, the more historical data of the microseismic event at the previous moment is brought into the current state. The formula is:

[0186]

[0187] The reset gate is used to control the degree of ignoring the historical data information of microseismic events at the previous moment. The smaller the reset gate value, the more historical data information of microseismic events at the previous moment is ignored. The formula is:

[0188]

[0189] in is the input at the current moment, is the hidden state at the previous moment, To update the gate parameter matrix to be trained, is the parameter matrix to be trained for the reset gate, and σ is the sigmoid activation function.

[0190] When calculating the hidden state at the current moment, a candidate state is first calculated And consider resetting the gate value. Candidate state The formula is:

[0191]

[0192] in is the parameter matrix to be trained, is the reset gate output, is the input at the current moment, It is the output at the previous moment.

[0193] Finally, the current output formula of the model is:

[0194]

[0195] in is the update gate output, is the output of the previous moment, is the update gate output, Candidate status.

[0196] Reverse GRU calculation process:

[0197]

[0198]

[0199]

[0200]

[0201] Output the forward GRU hidden layer With the reverse GRU hidden layer output Superposition is performed to obtain the BiGRU hidden layer output h t for:

[0202]

[0203] 4.3) Use Bayesian optimization algorithm to optimize the hyperparameters of BWF and BiGRU;

[0204] The hyperparameters to be optimized for the BWF-BiGRU-BO model are as follows: Figure 3 As shown. The Bayesian optimization algorithm process is:

[0205] 1. According to the existing parameter adjustment history H = (x 1:t ,f(x 1:k )) , establish a probability distribution model M;

[0206] 2. Select the next hyperparameter x based on the acquisition function t+1 ;

[0207] 3. Translate the new observation (x t+1,f(x t+1 )) added to H;

[0208] 4. Repeat steps 1 to 3 until the maximum number of iterations is reached.

[0209] Step 2 in the Bayesian optimization algorithm is specifically expressed as:

[0210]

[0211]

[0212] in It represents the best y value in the parameter adjustment history library after t iterations, EI represents the conditional expectation, M t Represents the probability distribution model generated based on the existing t observations.

[0213] 5) Input the microseismic training set data into the prediction model in step 4 and train it

[0214] 5.1) Optimize the hyperparameters of the prediction model in step 4;

[0215] Specify the Bayesian optimization algorithm objective function and continuously reduce the objective function to find the optimal hyperparameter value of the model. The objective function is expressed as:

[0216]

[0217] Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, high, middle, and low represent high-frequency, medium-frequency, and low-frequency prediction models respectively, and n represents the number of samples.

[0218] 5.2) Load the BiGRU weight matrix of the prediction model in step 4 of the optimal hyperparameter training to construct the mapping relationship between the input microseismic sequence data and the output value;

[0219] The Adam optimization algorithm is used to train the BiGRU weight matrix, and the loss function is:

[0220]

[0221] Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, and n represents the number of samples.

[0222] 6) Input the microseismic test data into the trained model and use the evaluation indicators to verify the model effect in step 5;

[0223] The evaluation indicators suitable for regression prediction problems are selected: Absolute Percentage Error (APE), Root Mean Square Error (RMSE) and Mean Absolute Percentage Error (MAPE).

[0224] Absolute percentage error (APE) is the difference between the true value and the predicted value of the microseismic event time series data, and then divided by the true value. It is mainly used to measure the relative error between the true value and the predicted value of each sample point in the microseismic event time series data. The smaller the value, the higher the accuracy of the model. The absolute percentage error (APE) calculation formula is:

[0225]

[0226] The root mean square error (RMSE) is mainly used to measure the average deviation between the true value and the predicted value of multiple sample points in the microseismic event time series. The smaller the value, the higher the accuracy of the model. The formula for calculating the root mean square error (RMSE) is:

[0227]

[0228] The mean absolute percentage error is the sum of the absolute percentage errors of each sample point and the average value. The smaller the value, the higher the accuracy of the model. The calculation formula of MAPE is as follows:

[0229]

[0230] where y i is the true value of the microseismic event time series data, is the predicted value of microseismic event time series data, and N is the number of samples.

[0231] When the MAPE value is less than 10%, the model training is considered complete; when the MAPE value is greater than 10%, it is necessary to return to step 5), increase the number of iterations of the Bayesian optimization algorithm, re-optimize the hyperparameters, and retrain the model.

[0232] 7) The actual microseismic sequence data obtained: earthquake time interval, location coordinates X, Y, Z, and released energy are input into the prediction model trained in step 5 for quantitative prediction.

[0233] 7.1) Load the optimal hyperparameters for the Butterworth filter and input the microseismic event data into the Butterworth low-pass, band-pass, and high-pass filters, respectively. The Butterworth filter decomposes the raw data into low-frequency, medium-frequency, and high-frequency data.

[0234] 7.2) Normalize the microseismic event data for each frequency band and input them into the trained low-frequency, medium-frequency, and high-frequency bidirectional gated recurrent units (BiGRUs) to obtain the predicted values ​​for each frequency band.

[0235] 7.3) Denormalize the prediction results and add the prediction values ​​of the three frequency bands to obtain the final prediction results of the microseismic event data.

[0236] Example 1:

[0237] A coal mine operating at a depth of nearly 1000 meters experienced numerous high-energy microseismic events during underground mining, with noticeable tremors felt at the surface. We collected historical microseismic data from the mine's 105 working face, monitored by the mine's microseismic monitoring system from October 19, 2020, to June 27, 2021.

[0238] The present invention provides a method for quantitatively predicting long-term and short-term microseismic events in rock burst mines, which is characterized by comprising the following steps:

[0239] Step S1: Import historical microseismic data, including the time of occurrence (year, month, day, hour, minute, second), focal coordinates, and microseismic energy;

[0240] When a microseismic event occurs, the rock mass breaks, energy is released, and vibration waves are generated. After multiple sensors collect the vibration waves, the program calculates three elements: the time of the earthquake, the coordinates of the earthquake source, and the time series data of the microseismic energy. The data format is as follows: Figure 1 The earthquake time is composed of year / month / day hour:minute:second, and each time corresponds to a microseismic event. The source coordinates are composed of three components: X, Y, and Z, which are the relative offset of the coordinate origin defined by the mining area, in meters. The microseismic energy represents the amount of energy released by the microseismic event, in joules.

[0241] Step S2: preprocessing the microseismic data;

[0242] 2.1) Use formula (1) to process the earthquake time series data into time interval series data I;

[0243] 2.2) Analyze the autocorrelation of microseismic data. Calculate the autocorrelation of the earthquake occurrence time (year-month-day hour: minute: second), the focal coordinates, and the microseismic energy. Specifically,

[0244] When calculating the autocorrelation of earthquake occurrence time, first use formula (1) to obtain the earthquake occurrence time interval series data I t , will I t+h As the time interval sequence I t The new sequence is shifted by h units on the time axis; the autocorrelation analysis when the value of h is 0 to 24 is obtained in turn using formula (2) Figure 5 ;

[0245] 2.2.2) When calculating the autocorrelation of the earthquake source coordinates, the coordinates X, Y, and Z are used as the time series data X t 、Y t 、Z t , X t+h 、Y t+h 、Z t+h As coordinates X, Y, Z time series X t 、Y t 、Z t The new sequence is shifted by h units on the time axis; the autocorrelation analysis when the value of h is 0 to 24 is obtained in turn using formula (2) Figure 5 ;

[0246] 2.2.3) When calculating microseismic energy, the microseismic energy is used as time series data E t , E t+h As the microseismic energy time series E t The new sequence is shifted by h units on the time axis; the autocorrelation analysis when the value of h is 0 to 24 is obtained in turn using formula (2) Figure 5 .

[0247] Depend on Figure 5 It can be found that the autocorrelation value points of the time series of time interval, coordinate X, coordinate Y, and coordinate Z are not all located within the shaded area, so they have good autocorrelation and can proceed to step 3. The autocorrelation value points of the energy release time series data are all located within the shaded area except for the autocorrelation points with the time interval of 0. Therefore, formula (3) is used to perform data enhancement on this series data.

[0248] Step S3: using the processed microseismic data to create a microseismic training dataset and a microseismic test dataset;

[0249] 3.1) Each column of the microseismic data obtained after preprocessing in step 2 is normalized using formula (4);

[0250] 3.2) The corresponding relationship between the microseismic event precursor sequence and the future value is constructed for the time interval I, coordinate X, coordinate Y, coordinate Z and energy E of the microseismic data to obtain a supervised data set. The construction rules are formulas (5)(6)(7)(8)(9)(10)(11)(12)(13)(14);

[0251] 3.3) The first 70% of the constructed dataset is used as the training dataset, and the last 30% is used as the test dataset;

[0252] Step S4: constructing a BWF-BiGRU-BO microseismic prediction model based on Butterworth filter BWF, bidirectional gated recurrent unit BiGRU and Bayesian optimization BO;

[0253] Butterworth filter BWF is used to decompose the data, and the input is I obtained in step 3 xn 、X xn 、Y xn 、Z xn 、E xn ; Output is I xn 、X xn 、Y xn 、Z xn 、E xn The decomposed low-frequency, medium-frequency and high-frequency data; the bidirectional gated recurrent unit BiGRU is used to extract data features, and there are three BiGRUs of low frequency, medium frequency and high frequency, which are used to extract low-frequency data, medium frequency data and high frequency data features respectively; the low-frequency BiGRU input is low-frequency data, and the output is the low-frequency prediction value; the medium-frequency BiGRU input is medium-frequency data, and the output is the medium-frequency prediction value; the high-frequency BiGRU input is high-frequency data, and the output is the high-frequency prediction value; the three-band prediction values ​​are added together to form the final prediction result, and the final prediction result needs to be combined with the I in step 3 yn 、X yn 、Y yn 、Z yn 、E yn For comparison, the Adam optimization algorithm was used to backpropagate and update the BiGRU weight matrix. Bayesian optimization was used to optimize the hyperparameters of the BWF and BiGRU, and the trained models were used for prediction in subsequent steps.

[0254] Step S5: input the microseismic training set data into the prediction model in step 4 to train it;

[0255] 5.1) Optimize the hyperparameters of the prediction model in step 4;

[0256] Formula (29) is used as the objective function of the Bayesian optimization algorithm, and the objective function is continuously reduced to find the optimal hyperparameter value of the model.

[0257] 5.2) Load the optimal hyperparameters and use the Adam optimization algorithm to train the BiGRU weight matrix of the prediction model in step 4 to construct the mapping relationship between the input microseismic sequence data and the output value. The loss function is shown in Equation (30);

[0258] Step S6: Input the microseismic test set data into the trained model and use the evaluation indicators (31), (32), and (33) to verify the model effect;

[0259] When the MAPE value is less than 10%, the model training is considered complete; when the MAPE value is greater than 10%, it is necessary to return to step 5), increase the number of iterations of the Bayesian optimization algorithm, re-optimize the hyperparameters, and retrain the model.

[0260] The training set data of time interval, coordinate X, coordinate Y, coordinate Z and released energy are input into the corresponding and trained prediction model respectively, and the evaluation index results are as follows:

[0261] The RMSE of the predicted time interval is 0.260785 and the MAPE is 3.168639%;

[0262] The RMSE of the predicted coordinate X is 24.730886 and the MAPE is 0.0045%;

[0263] The RMSE of the predicted coordinate Y is 39.498736 and the MAPE is 0.0005%;

[0264] The RMSE of the predicted coordinate Z is 9.899231 and the MAPE is 7.7258%;

[0265] The RMSE of the predicted released energy is 0.086791 and the MAPE is 3.9277%.

[0266] Since the MAPE values ​​are all less than 10%, the trained models meet the prediction requirements.

[0267] Step S7: The actual microseismic sequence data, including the time interval between earthquakes, the position coordinates X, Y, Z, and the released energy, are input into the prediction model trained in step 5 for quantitative prediction.

[0268] 7.1) Load the optimal hyperparameters for the Butterworth filter and input the microseismic event data into the Butterworth low-pass, band-pass, and high-pass filters, respectively. The Butterworth filter decomposes the raw data into low-frequency, medium-frequency, and high-frequency data.

[0269] 7.2) Normalize the microseismic event data for each frequency band and input them into the trained low-frequency, medium-frequency, and high-frequency bidirectional gated recurrent units (BiGRUs) to obtain the predicted values ​​for each frequency band.

[0270] 7.3) Denormalize the prediction results and add the prediction values ​​of the three frequency bands to obtain the final prediction results of the microseismic event data.

[0271] Comparison between the actual and predicted time intervals between two microseismic events at a certain mine working face Figure 6 As shown;

[0272] Comparison between the actual value and predicted value of the microseismic source coordinate X of a mine working face Figure 7 As shown;

[0273] Comparison between the actual value and predicted value of the microseismic source coordinate Y of a mine working face Figure 8 As shown;

[0274] Comparison between the actual value and predicted value of the microseismic source coordinate Z of a mine working face Figure 9 As shown;

[0275] Comparison between the actual value and predicted value of microseismic energy released by a mine working face Figure 10 shown.

[0276] The hyperparameter optimization results of the BWF-BiGRU-BO model are as follows Figure 11 shown.

Claims

1. A quantitative prediction method for the three spatial and temporal elements of rock burst mine microseismic events, characterized by: The steps are: 1) Import historical microseismic data, including earthquake time, source coordinates, and microseismic energy; 2) Preprocessing of microseismic data; 2.1) Using the difference algorithm to process the earthquake time series data into time interval series data I; I=t next -t now (1) where t next is the time when the next microseismic event occurs, t now is the time when the current microseismic event occurs; 2.2) Analyze the autocorrelation of microseismic data. Calculate the autocorrelation of the earthquake occurrence time, source coordinates, and microseismic energy data. Specifically: Input the time series data of the time interval obtained by formula 1 into the autocorrelation function for calculation to obtain the autocorrelation analysis graph of the time interval I; input the time series data of coordinates X, coordinates Y, coordinates Z and energy E into the autocorrelation function to obtain the autocorrelation analysis graphs corresponding to coordinates X, coordinates Y, coordinates Z and energy E respectively; The autocorrelation function is independent of any specific time t and is a function of the time interval h, expressed as ρ h express: Among them, h is the time interval, cov(X t ,X t+h ) represents the variable sequence X t ,X t+h The covariance of t σ t+h represents the standard deviation; The horizontal axis of the autocorrelation analysis graph represents the time interval h, and the vertical axis represents the autocorrelation function value. When a point in the graph is outside the shaded area, it is considered that the time series data has autocorrelation of the time interval h; when all the points in the graph are within the shaded area, it is considered that the time series data has weak autocorrelation data and the autocorrelation needs to be enhanced. For weak autocorrelation data, the logarithm method is used to enhance the autocorrelation. The formula is: data=lg(time_series) (3) Where time_series is the time series of microseismic events; When there is a point outside the shaded area in the autocorrelation graph, it is considered that this method can be used and step 3 is performed; when all the points in the graph are within the shaded area, it is necessary to enhance the data autocorrelation using formula (3); if there is still no point outside the shaded area in the autocorrelation graph of the enhanced data, then the data is not suitable for prediction using this method and the process ends; 3) Using the processed microseismic data to create microseismic training datasets and microseismic test datasets; 4) Construct a BWF-BiGRU-BO microseismic prediction model based on Butterworth filter (BWF), bidirectional gated recurrent unit (BiGRU), and Bayesian optimization (BO); 5) Input the microseismic training set data into the prediction model in step 4 to train it; 6) Input the microseismic test data set into the model trained in step 5) and use the evaluation index to verify the effect of the model in step 5. If the effect meets the expectations, proceed to step 7). Otherwise, increase the number of iterations and retrain the model through step 5). 7) Input the actual microseismic sequence data into the prediction model trained in step 5 for quantitative prediction.

2. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: In the step 1), the historical microseismic data are specifically as follows: when a microseismic event occurs, the rock mass fractures, energy is released, and vibration waves are generated. After multiple sensors collect the vibration waves, the program calculates three elements based on the arrival time information: the time of earthquake occurrence, the source coordinates, and the time series data of the microseismic energy; the time of earthquake occurrence is composed of year / month / day hour: minute: second, and each moment corresponds to a microseismic event; the source coordinates are composed of three components of X, Y, and Z, which are the relative offset of the coordinate origin defined by the mining area, and the unit is meter; the microseismic energy represents the amount of energy released by this microseismic event, and the unit is joule.

3. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: When analyzing the autocorrelation of microseismic data: 2.2.1) When calculating the autocorrelation of earthquake occurrence time, first use formula (1) to obtain the earthquake occurrence time interval series data I t , will I t+h As the time interval sequence I t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24; 2.2.2) When calculating the autocorrelation of the earthquake source coordinates, the coordinates X, Y, and Z are used as the time series data X t 、Y t , Z t , X t+h 、Y t+h , Z t+h As coordinates X, Y, Z time series X t 、Y t , Z t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis graph when the value of h is 0 to 24; 2.2.3) When calculating microseismic energy, the microseismic energy is used as time series data E t , E t+h As the microseismic energy time series E t Shift the new sequence by h units on the time axis; use formula (2) to obtain the autocorrelation analysis diagram when the value of h is 0 to 24.

4. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: The specific method of step 3) is: 3.1) Each column of the microseismic data obtained after preprocessing in step 2) is normalized using the maximum and minimum normalization method. The formula is: where x max is the maximum value of the sample data, x min is the minimum value of the sample data, and x is the value to be normalized; 3.2) For the time interval I, coordinates X, Y, Z, and energy E of the microseismic data, the corresponding relationship between the microseismic event precursor sequence and the future value is constructed to obtain a supervised dataset. The construction rules are as follows: I xn ={i n ,i n+1 ,i n+2 ,...,i time_step+n-1 } (5) I yn =i time_step+n (6) Among them I xn Indicates the nth model input data, I yn Indicates the label corresponding to the nth model input data, i n Represents the value of the time interval series data of the microseismic event on the nth day, time_step is a fixed constant; X xn ={x n ,x n+1 ,x n+2 ,...,x time_step+n-1 } (7) X yn =x time_step+n (8) where X xn Represents the nth model input data, X yn Indicates the label corresponding to the nth model input data, x n Represents the value of the time series data of the microseismic event coordinate X on day n. time_step is a fixed constant; AND xn ={and n ,and n+1 ,and n+2 ,...,and time_step+n-1 } (9) AND yn =and time_step+n (10) where Y xn Represents the nth model input data, Y yn Indicates the label corresponding to the nth model input data, y n Represents the value of the time series data of the microseismic event coordinate Y on day n. time_step is a fixed constant; WITH xn ={z n ,With n+1 ,With n+2 ,...,With time_step+n-1 } (11) WITH yn =z time_step+n (12) where Z xn Represents the nth model input data, Z yn Indicates the label corresponding to the nth model input data, z n Represents the value of the time series data of the microseismic event coordinate Z on day n, time_step is a fixed constant; AND xn ={and n ,And n+1 ,And n+2 ,...,And time_step+n-1 } (13) AND yn =and time_step+n (14) Among them E xn Indicates the nth model input data, E yn Indicates the label corresponding to the nth model input data, e n Represents the value of the time series data of the microseismic event energy E on the nth day, and time_step is a fixed constant; 3.3) The first 70% of the constructed dataset is used as the training dataset, and the last 30% is used as the test dataset.

5. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: The microseismic prediction model in step 4) is specifically: Model structure: Butterworth filter BWF is used to decompose the data, and the input is I obtained in step 3 xn 、X xn 、Y xn , Z xn 、E xn ; Output is I xn 、X xn 、Y xn , Z xn 、E xn Decomposed low-frequency, medium-frequency and high-frequency data; Bidirectional gated recurrent unit BiGRU is used to extract data features. There are three BiGRUs for low frequency, medium frequency and high frequency, which are used to extract low-frequency data, medium-frequency data and high-frequency data features respectively; The input of low-frequency BiGRU is low-frequency data, and the output is low-frequency predicted value; The input of medium-frequency BiGRU is medium-frequency data, and the output is medium-frequency prediction value; the input of high-frequency BiGRU is high-frequency data, and the output is high-frequency prediction value; the three-band prediction values ​​are added together to form the final prediction result, which needs to be combined with the I in step 3. yn 、X yn 、Y yn , Z yn 、E yn For comparison, the Adam optimization algorithm is used to back-propagate and update the BiGRU weight matrix. Bayesian optimization is used to optimize the hyperparameters of BWF and BiGRU, and the trained model is used for prediction in subsequent steps. 4.1) Decompose the data using Butterworth filter (BWF); The low-pass, band-pass, and high-pass filters of the Butterworth filter BWF are used to decompose the microseismic event time series data into three different frequency bands: low frequency, medium frequency, and high frequency. The system function of the Butterworth filter is: Where N is the filter order and Ω is the normalized cutoff frequency; The filtering process is the process of solving a linear difference equation with constant coefficients, specifically: Among them, the x(n) sequence is the signal sequence before filtering, a k ,b m is the system array of the denominator and numerator of the system function, y(n) is the filtered signal sequence; 4.2) Use BiGRU to extract the features of the three frequency bands respectively; BiGRU is implemented by two layers of GRU structure through propagation in different directions. Finally, the output results of the two layers of GRU hidden layers are merged into one result as the BiGRU hidden layer output. The two layers of GRU are described as forward GRU and reverse GRU respectively. Forward GRU calculation process: The gated recurrent unit GRU contains two gate structures, namely the update gate z t and reset gate r t The update gate is used to control the extent to which the historical data of the microseismic event at the previous moment is brought into the current state. The larger the value of the update gate, the more historical data of the microseismic event at the previous moment is brought into the current state. The formula is: The reset gate is used to control the degree of ignoring the historical data information of microseismic events at the previous moment. The smaller the reset gate value is, the more historical data information of microseismic events at the previous moment is ignored. The formula is: in is the input at the current moment, is the hidden state at the previous moment, To update the gate parameter matrix to be trained, is the parameter matrix to be trained for the reset gate, σ is the sigmoid activation function; When calculating the hidden state at the current moment, a candidate state is first calculated And consider resetting the gate value; candidate state The formula is: in is the parameter matrix to be trained, is the reset gate output, is the input at the current moment, is the output at the previous moment; Finally, the current output formula of the model is: in is the update gate output, is the output of the previous moment, is the update gate output, is a candidate status; Reverse GRU calculation process: Output the forward GRU hidden layer With the reverse GRU hidden layer output Superposition is performed to obtain the BiGRU hidden layer output h t for: 4.3) Use Bayesian optimization algorithm to optimize the hyperparameters of BWF and BiGRU; The Bayesian optimization algorithm process is: 4.3.1) According to the existing parameter adjustment history H = (x 1:t ,f(x 1:k )) , establish a probability distribution model M; 4.3.2) Select the next hyperparameter x based on the acquisition function t+1 ; in It represents the best y value in the parameter adjustment history library after t iterations, EI represents the conditional expectation, M t Represents the probability distribution model generated based on the existing t observation value; 4.3.3) The new observation (x t+1 ,f(x t+1 )) added to H; 4.3.4) Repeat steps 4.3.1) to 4.3.3) until the maximum number of iterations is reached.

6. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: The specific method of step 5) is: 5.1) Optimize the hyperparameters of the prediction model in step 4: Specify the Bayesian optimization algorithm objective function and continuously reduce the objective function to find the optimal hyperparameter value of the model. The objective function is expressed as: Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, high, middle, and low represent high-frequency, medium-frequency, and low-frequency prediction models respectively, and n represents the number of samples; 5.2) Load the BiGRU weight matrix of the prediction model trained in step 4) with the optimal hyperparameters and construct a mapping relationship between the input microseismic sequence data and the output value; The Adam optimization algorithm is used to train the BiGRU weight matrix, and the loss function is: Where y(x) represents the true value, a L (x) represents the predicted value of neuron output, and n represents the number of samples.

7. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: The specific method of step 6) is: Select evaluation indicators suitable for regression prediction problems: absolute percentage error, root mean square error; The absolute percentage error (APE) is the difference between the true value and the predicted value of the microseismic event time series data, and then divided by the true value. It is mainly used to measure the relative error between the true value and the predicted value of each sample point in the microseismic event time series data. The smaller the value, the higher the accuracy of the model. The absolute percentage error (APE) calculation formula is: The root mean square error (RMSE) is mainly used to measure the average deviation between the true value and the predicted value of multiple sample points in the microseismic event time series. The smaller the value, the higher the accuracy of the model. The root mean square error (RMSE) calculation formula is: The mean absolute percentage error is the sum of the absolute percentage errors of each sample point and the average value. The smaller the value, the higher the accuracy of the model. The calculation formula of MAPE is as follows: where y i is the true value of the microseismic event time series data, is the predicted value of microseismic event time series data, and N is the number of samples; When the MAPE value is less than 10%, the model training is considered complete; when the MAPE value is greater than 10%, it is necessary to return to step 5), increase the number of iterations of the Bayesian optimization algorithm, re-optimize the hyperparameters, and retrain the model.

8. The method for quantitatively predicting the temporal and spatial intensity three elements of rock burst mine microseismic events according to claim 1 is characterized in that: The specific method of step 7) is: 7.1) Load the optimal hyperparameters of the Butterworth filter and input the microseismic event data into the Butterworth low-pass, band-pass, and high-pass filters respectively; The Butterworth filter decomposes the original data as a signal into low-frequency, medium-frequency, and high-frequency data; 7.2) Normalize the microseismic event data for each frequency band and input them into the trained low-frequency, medium-frequency, and high-frequency bidirectional gated recurrent units (BiGRUs) to obtain the predicted values ​​for each frequency band; 7.3) Denormalize the prediction results and add the prediction values ​​of the three frequency bands to obtain the final prediction results of the microseismic event data.

Images