Rock mass abnormal vibration event monitoring and rock mass instability risk assessment method

Abstract

The invention discloses a rock mass abnormal vibration event monitoring and rock mass instability risk assessment method, which comprises the steps that a detection model for detecting rock mass vibration signals is constructed, the detection model comprises an encoder and a decoder connected with the encoder, the encoder and the decoder are respectively formed by connecting a plurality of layers of GRU networks in sequence, constructing a training set to train the detection model, and constructing a dynamic control chart based on exponentially weighted moving average; and abnormal vibration event judgment based on the gating line on the dynamic control chart. According to the method, subtle and early time sequence mode abnormity in a microseismic signal sequence is captured through the detection model, an early warning signal can be sent out before a macroscopic statistical index is obviously changed, the early warning advancement is high, the time window from early warning to catastrophe occurrence is obviously prolonged, and precious time is won for disaster prevention and control.

Description

Technical Field

[0001] This invention relates to the fields of mine safety, underground engineering and geological disaster monitoring and early warning technology, and in particular to a method for anomaly detection and rock mass instability risk assessment based on the temporal characteristics of microseismic signals. Background Technology

[0002] In geotechnical engineering fields such as mining, water conservancy tunnels, and slope engineering, the stability of rock masses is directly related to production safety and personnel lives. Microseismic monitoring technology, as an effective regional and dynamic monitoring method, has been widely used in rock mass stability assessment. This technology uses a sensor network deployed within the rock mass to capture elastic wave signals released during rock fracturing, thereby locating microseismic events and analyzing their activity. However, existing methods for assessing rock mass instability risk based on microseismic monitoring mainly suffer from the following problems and drawbacks:

[0003] (1) Reliance on empirical statistical indicators, resulting in strong early warning lag: Most current methods rely on empirical threshold judgments of macroscopic statistical parameters of microseismic events (such as the number of events, cumulative energy, apparent volume, b-value, etc.). These statistical indicators are macroscopic summaries of events over a period of time, making it difficult to capture the subtle, early temporal dynamic characteristics of the microfracture evolution process before rock mass instability. This leads to early warning signals often appearing near the instability critical point, resulting in a short early warning window and insufficient time for decision-making and response in disaster prevention and control.

[0004] (2) Single feature extraction dimension, ignoring the inherent temporal correlation of signals: Existing technologies usually treat each microseismic event as an independent individual, or only focus on its spatial clustering, while ignoring the dynamic evolution and inherent correlation of the microseismic event sequence on the time axis. The evolution of rock mass from stability to instability is a nonlinear and non-stationary process, and its microseismic activity will show specific patterns in time (such as cluster bursts, quiet periods, frequency component migration, etc.). These key precursor information have not been fully explored in traditional methods.

[0005] (3) Insufficient automation and intelligence in risk assessment: The existing risk assessment process relies heavily on expert experience for manual interpretation and lacks an end-to-end, automated intelligent assessment model. This makes the assessment results highly subjective and difficult to achieve standardized and rapid risk warnings in different engineering scenarios. Summary of the Invention

[0006] In view of this, and in response to the defects and deficiencies of the aforementioned downstream technologies, this invention provides a method for monitoring abnormal rock mass vibration events and assessing rock mass instability risks, in order to solve the problem of early, accurate, and automated early warning of rock mass instability processes, thereby effectively extending the early warning window and improving the objectivity and reliability of risk assessment.

[0007] The method for monitoring abnormal vibration events in rock mass according to the present invention includes:

[0008] A detection model for detecting rock mass vibration signals is constructed. The detection model includes an encoder and a decoder connected to the encoder. The encoder and decoder are each composed of several layers of GRU networks connected in sequence.

[0009] The process of constructing a training set and training the detection model includes:

[0010] 1) Several three-component accelerometers are deployed in the monitored rock mass, and the three-component accelerometers are connected to a ground data processing center. The ground data processing center processes the signals output by the three-component accelerometers, and the processing includes:

[0011] a) Noise reduction processing is performed on the output signal of the three-component accelerometer;

[0012] b) Extract microseismic events from the denoised signal and calculate the arrival times of the first waves of the microseismic events, including:

[0013] The characteristic function CF(t) is defined as follows:

[0014]

[0015] Where u(t) is the denoised signal. Let α represent the Hilbert transform, and α be the weighting coefficient. The denoised signal is sampled through both short and long time windows, and then the ratio of STA to LTA is calculated from the sampled data.

[0016]

[0017] Where: N s N represents the total number of sampling points within the short time window. l This represents the total number of sampling points within the long time window. Let i be the current sampling time point, i be the sequence index of the sampling point within the time window, and Δt be the sampling time interval; when the ratio STA(t) / LTA(t) exceeds the preset threshold T... h If the time is t, then it is determined that a microseismic event has been extracted at the current time t;

[0018] Picking up P-waves: Calculating the Akaike information criterion within a preset time window that includes the time point t of the microseismic event:

[0019]

[0020] In the formula: k is the index of the sampling point within the time window, and N is the total number of sampling points within the time window; This represents the subsequence consisting of the 1st to the kth sampling points within the time window. Let represent the subsequence consisting of the (k+1)th to the Nth sampling points within the time window, and var(⋅) represents calculating the variance; the minimum point k of the AIC function. min That is, the time of arrival of the P wave;

[0021] 2) Constructing physical feature vectors of microseismic events :

[0022]

[0023] in: This represents the absolute timestamp of the i-th microseismic event;

[0024] Let represent the energy of the i-th microseismic event, and its calculation formula is as follows:

[0025]

[0026] In the formula: It is the arrival time of the P-wave of the i-th microseismic event. For the first The duration of the signal for a microseismic event It is the index of the sampling point corresponding to the arrival time of the P wave. It is the duration of the signal. The corresponding total number of sampling points, It is the sampling point loop variable. Let Δt be the amplitude of the nth sampling point, and Δt be the sampling time interval.

[0027] It represents the signal complexity, specifically the ratio of the number of zero-crossings to the peak value of the signal envelope;

[0028] It is the dominant frequency of the signal of the i-th microseismic event, and its calculation formula is as follows:

[0029]

[0030] In the formula: It is the frequency spectrum function of the signal after Fourier transform, where It is a frequency variable;

[0031] It is the centroid of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows:

[0032]

[0033] Let V be the variance of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows:

[0034]

[0035] The wavelet energy entropy of the signal for the i-th microseismic event is calculated using the following formula:

[0036]

[0037] In the formula: E m Let be the energy of the m-th wavelet scale band;

[0038] 3) Reconstruct the phase space of the energy time series of each microseismic event, and calculate the Lyapunov exponent based on the reconstructed phase space;

[0039] 4) Constructing the input features of the detection model :

[0040]

[0041] in, It is the transpose of the physical feature vector of the microseismic event corresponding to time t. The Lyapunov exponent is calculated by reconstructing the phase space based on the energy time series of the microseismic event corresponding to time t.

[0042] Construct time series samples S in ={Z t−L+1 ,…,Z t} where L is the length of the time series samples. A training set is constructed using time series samples obtained during the stable period of the rock mass. The detection model is trained using samples from the training set, and the training objective is to minimize the reconstruction error. :

[0043]

[0044] In the formula, It represents the i-th actual value in the input feature sequence. It is the corresponding reconstructed value output by the detection model;

[0045] Construct a dynamic control chart based on exponentially weighted moving averages:

[0046] a) Use the trained detection model to process all time series samples S from the start of detection to the current time t. in Each time series sample S is obtained. in Reconstruction error ;

[0047] b) Establish a dynamic control chart based on an exponentially weighted moving average, wherein the upper gate control line of the dynamic control chart is... Calculated using the following formula:

[0048]

[0049] In the formula: μ0 is the mean of all reconstruction errors obtained in step a), σ0 is the standard deviation of all reconstruction errors obtained in step a), and k is the control limit width coefficient. It is a smoothing factor;

[0050] Determination of abnormal vibration events based on the gating line in the dynamic control chart:

[0051] Calculate the exponentially weighted moving average value at the current time t. :

[0052]

[0053] exist hour, Initialize to the average reconstruction error on the training set;

[0054] like If the vibration event is detected, it is determined that an abnormal vibration event has occurred; otherwise, it is determined that no abnormal vibration event has occurred.

[0055] Furthermore, step 2) of the noise reduction processing of the monitoring signal output by the three-component accelerometer includes:

[0056] First, wavelet packet decomposition is performed on the raw signal output from the three-component accelerometer to obtain coefficients for different frequency bands; then, the threshold λ is used... j The coefficients for each frequency band are processed; the threshold λ j For the j-th layer, the definition is: ;

[0057] Where: σ j For the noise standard deviation estimate of the wavelet coefficients of the j-th layer, N j denoted as the number of wavelet coefficients in the j-th layer; for the remaining low-frequency trend term, EMD is further applied to decompose the signal into several intrinsic mode functions, and the intrinsic mode function components related to noise are removed.

[0058] Furthermore, the rock mass instability risk assessment method of the rock mass abnormal vibration event monitoring method includes calculating and monitoring the rock mass instability risk index. :

[0059]

[0060] in: Normalize this item;

[0061] It is a normalized outlier score. It is through the reconstruction error at the current moment. The result was obtained after normalization.

[0062] It is the time derivative of the cumulative energy release rate. It is the energy released cumulatively;

[0063] It is the time derivative of the microseismic event occurrence rate. It is the number of microseismic events identified per unit of time;

[0064] Let be the weighting coefficient, satisfying ;

[0065] Will The current instability risk level of the rock mass is obtained by comparing it with the risk level threshold.

[0066] Furthermore, the risk level threshold is determined by statistical learning of historical instability cases; or by constructing a numerical model of the monitored rock mass and then performing rock mass instability inversion simulation.

[0067] The beneficial effects of this invention are:

[0068] 1. This invention captures subtle, early-stage temporal pattern anomalies in microseismic signal sequences through a detection model, enabling early warning signals to be issued before significant changes in macroscopic statistical indicators occur. The early warning is highly advanced and significantly extends the time window from warning to disaster occurrence, thus gaining valuable time for disaster prevention and control.

[0069] 2. This invention abandons the traditional approach of simply viewing microseismic events independently, and innovatively regards microseismic activity as a dynamic temporal process. By using LSTM-Autoencoder to automatically learn its high-dimensional and nonlinear temporal evolution law, it can more profoundly reveal the intrinsic mechanism of rock mass damage accumulation and fracture development, thereby improving the accuracy and reliability of risk assessment.

[0070] 3. This invention automates the entire process from data preprocessing, feature extraction, anomaly detection to risk assessment, greatly reducing reliance on expert experience, lowering the risk of subjective misjudgment, and making risk assessment more standardized and intelligent, making it easier to promote and apply in different engineering sites.

[0071] 4. This invention does not completely replace traditional parameters, but integrates the time-series anomaly patterns discovered by deep learning with traditional macro parameters to construct a comprehensive risk assessment index. This makes the assessment results include both micro-level dynamic precursors and macro-level activity trends, making the assessment system more scientific and comprehensive. Attached Figure Description

[0072] Figure 1 To reduce noise in sensor data.

[0073] Figure 2 This is a time-lapse graph of microseismic events.

[0074] Figure 3 This is a monitoring chart for the rock mass instability risk index.

[0075] Figure 4 This is a graph showing the time of abnormal vibration and the detection of rock mass instability.

[0076] Figure 5 This is a diagram for assessing the risk of rock mass instability. Detailed Implementation

[0077] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0078] The method for monitoring abnormal vibration events in rock mass according to the present invention includes:

[0079] A detection model for detecting rock mass vibration signals is constructed. The detection model includes an encoder and a decoder connected to the encoder. The encoder and decoder are each composed of several layers of GRU (Gate Recurrent Unit) networks connected in sequence.

[0080] The process of constructing a training set and training the detection model includes:

[0081] 1) Several three-component accelerometers are deployed in the monitored rock mass, and the three-component accelerometers are connected to a ground data processing center. The ground data processing center processes the signals output by the three-component accelerometers, and the processing includes:

[0082] a) Noise reduction processing is performed on the monitoring signal output by the three-component accelerometer.

[0083] The raw signal u collected by the sensor raw (t) contains real microseismic events s(t) and complex background noise n(t,θ), where θ represents time- and space-related noise parameters.

[0084] Signal model:

[0085] Traditional fixed filters perform poorly in non-stationary noise environments. This embodiment employs an adaptive noise reduction method based on wavelet thresholding and empirical mode decomposition (EMD).

[0086] First, regarding u raw (t) Perform multi-level wavelet packet decomposition to obtain coefficients W for different frequency bands. j,k , where j is the decomposition level and k is the coefficient index.

[0087] Subsequently, using the threshold λ j The coefficients for each frequency band are processed. The threshold λ... j For the j-th layer, the definition is:

[0088] Where: σ j For the noise standard deviation estimate of the wavelet coefficients of the j-th layer, N j Let be the number of wavelet coefficients in the j-th layer. For the remaining low-frequency trend term, EMD is further applied to decompose the signal into several intrinsic mode functions (IMFs) and remove the IMF components related to noise.

[0089] b) Extract microseismic events from the denoised signal and calculate the arrival times of the first waves of the microseismic events, including:

[0090] The characteristic function CF(t) is defined as follows:

[0091]

[0092] Where u(t) is the denoised signal. Let α represent the Hilbert transform, and α be the weighting coefficient. The denoised signal is sampled through both short and long time windows, and then the ratio of STA to LTA is calculated from the sampled data.

[0093]

[0094] Where: N s N represents the total number of sampling points within the short time window. l This represents the total number of sampling points within the long time window. Let i be the current sampling time point, i be the sequence index of the sampling point within the time window, and Δt be the sampling time interval; when the ratio STA(t) / LTA(t) exceeds the preset threshold T... h If the time is t, then it is determined that a microseismic event has been extracted at the current time t;

[0095] Picking up P-waves: Calculating the Akaike information criterion within a preset time window that includes the time point t of the microseismic event:

[0096]

[0097] In the formula: k is the index of the sampling point within the time window, and N is the total number of sampling points within the time window; This represents the subsequence consisting of the 1st to the kth sampling points within the time window. Let represent the subsequence consisting of the (k+1)th to the Nth sampling points within the time window, and var(⋅) represents calculating the variance; the minimum point k of the AIC function. min This corresponds to the arrival time of the P wave.

[0098] 2) Constructing physical feature vectors of microseismic events :

[0099]

[0100] in: This represents the absolute timestamp of the i-th microseismic event;

[0101] Let represent the energy of the i-th microseismic event, and its calculation formula is as follows:

[0102]

[0103] In the formula: It is the arrival time of the P-wave of the i-th microseismic event. For the first The duration of the signal for a microseismic event It is the index of the sampling point corresponding to the arrival time of the P wave. It is the duration of the signal. The corresponding total number of sampling points, It is the sampling point loop variable. Let Δt be the amplitude of the nth sampling point, and Δt be the sampling time interval.

[0104] It represents the signal complexity, specifically the ratio of the number of zero-crossings to the peak value of the signal envelope;

[0105] It is the dominant frequency of the signal of the i-th microseismic event, and its calculation formula is as follows:

[0106]

[0107] In the formula: It is the frequency spectrum function of the signal after Fourier transform, where It is a frequency variable;

[0108] It is the centroid of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows:

[0109]

[0110] Let V be the variance of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows:

[0111]

[0112] The wavelet energy entropy of the signal for the i-th microseismic event is calculated using the following formula:

[0113]

[0114] In the formula: E m Let be the energy of the m-th wavelet scale band.

[0115] 3) Reconstruct the phase space of the energy time series of each microseismic event, and calculate the Lyapunov exponent based on the reconstructed phase space. The reconstruction process is guided by Takens' Embedding Theorem. The time delay required for phase space reconstruction is determined using the mutual information method, and the embedding dimension required for phase space reconstruction is calculated using the spurious nearest neighbor method.

[0116] 4) Constructing the input features of the detection model :

[0117]

[0118] in, It is the transpose of the physical feature vector of the microseismic event corresponding to time t. The Lyapunov exponent is calculated by reconstructing the phase space based on the energy time series of the microseismic event corresponding to time t.

[0119] Construct time series samples S in ={Z t−L+1 ,…,Z t} where L is the length of the time series sample. A training set is constructed using time series samples obtained during the stable period of the rock mass. The detection model is trained using samples in the training set. The encoder encodes the entire sequence into a fixed-dimensional state vector h. t h t =GRU enc (S in ); the decoder uses h t Starting from the initial state, the input sequence is gradually reconstructed. The training objective is to minimize the reconstruction error. :

[0120]

[0121] In the formula, It represents the i-th actual value in the input feature sequence. It is the corresponding reconstructed value output by the detection model.

[0122] Construct a dynamic control chart based on exponentially weighted moving averages:

[0123] a) Use the trained detection model to process all time series samples S from the start of detection to the current time t. in Each time series sample S is obtained. in Reconstruction error .

[0124] b) Establish a dynamic control chart based on an exponentially weighted moving average, wherein the upper gate control line of the dynamic control chart is... Calculated using the following formula:

[0125]

[0126] In the formula: μ0 is the mean of all reconstruction errors obtained in step a), σ0 is the standard deviation of all reconstruction errors obtained in step a), and k is the control limit width coefficient. This is a smoothing factor.

[0127] Determination of abnormal vibration events based on the gating line in the dynamic control chart:

[0128] Calculate the exponentially weighted moving average value at the current time t. :

[0129]

[0130] exist hour, Initialize to the average reconstruction error on the training set;

[0131] like If the vibration event is detected, it is determined that an abnormal vibration event has occurred; otherwise, it is determined that no abnormal vibration event has occurred.

[0132] The present invention also discloses a method for assessing the risk of rock mass instability according to the rock mass abnormal vibration event monitoring method described in Example 1, which includes calculating and monitoring the rock mass instability risk index. :

[0133]

[0134] in: Normalize this item;

[0135] is the normalized anomaly score, reflecting the degree to which the current microseismic activity pattern deviates from the normal benchmark. is obtained by normalizing the reconstruction error at the current moment.

[0136] is the time derivative of the cumulative energy release rate, reflecting the degree to which the current microseismic activity pattern deviates from the normal benchmark. is the cumulatively released energy.

[0137] is the time derivative of the microseismic event occurrence rate, is the number of microseismic events identified per unit time.

[0138] is the weight coefficient, satisfying ; it can be determined through sensitivity analysis of historical catastrophe cases. For example, in the case of rock burst scenarios, the weight of can be increased to highlight the warning role of abnormal time series patterns.

[0139] Compare with the risk level threshold, and obtain the current instability risk level of the rock mass according to the comparison result.

[0140] The risk level threshold is determined through statistical learning of historical instability cases; or by constructing a numerical model of the monitored rock mass and then performing inverse simulation of rock mass instability.

[0141] In specific implementation, the instability risk level of the rock mass can be set to four levels:

[0142] Level I (low risk, blue): MRI(t) < M1, and the change is gentle. The rock mass is stable.

[0143] Level II (medium risk, yellow): M1 ≤ MRI(t) < M2. The rock mass enters the initial stage of instability and needs to be monitored intensively.

[0144] Level III (high risk, orange): M2 ≤ MRI(t) < M3, or the MRI index rises rapidly in the short term. The risk of rock mass instability is high, and a warning should be issued and some measures should be taken.

[0145] Level IV (extremely high risk, red): MRI(t) ≥ M3. The rock mass is in a critical instability state, and an alarm should be issued immediately and emergency measures should be taken.

[0146] The rock mass abnormal vibration event monitoring and rock mass instability risk assessment method described in the above embodiments was applied to a coal mine. The average burial depth of the coal seam in this mine is 396 meters, and it is determined to have a weak shock tendency. The working face length is 270m, the mining length is 4449.2m, the mining height is 3.7m, the coal release height is 2.6m, and the mining-to-release ratio is 1:0.7. To achieve accurate early warning of rock mass instability, a 16-channel microseismic monitoring system was deployed in the surrounding rock of the mining area, collecting continuous waveform data for 90 days from July 1, 2025 to September 30, 2024. During this period, the system captured 3074 valid microseismic events. The original microseismic signals inevitably contained various mine-specific noises, such as mechanical vibration and electrical interference. The noise reduction steps in the technical solution of this invention were implemented, using an adaptive noise reduction method based on wavelet thresholding and empirical mode decomposition (EMD) to process the signal. Figure 1 As shown, after noise reduction, the signal-to-noise ratio was improved by approximately 10 dB, effectively suppressing background noise while preserving the physical characteristics of the actual microseismic events, laying a solid foundation for subsequent high-precision analysis. Next, high-sensitivity microseismic event detection and accurate first-arrival P-wave acquisition were performed. Using a combination of short-time averaging / long-time averaging (STA / LTA) algorithms and the Akaike Information Criterion (AIC) function, microseismic events were automatically and accurately identified from continuous waveforms, and their first-arrival P-wave times were acquired. Figure 2 The method is visually demonstrated to achieve reliable event detection and accurate P-wave localization even in noisy environments.

[0147] After successfully identifying and extracting each microseismic event, a high-dimensional physical feature vector of the microseismic event was constructed. To comprehensively characterize the physical mechanism of rock mass fracturing, we extracted a multi-dimensional feature vector for each event, containing information in the time domain, frequency domain, and time-frequency domain, specifically including signal energy, complexity, dominant frequency, spectral centroid, spectral variance, and wavelet energy entropy. All events were sorted according to their timestamps to construct a high-dimensional time-series dataset. Figure 3 The evolution of some key features of the working face during the 90-day monitoring period is shown: Subfigures (a) and (b) reveal a clear fluctuating increasing trend in both the daily frequency of microseismic events and energy release; subfigure (c) shows a certain degree of decrease in the dominant event frequency in the later stages, possibly indicating an increase in the rupture scale; subfigure (e) is particularly crucial, showing a significant and sustained increase in the Lyapunov exponent (λ1), which characterizes the degree of system chaos, in the later stages of monitoring (around day 60), a common precursor to dynamic system instability. The evolution of these temporal characteristics indicates that the rock mass micro-fracture activity gradually shifted from an initial random and dispersed state to a later ordered and clustered state, fundamentally changing the dynamic behavior of the microseismic system and providing a rich data foundation for subsequent deep anomaly detection.

[0148] Based on the constructed temporal characteristics, we further performed phase space reconstruction and dynamic invariant extraction of the microseismic system. We selected the energy sequence of microseismic events as the key observations, and using the Takens embedding theorem, we determined the time delay τ using the mutual information method and the embedding dimension m using the spurious nearest neighbor method, successfully reconstructing a phase space that reflects the true dynamics of the system. Figure 4 As shown in (a), the model's reconstruction error of the input sequence increased significantly in the later stages of monitoring, indicating that the current microseismic activity pattern has deviated from the "normal" baseline. We used an active control chart (EWMA) for anomaly detection. Figure 4 (b) The normalized anomaly score clearly shows that the system began to consistently report anomalies from approximately day 65, compared to traditional methods (such as simple energy consumption rate or event frequency, see...). Figure 5 (c)) The response is more advanced and sensitive.

[0149] Finally, we modeled a comprehensive index for rock mass instability risk, classified risk levels, and provided visual early warnings. We constructed a comprehensive rock mass instability risk index (MRI) by weighting and fusing the temporal pattern anomaly score (A'(t)) detected by the GRU autoencoder, the traditional microseismic energy release rate (dE_cum / dt), and the event occurrence rate (dN / dt). Figure 5 (a) The evolution of the MRI index over 90 days was dynamically displayed, and it was divided into four risk levels—low, medium, high, and very high—based on pre-set thresholds (M1=0.3, M2=0.6, M3=0.8). The graph clearly shows that the MRI index began to exceed the medium-risk threshold (yellow) around day 65 and entered the high-risk zone (orange) around day 72, with the system continuously issuing warnings. On day 85, a significant rock mass instability event was indeed recorded on-site (marked by a red inverted triangle in the graph). The warning lead time of this method reached 18 days, far exceeding that of the traditional energy method (approximately 5 days) and the statistical b-value method (approximately 8 days). Figure 4 (d) and Figure 5 As shown in (a). Figure 5 (b) The pie chart further statistically analyzes the proportion of each risk level, showing that the system can effectively distinguish between low-risk conditions most of the time and maintain focus on high-risk conditions. This case demonstrates that the present invention, by deeply integrating the temporal dynamic characteristics of microseismic signals with deep learning intelligent detection, achieves early, accurate, and quantitative assessment of coal mine rock mass instability risks, buying valuable time for on-site pressure relief and mitigation measures, and avoiding potential major safety accidents and economic losses.

[0150] The above description is merely an embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention's specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

[0151] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for monitoring abnormal vibration events in rock mass, characterized in that: include: A detection model for detecting rock mass vibration signals is constructed. The detection model includes an encoder and a decoder connected to the encoder. The encoder and decoder are each composed of several layers of GRU networks connected sequentially. The process of constructing a training set and training the detection model includes: 1) Several three-component accelerometers are deployed in the monitored rock mass, and the three-component accelerometers are connected to a ground data processing center. The ground data processing center processes the signals output by the three-component accelerometers, and the processing includes: a) Noise reduction processing is performed on the output signal of the three-component accelerometer; b) Extract microseismic events from the denoised signal and calculate the arrival times of the first waves of the microseismic events, including: The characteristic function CF(t) is defined as follows: Where u(t) is the denoised signal. Let α represent the Hilbert transform, and α be the weighting coefficient. The denoised signal is sampled through both short and long time windows, and then the ratio of STA to LTA is calculated from the sampled data. Where: N s N represents the total number of sampling points within the short time window. l This represents the total number of sampling points within the long time window. Let i be the current sampling time point, i be the sequence index of the sampling point within the time window, and Δt be the sampling time interval; when the ratio STA(t) / LTA(t) exceeds the preset threshold T... h If the time is t, then it is determined that a microseismic event has been extracted at the current time t; Picking up P-waves: Calculating the Akaike information criterion within a preset time window that includes the time point t of the microseismic event: In the formula: k is the index of the sampling point within the time window, and N is the total number of sampling points within the time window; This represents the subsequence consisting of the 1st to the kth sampling points within the time window. The expression represents the subsequence consisting of the (k+1)th to the Nth sampling points within the time window, and var(⋅) represents calculating the variance; the minimum point k of the AIC function. min That is, the time of arrival of the P wave; 2) Constructing physical feature vectors of microseismic events : in: This represents the absolute timestamp of the i-th microseismic event; Let represent the energy of the i-th microseismic event, and its calculation formula is as follows: In the formula: It is the arrival time of the P-wave of the i-th microseismic event. For the first The duration of the signal for a microseismic event It is the index of the sampling point corresponding to the arrival time of the P wave. It is the duration of the signal. The corresponding total number of sampling points It is the sampling point loop variable. Let Δt be the amplitude of the nth sampling point, and Δt be the sampling time interval. It represents the signal complexity, specifically the ratio of the number of zero-crossings to the peak value of the signal envelope; It is the dominant frequency of the signal of the i-th microseismic event, and its calculation formula is as follows: In the formula: It is the frequency spectrum function of the signal after Fourier transform, where It is a frequency variable; It is the centroid of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows: Let V be the variance of the signal spectrum of the i-th microseismic event, and its calculation formula is as follows: The wavelet energy entropy of the signal for the i-th microseismic event is calculated using the following formula: In the formula: E m The energy of the m-th wavelet scale band; 3) Reconstruct the phase space of the energy time series of each microseismic event, and calculate the Lyapunov exponent based on the reconstructed phase space; 4) Constructing the input features of the detection model : in, It is the transpose of the physical feature vector of the microseismic event corresponding to time t. The Lyapunov exponent is calculated by reconstructing the phase space based on the energy time series of the microseismic event corresponding to time t. Construct time series samples S in ={Z t−L+1 ,…,Z t } where L is the length of the time series samples. A training set is constructed using time series samples obtained during the stable period of the rock mass. The detection model is trained using samples from the training set, and the training objective is to minimize the reconstruction error. : In the formula, It represents the i-th actual value in the input feature sequence. It is the corresponding reconstructed value output by the detection model; Construct a dynamic control chart based on exponentially weighted moving averages: a) Use the trained detection model to process all time series samples S from the start of detection to the current time t. in Each time series sample S is obtained. in Reconstruction error ; b) Establish a dynamic control chart based on an exponentially weighted moving average, wherein the upper gate control line of the dynamic control chart is... Calculated using the following formula: In the formula: μ0 is the mean of all reconstruction errors obtained in step a), σ0 is the standard deviation of all reconstruction errors obtained in step a), and k is the control limit width coefficient. It is a smoothing factor; Determination of abnormal vibration events based on the gating line in the dynamic control chart: Calculate the exponentially weighted moving average value at the current time t. : exist hour, Initialize to the average reconstruction error on the training set; like If the vibration event is detected, it is determined that an abnormal vibration event has occurred; otherwise, it is determined that no abnormal vibration event has occurred.

2. The method for monitoring abnormal vibration events in rock mass according to claim 1, characterized in that: Step 2) The noise reduction processing of the monitoring signal output by the three-component accelerometer includes: First, wavelet packet decomposition is performed on the raw signal output from the three-component accelerometer to obtain coefficients for different frequency bands; then, the threshold λ is used... j The coefficients for each frequency band are processed; the threshold λ j For the j-th layer, the definition is: ; In the formula: σ j For the noise standard deviation estimate of the wavelet coefficients of the j-th layer, N j denoted as the number of wavelet coefficients in the j-th layer; for the remaining low-frequency trend term, EMD is further applied to decompose the signal into several intrinsic mode functions, and the intrinsic mode function components related to noise are removed.

3. The rock mass instability risk assessment method according to claim 1, characterized in that: This includes calculating and monitoring the rock mass instability risk index. : in: Normalize this item; It is a normalized outlier score. It is through the reconstruction error at the current moment. The result was obtained after normalization. It is the time derivative of the cumulative energy release rate. It is the energy released cumulatively; It is the time derivative of the microseismic event occurrence rate. It is the number of microseismic events identified per unit of time; Let be the weighting coefficient, satisfying ; Will The current instability risk level of the rock mass is obtained by comparing it with the risk level threshold.

4. The rock mass instability risk assessment method according to claim 1, characterized in that: The risk level threshold is determined by statistical learning from historical instability cases; or by constructing a numerical model of the monitored rock mass and then performing rock mass instability inversion simulation.

Images