Method and system for identifying micro-seismic signal mode of surrounding rock rock burst risk

Abstract

The invention discloses a method and a system for identifying a micro-seismic signal mode of a surrounding rock rockburst risk, particularly relates to the technical field of deep underground engineering rockburst monitoring, and is used for solving the problem that an existing micro-seismic early warning model lags in identifying a rock mass damage evolution critical state. The method comprises the following steps: collecting a surrounding rock micro-seismic event signal in real time, intercepting an event sequence based on a sliding time window, and calculating a space-time fractal dimension; extracting an event space density field in the window to calculate an energy gradient, and marking a strong interaction cluster of which the energy gradient exceeds a threshold value; storing the fractal dimension change sequence in real time; performing characteristic frequency band filtering on the signals in the window, calculating a ratio of an energy ratio entropy descent rate to a fractal dimension descent rate, and determining a damage acceleration state when the ratio is a supercritical coefficient; generating a rockburst high-risk instruction when the strong interaction cluster and the damage acceleration state coexist and the fractal dimension descent rate is supercritical threshold; the multi-dimensional collaborative identification of the rockburst precursor is realized, and the early warning accuracy and timeliness are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of rockburst monitoring technology in deep underground engineering, and more specifically, to a method and system for identifying microseismic signal patterns of rockburst risk in surrounding rock. Background Technology

[0002] In the field of rockburst risk monitoring in deep underground engineering, microseismic monitoring technology has become a core means of safety status assessment. Existing technologies generally employ early warning models based on energy thresholds and frequency statistics. These models acquire microseismic signals generated by surrounding rock fracturing in real time and trigger alarms based on preset fixed thresholds. This type of method is widely deployed in engineering practice in scenarios such as mines and tunnels. Its technological foundation relies on the assumption of a linear correlation between source parameters (such as energy level and cumulative event frequency) and rock mass instability risk, and a standardized monitoring process has been established.

[0003] Existing early warning models lack the ability to effectively identify abrupt changes in the dynamic state during the rock mass damage evolution process. Since they do not consider the inherent dynamic mode transformation characteristics of the microseismic event sequence, traditional methods cannot capture the essential laws of microseismic activity under critical conditions when the rock mass transitions from steady-state rupture to unstable failure. This results in risk assessment results lagging behind the actual damage state of the rock mass, leading to a significant increase in the false alarm rate under high-risk conditions. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a microseismic signal pattern recognition method and system for the risk of rockburst in surrounding rock to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A method for identifying microseismic signal patterns in surrounding rock burst risk includes:

[0007] S1. Real-time acquisition of microseismic event signals generated by the surrounding rock within the target area;

[0008] S2. Based on a preset sliding time window, extract the microseismic event sequence and calculate the fractal dimension of the spatiotemporal distribution of the microseismic events within the current sliding time window;

[0009] S3. Calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster.

[0010] S4. Update and store the sequence of changes in fractal dimension over time in real time;

[0011] S5. Perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient.

[0012] S6. When a strong interactive cluster exists and the damage acceleration state is confirmed, detect the rate of decrease of the fractal dimension over time. When the rate of decrease exceeds the preset critical rate threshold, generate a high-risk rockburst command.

[0013] Furthermore, real-time acquisition of microseismic event signals generated by the surrounding rock within the target area, including:

[0014] A microseismic sensor array is deployed on the surface of the surrounding rock in the target area, and the original microseismic waveform signal is acquired through the microseismic sensor array.

[0015] The original micro-vibration waveform signal is denoised to eliminate equipment vibration noise and environmental electromagnetic interference noise, resulting in a denoised micro-vibration waveform signal.

[0016] Microseismic events are identified based on the denoised microseismic waveform signals, and the occurrence time and spatial coordinates of each microseismic event are determined.

[0017] Furthermore, based on a preset sliding time window, the sequence of microseismic events is extracted, and the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window is calculated, including:

[0018] Based on the preset sliding time window and the occurrence time of the micro-seismic event, the micro-seismic events within the current sliding time window are extracted to form the micro-seismic event sequence of the current sliding time window;

[0019] Extract the spatial coordinates of each microseismic event in the current sliding time window's microseismic event sequence;

[0020] A spatial point set of microseismic events is constructed based on the extracted spatial location coordinates;

[0021] The fractal dimension of the spatial point set of microseismic events was calculated using the box counting method.

[0022] Furthermore, the fractal dimension of the spatial point set of microseismic events calculated using the box counting method specifically includes: dividing the target area into multiple grid cells, counting the number of grid cells containing at least one microseismic event, and calculating the fractal dimension based on the relationship between the number of grid cells and the size of the grid cells.

[0023] Furthermore, by extracting the spatial density field of microseismic events within the current sliding time window, the energy gradient between adjacent microseismic events is calculated. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster, including:

[0024] Obtain the spatial coordinates and energy values ​​of all microseismic events within the current sliding time window;

[0025] Based on the spatial location coordinates of microseismic events, a spatial density field is generated using a kernel density estimation algorithm.

[0026] For any two microseismic events whose spatial distance is less than the preset neighborhood radius, calculate the energy gradient between the pair of microseismic events;

[0027] When the energy gradient between the pair of microseismic events exceeds the energy gradient threshold, the pair of microseismic events is marked as a strong interaction event.

[0028] A set of microseismic events that are spatially interconnected and contain at least one strong interactive event is labeled as a strong interactive cluster.

[0029] Furthermore, the energy gradient is the absolute value of the energy difference between two microseismic events divided by the spatial distance between the pair of microseismic events.

[0030] Furthermore, the fractal dimension is updated and stored in real time as a sequence of changes over time, including:

[0031] After each sliding time window update, obtain the fractal dimension value calculated for the current sliding time window;

[0032] The data record is formed by associating the end time of the current sliding time window with the corresponding fractal dimension value.

[0033] The data records are appended to the storage structure of the fractal dimension change sequence over time;

[0034] The updated fractal dimension is persistently stored over time using a time series database.

[0035] Furthermore, the microseismic event signal within the sliding time window is bandpass filtered using a preset characteristic frequency band. The energy proportion entropy value of the filtered signal is calculated. When the ratio of the rate of decrease of the energy proportion entropy value to the rate of decrease of the fractal dimension exceeds a preset critical coefficient, the accelerated damage state is confirmed, including:

[0036] The microseismic event signal within the current sliding time window is subjected to bandpass filtering of a preset characteristic frequency band to obtain the filtered microseismic signal;

[0037] The filtered microseismic signal is divided into multiple sub-time periods according to time, and the signal energy value of each sub-time period is calculated.

[0038] Calculate the energy percentage entropy value based on the signal energy values ​​of all sub-time periods;

[0039] Get the rate of decrease of the fractal dimension corresponding to the current sliding time window;

[0040] Calculate the rate of decrease of the entropy value of the energy percentage in the current sliding time window;

[0041] When the ratio of the rate of decrease of the energy percentage entropy to the rate of decrease of the fractal dimension exceeds a preset critical coefficient, the state of accelerated damage is confirmed.

[0042] Furthermore, when a strong interactive cluster exists and the damage acceleration state is confirmed, the rate of decrease of the fractal dimension over time is detected. When the rate of decrease exceeds a preset critical rate threshold, a high-risk rockburst command is generated, including:

[0043] Determine if a strong interaction cluster exists within the current sliding time window;

[0044] Determine whether the accelerated damage state has been confirmed;

[0045] When a strong interactive cluster exists and the damage acceleration state is confirmed, extract the latest consecutive preset number of fractal dimension values ​​from the fractal dimension change sequence over time.

[0046] The rate of decrease of the fractal dimension over time is calculated based on the latest continuously preset number of fractal dimension values;

[0047] When the rate of decrease of the fractal dimension over time exceeds a preset critical rate threshold, a high-risk rockburst warning command is generated and output to the safety monitoring system.

[0048] On the other hand, the present invention provides a microseismic signal pattern recognition system for rockburst risk in surrounding rock, comprising the following modules:

[0049] The microseismic acquisition module is used to acquire microseismic event signals generated by the surrounding rock within the target area in real time.

[0050] The fractal calculation module is used to extract microseismic event sequences based on a preset sliding time window and calculate the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window.

[0051] The cluster marking module is used to calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster.

[0052] The sequence storage module is used to update and store the sequence of changes in fractal dimension over time in real time;

[0053] The damage confirmation module is used to perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient.

[0054] The early warning generation module is used to detect the rate of decrease of the fractal dimension over time when there is a strong interactive cluster and the damage acceleration state is confirmed. When the rate of decrease exceeds the preset critical rate threshold, a high-risk rockburst command is generated.

[0055] Compared with the prior art, the present invention has the following beneficial effects:

[0056] 1. By employing a three-level coupling mechanism of dynamic monitoring of fractal dimension, identification of strong interactive clusters, and collaborative analysis of energy entropy, the limitations of traditional linear early warning models are overcome. Based on a sliding time window, the spatiotemporal fractal dimension of microseismic events is continuously calculated to quantify the degree of rock mass damage accumulation in real time. Through spatial density field energy gradient detection, strong interactive clusters within stress concentration areas are accurately located. By combining the nonlinear correlation determination of characteristic frequency band energy entropy and fractal dimension decrease rate, the dynamic abrupt change of rock mass from steady-state rupture to critical state of unstable failure is identified.

[0057] 2. Through a triple verification mechanism of spatiotemporal dimension coupling (fractal dimension sequence), spatial interaction intensity (strong interaction cluster), and energy release mode (band energy entropy), the timeliness of early warning is significantly improved. Strong interaction cluster markers reveal local stress concentration sources, the decrease in fractal dimension reflects the overall trend of damage expansion, and the sharp drop in energy entropy captures the anomaly of energy release mode. The three work together to form a three-dimensional perception of rockburst precursors, making the early warning window earlier than traditional methods, reducing the missed rate of high-risk conditions, and providing a critical time buffer for support decisions and personnel evacuation. Attached Figure Description

[0058] Figure 1 This is a flowchart of a microseismic signal pattern recognition method for rockburst risk in surrounding rock according to the present invention;

[0059] Figure 2 This is a schematic diagram of the structure of a microseismic signal pattern recognition system for rockburst risk in surrounding rock according to the present invention. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0061] Example 1: Figure 1 This invention presents a method for microseismic signal pattern recognition of rockburst risk in surrounding rock, comprising:

[0062] S1. Real-time acquisition of microseismic event signals generated by the surrounding rock within the target area;

[0063] S2. Based on a preset sliding time window, extract the microseismic event sequence and calculate the fractal dimension of the spatiotemporal distribution of the microseismic events within the current sliding time window;

[0064] S3. Calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster.

[0065] S4. Update and store the sequence of changes in fractal dimension over time in real time;

[0066] S5. Perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient.

[0067] S6. When a strong interactive cluster exists and the damage acceleration state is confirmed, detect the rate of decrease of the fractal dimension over time. When the rate of decrease exceeds the preset critical rate threshold, generate a high-risk rockburst command.

[0068] S1. Real-time acquisition of microseismic event signals generated by the surrounding rock within the target area, specifically implemented as follows:

[0069] During the deployment of the microseismic sensor array on the surrounding rock surface of the target area, the spatial distribution of the sensors is determined based on the three-dimensional geological structural characteristics of the monitoring area, ensuring that the array covers the tunnel roof area, the tunnel sidewalls area, and the area in front of the tunnel excavation. The microseismic sensor array consists of multiple broadband accelerometers, each of which is fixed to the borehole wall by corrosion-resistant metal expansion bolts. The sensor spacing is dynamically adjusted according to the degree of rock joint development: for example, in the fractured rock zone area, the sensor spacing is set to 8 meters; in the intact rock zone area, the sensor spacing is set to 15 meters. When acquiring raw microseismic waveform signals through the microseismic sensor array, the sampling frequency of the signal acquisition system is set to 5000 Hz, with a dynamic range of not less than 120 dB, fully covering the rock fracture characteristic frequency band from 10 Hz to 2000 Hz.

[0070] The noise reduction process for the original microseismic waveform signal involves two core steps: First, a band-stop filter based on the equipment's vibration characteristic frequency is used to eliminate equipment vibration noise. The dominant frequency band of the equipment vibration noise is determined by collecting background signals under no-load operation conditions. For example, for a rock drill, its dominant vibration frequency is determined to be around 85 Hz through spectrum analysis, corresponding to a band-stop filter with a center frequency of 85 Hz and a bandwidth of 20 Hz. Second, adaptive noise cancellation technology is used to suppress environmental electromagnetic interference noise. Specifically, pure interference signals collected by a dedicated reference sensor 2 meters away from the power cable are used as noise samples. The least mean square algorithm is used to iteratively update the coefficients of the 128th-order filter, with an iteration step size set to 0.01. The convergence condition is set to the mean square error change being less than 0.001 for 10 consecutive iterations. The signal-to-noise ratio of the denoised microseismic waveform signal is improved to over 25 dB, forming the denoised microseismic waveform signal. For residual noise, wavelet thresholding is further applied: the db4 wavelet basis function is selected to decompose the signal into 5 levels, and the detail coefficients of each level are processed by a soft thresholding function. The threshold is set to 3 times the estimated standard deviation of the noise in each level.

[0071] When identifying microseismic events based on denoised microseismic waveform signals, a dual-threshold event detection algorithm is employed: a short-time analysis window of 2 milliseconds and a long-time analysis window of 50 milliseconds are set. A valid event trigger point is determined when the ratio of the average signal energy within the short-time window to the average energy within the long-time window exceeds a dynamic threshold of 4.0. Event validity verification requires simultaneous fulfillment of two conditions: the signal energy must be continuously more than three times higher than the background noise energy within 80 milliseconds after the trigger point, and at least three non-collinear sensors in the array must synchronously detect the event waveform. The occurrence time of the microseismic event is determined as the precise timestamp corresponding to the rising edge of the waveform one-quarter before the trigger point. Spatial location coordinates are calculated using the time difference of arrival (TDOA) method: a set of location equations is established using the first arrival time differences of P-waves recorded by at least four sensors in the microseismic sensor array, and the source coordinates are solved by combining this with a pre-calibrated rock mass wave velocity field model. The wave velocity field model was calibrated through active source experiments: Explosives were detonated at a blast point with known three-dimensional coordinates, for example, using 0.5 kg of emulsion explosive in a single blast. The arrival times of the P-waves from each sensor were recorded. A grid search inversion algorithm was used to calculate the optimal wave velocity distribution. The calibrated wave velocity values ​​for sandstone areas ranged from 3200 m / s to 3800 m / s, and for mudstone areas, they ranged from 2800 m / s to 3200 m / s. During the positioning calculation, the time synchronization system used a Global Positioning System (GPS) timing module, achieving a clock synchronization accuracy of 0.1 milliseconds between each acquisition station. When the initial positioning residual exceeded 0.5 milliseconds, an iterative reweighted optimization algorithm was initiated: the theoretical travel time of each sensor was calculated based on the initial solution, and the equations were resolved using the reciprocal of the absolute value of the travel time residual as weighting coefficients. The iteration termination condition was set to a residual change of less than 0.05 milliseconds or the number of iterations reaching 10. Elevation constraints are achieved using sensor pairs deployed on the roof and floor of the tunnel. When the travel time difference between the roof and floor sensors is less than 0.3 milliseconds, the source depth is forcibly constrained within the tunnel profile. Verification of the location results requires that the event-triggered sensor combination include at least one roof sensor and one floor sensor; otherwise, the event record is discarded. The final output is a dataset containing the precise occurrence time and three-dimensional spatial coordinates of the microseismic event.

[0072] S2. Based on a preset sliding time window, extract the microseismic event sequence and calculate the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window. Specifically, this is implemented as follows:

[0073] When extracting microseismic event sequences according to a preset sliding time window, the length of the sliding time window is set based on the characteristic cycle of rock mass stress adjustment. For example, it is set to 8 hours in the tunnel excavation face area and 24 hours in the goaf-affected area. The window sliding step size is set to one-quarter of the window length. When the cumulative number of newly added microseismic events reaches 10% of the total number of events in the current window, a window update operation is triggered. When extracting microseismic events within the current sliding time window, the precise occurrence time data field of the microseismic events output in step S1 is used for filtering: all events whose occurrence time is greater than or equal to the start time of the current window and less than the end time of the current window are included in the sequence, forming the microseismic event sequence of the current sliding time window. The sequence data storage structure includes a unique event number field, a millisecond-level occurrence time field, spatial location coordinate X-axis component values, spatial location coordinate Y-axis component values, spatial location coordinate Z-axis component values, and an energy value field, where the spatial location coordinate values ​​are uniformly in metric units.

[0074] When extracting the spatial coordinates of each microseismic event in the current sliding time window's microseismic event sequence, the X-axis, Y-axis, and Z-axis coordinate values ​​are completely read from the sequence storage structure to form a three-dimensional coordinate triplet dataset. The coordinate values ​​are preserved to two decimal places. The coordinate reference system adopts the tunnel engineering local coordinate system, with its origin defined as the center point of the tunnel entrance. The positive X-axis points to geographic north, the positive Y-axis points to geographic east, and the positive Z-axis points vertically upwards towards the roof. Coordinate transformation is performed when necessary: ​​when the original coordinates are in the geodetic coordinate system, they are converted to the local coordinate system using a pre-calibrated rotation and translation matrix. The matrix parameters are calibrated using measured data from four known control points, and the transformation residual is controlled within 0.1 meters.

[0075] When constructing a spatial point set of microseismic events based on extracted spatial coordinates, each coordinate triplet is converted into a three-dimensional spatial point geometric object, and all point objects constitute an unordered spatial point set. The spatial range of the point set is limited to the minimum bounding cube space of the target area. The boundary of this cube is extended by a safety buffer distance of 2 meters from the original area. For example, when the actual span of the tunnel is 15 meters, the cube size is set to 19 meters multiplied by 19 meters multiplied by the actual height of the tunnel. The point set data structure adopts a quadtree spatial index structure for organization and management, supporting efficient neighborhood query operations. Abnormal data processing is performed during the construction process: when the spatial coordinates of a single point exceed the 20-meter range of the tunnel design boundary, it is marked as a location anomaly event and removed; when the point set size exceeds 1000 points, a uniform thinning algorithm is initiated to retain 500 representative points. The thinning rule is to divide the space into a cubic grid with a side length of 5 meters, and each grid cell randomly retains at most 3 event points.

[0076] When calculating the fractal dimension of the spatial point set of microseismic events using the box counting method, first determine the grid cell size sequence. The size values are generated according to the rule of geometric progression with a common ratio of 0.5 and decreasing in integer powers of 2. The starting size is taken as half of the longest side size of the circumscribed cube of the point set, and the ending size is taken as half of the average spacing value of microseismic events. The intermediate sizes are generated according to a geometric sequence with a common ratio of 0.5. For example, when the size of the roadway area is 50 meters by 30 meters by 8 meters, a six-level sequence of sizes 25 meters, 12.5 meters, 6.25 meters, 3.125 meters, 1.5625 meters, and 0.78125 meters is generated. Divide the target area into cube grid cells of equal size along the X, Y, and Z axes, and count the number of grid cells containing at least one microseismic event. The grid attribution determination rule is: define the grid cell as a semi-open and semi-closed interval spatial range. Xmin represents the minimum boundary coordinate value of the current grid cell in the X-axis direction, and Xmax represents the maximum boundary coordinate value in the X-axis direction; Ymin and Ymax correspond to the Y-axis boundaries; Zmin and Zmax correspond to the Z-axis boundaries. When the microseismic event coordinates (Xe, Ye, Ze) satisfy Xmin ≤ Xe < Xmax and Ymin ≤ Ye < Ymax and Zmin ≤ Ze < Zmax, it is counted into this cell (for example, when the X-axis range of the grid is [10.00 meters, 15.00 meters), the coordinate 12.30 meters is counted in while 15.00 meters is not). Record the number of valid grid cells corresponding to each grid size.

[0077] The process of calculating the fractal dimension based on the relationship between the number of grid cells and the grid cell size is as follows: Take the logarithm to the base 10 of the reciprocal value of the grid cell size as the abscissa value, and take the logarithm of the number of grid cells as the ordinate value to form a data point sequence in a double logarithmic coordinate system. Use the least squares method to fit a straight line: Calculate the mean of the abscissa and the mean of the ordinate. The slope calculation formula is the sum of the products of the deviation of the abscissa from the mean and the deviation of the ordinate from the mean divided by the sum of the squares of the deviation of the abscissa from the mean. Exclude outlier points with too large or too small sizes during fitting: When the grid size is greater than one-fourth of the maximum span of the point set, discard this data point; when the grid size is less than half of the average event spacing, discard this data point. The fractal dimension value is equal to the absolute value of the slope of the fitted straight line, and the calculation result is reserved to three decimal places. The calculation result verification mechanism includes: When the standard deviation of the slope of the fitted straight line exceeds 0.05, readjust the grid size sequence and calculate again; when the difference between the results of two independent calculations exceeds 0.1, return to check the quality of the original positioning data. When the fractal dimension value exceeds 3.0, it is forced to be corrected to 3.0, and when it is less than 1.0, it is corrected to 1.0 to ensure compliance with the theoretical limit range of the spatial distribution of rock mass fractures. The calculation result is real-time associated with the end time of the time window and deposited into the time series database.

[0078] The boundary handling rules in mesh generation are as follows: when the target region size is not divisible by the mesh size, the region boundary is extended to the nearest integer multiple of the mesh size. The extended region does not contain actual event points and does not affect the mesh count statistics. The quality control index for fractal dimension calculation requires a fitted correlation coefficient squared greater than 0.98; otherwise, the mesh size sequence is regenerated or anomalies in the spatial distribution of the point set are checked. Memory usage is monitored in real time during the calculation. When the point set size exceeds 5000 points, it automatically switches to a stratified sampling algorithm, maintaining the original spatial distribution characteristics of the sampling rate. After each window slide, the system automatically generates a graph showing the relationship between the mesh size sequence and the number of meshes, allowing engineers to verify the validity of the calculation process.

[0079] S3. Calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster. Specifically, the implementation is as follows:

[0080] When acquiring the spatial coordinates and energy values ​​of all microseismic events within the current sliding time window, complete data fields are read from the microseismic event sequence storage structure generated in step S2. Spatial coordinates include three data items: X-axis component values, Y-axis component values, and Z-axis component values, all in metric units with numerical precision retained to two decimal places. Energy values ​​are measured in joules. The amplitude of the original sensor output signal is converted to energy value using a calibration formula. The calibration coefficient is determined through laboratory blasting tests: for example, detonating 0.2 kg of emulsion explosive at a known location, recording the peak value of the sensor output voltage, and establishing the energy conversion relationship. Validity checks are performed during data acquisition: energy values ​​less than 10 joules are considered background noise events and discarded; spatial coordinates exceeding 20 meters of the tunnel design boundary are marked as abnormal location events and excluded.

[0081] When generating the spatial density field based on the spatial location coordinates of microseismic events, a kernel density estimation algorithm is used. The specific operation process is as follows: First, the kernel function type is determined to be a three-dimensional Gaussian kernel function, and the kernel function expression is an exponential decay function of spatial distance. The kernel function bandwidth parameter is adaptively calculated according to the Silverman criterion: the standard deviation of the position coordinates of all microseismic events in the XYZ directions is calculated, and the median of the standard deviation is multiplied by a coefficient of 1.06, and then multiplied by the negative 1 / 5 power of the total number of events as the baseline bandwidth. For example, when there are 200 events and the median standard deviation is 5.2 meters, the baseline bandwidth is calculated as 5.2 multiplied by 1.06 multiplied by 200 to the power of -0.2, approximately equal to 3.8 meters. The target area is divided into cubic grid nodes with a side length of 0.5 meters. The density value of each grid node is equal to the weighted sum of the contributions of all microseismic events to that node. The contribution value of a single event is calculated as the value of the kernel function at the distance from the event location to the grid node. The density field output is a three-dimensional matrix data format, with the matrix index corresponding to the spatial grid coordinates.

[0082] The operation of calculating the energy gradient between adjacent microseismic events includes the following steps: Define a preset neighborhood radius of 5 meters, which is set based on twice the average spacing between microseismic events. Iterate through all combinations of microseismic events, and perform the calculation when the spatial distance between two events is less than the preset neighborhood radius. The spatial distance is calculated using the Euclidean distance formula: the distance equals the square root of the sum of the squares of the differences in the X-axis, Y-axis, and Z-axis coordinates. The energy gradient is calculated as the absolute value of the energy difference between the events divided by the spatial distance. For example, if event A has an energy value of 150 joules, event B has an energy value of 80 joules, and the spatial distance is 3.2 meters, the energy gradient is the absolute value 150 - 80 divided by 3.2, which equals 21.875 joules per meter. The calculation result is rounded to two decimal places.

[0083] The rule for marking strongly interactive events is as follows: An energy gradient threshold of 15 joules per meter is set, determined based on historical data statistics—100 sets of rock mass fracture case data are selected, and the 90th percentile value of the energy gradient distribution during the stable fracture stage is calculated. When the calculated energy gradient of an event pair is greater than or equal to the threshold, both events are marked as strongly interactive events, and a bidirectional association is established for the event pair. The marked data storage structure includes an event number, associated event number, and energy gradient value fields.

[0084] The process of constructing a strong interaction cluster is as follows: First, an empty cluster set is created. Then, all marked strong interaction events are traversed, and an event association network is established using a depth-first search algorithm. Specifically, starting with any unvisited strong interaction event, all unvisited events in its associated event list are searched, recursively expanding until no new events are added. A single cluster is defined as a set of spatially connected events. The connectivity criterion is that there exists a path between any two events in the set consisting of a pair of strong interaction events, and the distance between adjacent events in the path does not exceed twice the preset neighborhood radius, i.e., 10 meters. A cluster is considered a valid strong interaction cluster when it contains at least three microseismic events. The cluster's spatial extent is calculated as the smallest bounding cube of the coordinates of all events within the cluster. The output data includes the cluster number, the list of event members, and the spatial extent coordinates.

[0085] The cluster validity verification mechanism includes: checking for anomalies in energy value calibration data when the average energy gradient of events within the cluster is below a threshold; and checking the reasonableness of the neighborhood radius setting when the cluster spatial volume exceeds 1000 cubic meters. Before outputting results, a spatial aggregation degree verification is performed: the ratio of the average distance between events within the cluster to the expected value of the random distribution is calculated; a ratio less than 0.6 confirms a valid spatial aggregation pattern. The final output data is associated and stored in the rock mass stability analysis database using time window identifiers.

[0086] The bandwidth optimization process of the kernel density estimation algorithm is as follows: After initial bandwidth calculation, cross-validation adjustment is performed. The target region is divided into training and validation sets, and the optimization objective is to maximize the log-likelihood value of the validation set, searching for the optimal value within a 20% fluctuation range above and below the baseline bandwidth. The grid size is set according to the bandwidth value: when the bandwidth is greater than 3 meters, the grid size is one-fifth of the bandwidth; when the bandwidth is less than 3 meters, the grid size is a fixed value of 0.5 meters. For density field visualization, a linear interpolation algorithm is used to generate isosurfaces, and the spacing between isosurfaces is dynamically determined based on the density value distribution range.

[0087] Anomaly handling in energy gradient calculation includes: when the spatial distance between two events is less than 0.3 meters in positioning accuracy, they are considered duplicate positioning events and merged; when the energy difference exceeds twice the maximum recorded value, a sensor calibration status check is triggered. The real-time update mechanism for strong interaction event markers is: when a new event is added to the current time window, only its correlation with existing strong interaction events is calculated to avoid full duplicate calculations.

[0088] The cluster boundary is determined as follows: the standard deviation of the coordinates of all events within the cluster is calculated, and an ellipsoidal spatial range is established with the mean location as the center and three times the standard deviation as the boundary. Cluster spatial connectivity is verified using the Delaunay triangulation algorithm: a connected structure is confirmed when there is a triangular network connecting any two points within the cluster and all side lengths are less than twice the radius of their neighboring regions. The cluster merging rule is as follows: when the overlap of the outer cubes of two clusters exceeds 30%, the correlation of events within the overlapping area is checked; if strong interactive event pairs exist, the clusters are merged.

[0089] S4. Real-time update and storage of the fractal dimension change sequence over time, specifically implemented as follows:

[0090] Upon completion of each sliding time window update operation, the fractal dimension value acquisition process is automatically triggered, retrieving the fractal dimension value corresponding to the current sliding time window from the fractal dimension calculation result storage area of ​​step S2. This value is a double-precision floating-point data type, retaining three decimal places of precision, with a value range constrained between 1.000 and 3.000. The acquisition process includes data validity verification: if the fractal dimension value exceeds the theoretical range, it is re-verified; if three consecutive acquisition failures occur, a system abnormal alarm is triggered and an error log is recorded. The read operation is implemented using memory-mapped file technology, with access latency controlled within 10 milliseconds.

[0091] When associating the end time of the current sliding time window with the corresponding fractal dimension value to form a data record, the timestamp uses an international standard time format accurate to the millisecond level. The calculation rule for the end time of the time window is as follows: the initial window end time is set to the system startup time plus the window length value, and the end time of subsequent windows is accumulated according to the window sliding step size. The data record structure is defined as a binary tuple, containing a timestamp field and a fractal dimension value field. For example, a typical record is represented as (2023-08-15T14:30:45.250Z, 2.375). The record generation process performs time consistency verification: when the deviation between the current time and the window end time exceeds 10% of the sliding step size, the time base is automatically corrected and a data flag is marked.

[0092] When appending data records to a fractal dimension-dependent time-varying sequence storage structure, a circular buffer data management mechanism is employed. The buffer capacity is set to store 90 days of continuous data. When the data volume reaches 95% of the buffer capacity, data archiving and transfer are automatically initiated. The append operation includes a sequence continuity check: verifying that the timestamp of the new record is strictly greater than the timestamp of the most recent record; when timestamp rollback is abnormal, a data conflict resolution algorithm is activated to retain the valid record with the largest timestamp. The storage structure index uses a B+ tree organization, supporting fast retrieval by time range. The index key is the timestamp field, and the leaf nodes store the physical address of the data record.

[0093] When persistently storing the updated fractal dimension over time using a time-series database, a dual-write mechanism ensures data reliability. The primary storage uses a columnar time-series database, with table structures including a timestamp primary key column, a fractal dimension value column, and a data quality flag column. The secondary storage uses a transaction log file to record all data change operations. Persistence operations are performed asynchronously after data appending, with a write batch size of 100 records and a write timeout threshold of 5 seconds. Database connection parameters are managed through an encrypted configuration file, and the connection pool size is dynamically adjusted based on the data update frequency; for example, when the window sliding step is 2 hours, the connection pool maintains 5 active connections.

[0094] The data persistence process includes integrity checks: after each write, a record count verification is performed to ensure that the number of records in the memory buffer matches the number of records stored in the database; a data checksum calculation is performed weekly to compare the data consistency between memory and storage media. The storage system implements a regular backup strategy: incremental backups are performed daily at midnight, and a full backup is performed every Sunday, with a backup data retention period of three years. Data compression uses a time-series dedicated compression algorithm, applying delta-delta encoding combined with simple 8-bit compression to fractal dimension numerical sequences, achieving a compression rate of up to 20% of the original data size.

[0095] The historical data query interface is implemented as follows: Define time range query parameters; during retrieval, first search for recent data in the memory buffer; if no match is found, query the database. Query response time requirements are: less than 500 milliseconds for data within three months, and less than 3 seconds for data throughout the year. The data visualization service pushes sequence updates in real time via the WebSocket protocol, with the push frequency synchronized with the window sliding step size; for example, with a step size of 2 hours, the latest data point is pushed every two hours. Data access permissions are subject to hierarchical control: raw data is only accessible to authorized engineers, while aggregated analysis results are available to monitoring personnel.

[0096] Storage system monitoring metrics include: data ingestion latency, storage space utilization, and query response time. When the ingestion latency exceeds 50% of the sliding step size, the database compute nodes are automatically expanded; when storage space utilization exceeds 85%, the old data archiving and cleanup process is triggered. The system generates a daily data health report, including parameters such as record integrity statistics, outlier ratios, and storage efficiency metrics. All configuration parameters are centrally maintained through the management console, and parameter modifications are logged in an audit log with version rollback support.

[0097] S5. Perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient. Specifically, the implementation is as follows:

[0098] When performing bandpass filtering on microseismic event signals within the current sliding time window using a preset characteristic frequency band, the characteristic frequency band range is first determined to be 100 Hz to 300 Hz. This range is set based on the spectral characteristics of rock mass fracture: typical fracture signals are collected through laboratory rock sample failure tests, and the mean of the energy concentration frequency band plus or minus two standard deviations is extracted as the boundary value. The filtering process uses a zero-phase finite-length unit impulse response filter. The filter order is adaptively calculated based on the sampling frequency: the sampling frequency is divided by the center frequency of the characteristic frequency band and then multiplied by a coefficient of 20. For example, when the sampling frequency is 5000 Hz, the calculated order is 5000 divided by 200 multiplied by 20, which equals 500. Signal preprocessing is performed before filtering: DC component removal, sensor frequency response deviation compensation, and elimination of 50 Hz power frequency interference. Each microseismic event signal is processed independently, and the output filtered signal retains the original sampling rate of 5000 Hz, with the data length consistent with the original event signal.

[0099] When the filtered microseismic signal is divided into multiple sub-time periods, the number of sub-time periods is dynamically determined based on the signal duration: periods shorter than 2 seconds are divided into 4 equal-length segments; periods between 2 and 10 seconds are divided into 8 segments; and periods longer than 10 seconds are divided into 2 additional segments for every 5 seconds of increase. Time division is accurate to the sampling point level. When the signal length is not divisible by an integer, the last segment is allowed to have a maximum reduction of 10% in the number of sampling points. The start time of each sub-time period is stored in association with the absolute occurrence time of the microseismic event, with time synchronization errors controlled within 0.1 milliseconds.

[0100] The signal energy value for each sub-time period is calculated using the square integral method: first, the signal amplitude of all sampling points within the sub-time period is squared; then, all squared values ​​are summed to obtain the original energy value; finally, this is multiplied by the sampling time interval of 0.0002 seconds to convert it into a physical energy value, with the unit uniformly set to microjoules. Signal quality checks are performed during the calculation process: when a sub-time period contains saturated sampling points, it is marked as a signal distortion segment and replaced with the average of adjacent time periods; when the signal-to-noise ratio is below 10 dB, wavelet denoising is initiated before further calculation. The energy values ​​of all sub-time periods are arranged chronologically to form an energy sequence, with the sequence index corresponding to the chronological order.

[0101] When calculating the energy percentage entropy value based on the signal energy values ​​of all sub-time periods, the proportion of energy in each sub-time period to the total energy is first calculated: the energy value of a single sub-time period is divided by the sum of the energy values ​​of all sub-time periods. The formula for calculating the energy percentage entropy value is the sum of the negative values ​​of the energy proportions of all sub-time periods multiplied by the base-2 logarithmic values ​​of those proportions, with the result rounded to four decimal places. Entropy calculation includes boundary handling: when the energy proportion of a certain sub-time period is less than 0.0001, it is calculated as 0.0001 to avoid logarithmic operation errors; when the total energy value is less than 1 microjoule, it is considered an invalid event and the entropy calculation is skipped. The physical meaning of entropy value characterizes the uniformity of energy distribution, with a numerical range between 0 and 10.

[0102] When obtaining the rate of decrease of the fractal dimension corresponding to the current sliding time window, data is extracted from the fractal dimension time series storage structure of step S4. The extraction rule is: take all fractal dimension values ​​within a time period three times the window length preceding the end time of the current time window to form a time series dataset. The rate of decrease is calculated using a linear regression method: with the time point as the independent variable and the fractal dimension value as the dependent variable, the slope of the fitted line is the rate of decrease, and the unit is defined as the change per day. Before calculation, data validity screening is performed: invalid fractal dimension points caused by abnormal positioning are removed; when there are fewer than 5 data points, the time range is expanded to five times the window length. The rate value is retained to three decimal places, with positive values ​​indicating an increase and negative values ​​indicating a decrease.

[0103] When calculating the rate of decrease of the energy percentage entropy value within the current sliding time window, an entropy dataset with the same time span is used. After calculating the energy percentage entropy value for all valid microseismic events within the current window, the events are sorted by their occurrence time to form a time series. The calculation rule for the rate of decrease is consistent with that for the fractal dimension rate of decrease: a linear regression is performed with the event occurrence time as the independent variable and the entropy value as the dependent variable, and the slope is measured in hourly changes. Special handling includes: when the number of events within the window is less than 10, a null value is returned and the event is not included in subsequent judgments; when the standard deviation of the entropy value sequence is less than 0.01, the event is considered to be in a constant state and the rate is not calculated.

[0104] When the ratio of the rate of decrease in energy percentage entropy to the rate of decrease in fractal dimension exceeds a preset critical coefficient, an accelerated damage state is confirmed. The preset critical coefficient is set at 2.5, determined based on historical rockburst case statistics: data from 30 rockbursts occurring 72 hours prior to their occurrence is analyzed, and the 95th percentile of the ratio of the rate of decrease in energy percentage entropy to the rate of decrease in fractal dimension is used. The ratio is calculated by dividing the absolute value of the rate of decrease in energy percentage entropy by the absolute value of the rate of decrease in fractal dimension, with the result rounded to two decimal places. The state confirmation includes additional conditions: the rate of decrease in fractal dimension must be less than -0.05 per day, and the rate of decrease in energy percentage entropy must be less than -0.2 per hour. The state labeling result, associated with the end of the time window, is stored in the disaster early warning database and pushed to the monitoring terminal in real time via a message queue.

[0105] The verification mechanism includes: when three consecutive windows trigger the damage acceleration marker, the warning level is raised to orange; when the ratio exceeds the critical coefficient but the fractal dimension does not continue to decrease, the sensor network calibration procedure is initiated. The system automatically generates a dual-rate trend overlay graph, with the horizontal axis representing time, the left vertical axis representing the fractal dimension decrease rate, and the right vertical axis representing the energy entropy decrease rate. The degree of separation between the two curve trends intuitively reflects the damage development status. The historical data backtracking function supports retrieving the dual-rate ratio curve for any time period, assisting engineers in verifying the accuracy of status judgments.

[0106] S6. When a strong interactive cluster exists and the damage acceleration state is confirmed, the rate of decrease of the fractal dimension over time is detected. When the rate of decrease exceeds a preset critical rate threshold, a high-risk rockburst command is generated. The specific implementation is as follows:

[0107] When determining whether a strongly interactive cluster exists within the current sliding time window, the record corresponding to the current time window identifier is retrieved from the strongly interactive cluster storage structure in step S3. The retrieval criteria are that the cluster validity flag is true and the deviation between the cluster generation timestamp and the current window end time is less than half of the window sliding step. The existence determination of a strongly interactive cluster is based on the following: when at least one valid cluster record is retrieved and the cluster contains three or more events, the existence status is returned; when the cluster record is empty or the total number of cluster events is less than three, the non-existence status is returned. The retrieval process performs data timeliness verification: cluster record timestamps exceeding twice the current window end time are automatically marked as expired data and excluded.

[0108] When determining whether an accelerated damage state has been confirmed, the latest record in the damage state database from step S5 is queried. The query rule is: retrieve the latest state marker record within a sliding step time period one window before the end of the current time window. Conditions for a confirmed state include: the accelerated damage state flag is true, the state confidence score is greater than 0.85, and the state duration exceeds six hours. The state database uses a time-partitioned table for storage, and the query response time is controlled within 100 milliseconds. If a record is missing due to data transmission delay, the state recalculation service is initiated and waits for a maximum of two sliding steps.

[0109] When a strong interactive cluster exists simultaneously and the damage acceleration state is confirmed, the latest consecutive preset number of fractal dimension values ​​are extracted from the fractal dimension time series storage structure of step S4. The preset number is fixed at 10 data points, and this value is set based on the duration analysis of rockburst precursor characteristics: statistical historical rockburst cases show that the continuous decrease in fractal dimension covers an average of eight time windows, plus two window margins to form the final set value. The extraction rule is: take the latest 10 valid records in reverse chronological order, and the time interval between records is strictly equal to the window sliding step size. Data validity filtering criteria include: excluding records with abnormal quality flags; when the fractal dimension value exceeds the range of 1.000 to 3.000, it is replaced with the average of the preceding and following records. The extraction operation is performed in batches using a database cursor, and the extraction time for a single extraction does not exceed 50 milliseconds.

[0110] When calculating the rate of decrease of the fractal dimension over time based on the latest continuously preset number of fractal dimension values, a weighted linear regression algorithm is used. The specific calculation process is as follows: the end time of the time window is used as the independent variable, and the time interval is in hours; the fractal dimension value is the dependent variable. A time weighting coefficient is set, with the latest data point having a weight of 1.0, and the weight decreasing by a factor of 0.1 for each previous point in reverse chronological order. The formula for calculating the rate of decrease is: the numerator is the weighted sum of the deviation of time from the mean multiplied by the deviation of the fractal dimension from the mean; the denominator is the weighted sum of the squares of the deviations of time from the mean; and the result is defined as the daily change in fractal dimension. Calculation example: when the latest 10 fractal dimensions continuously decrease from 2.85 to 2.15 over a time span of 48 hours, the calculated rate of decrease is -0.0146 per day. The rate value is rounded to four decimal places, and a negative value indicates a downward trend.

[0111] When the rate of decrease of the fractal dimension over time exceeds a preset critical rate threshold, a high-risk rockburst warning command is generated and output to the safety monitoring system. The preset critical rate threshold is set to -0.012 per day, determined through retrospective analysis of 35 rockburst accidents: the 99th percentile of the fractal dimension decrease rate within the 24 hours prior to the rockburst. The threshold comparison rule is: when the calculated decrease rate is less than or equal to the preset critical rate threshold (i.e., the absolute value of the negative value is larger), the warning command is triggered. Command generation includes multiple verifications: requiring that the rates at three consecutive calculation points are all below the threshold; and that the maximum size of the current window's highly interactive cluster exceeds one-third of the roadway cross-sectional size. The warning command data structure includes four parts: command code, risk level, time window identifier, and key parameter values. For example, the high-risk command code is "RB_ALERT_LEVEL3".

[0112] The command output to the safety monitoring system uses a standard application programming interface protocol, and the transmission process includes a three-way handshake confirmation mechanism. The output content is a structured data message containing the following fields: the command type field is fixed as "rockburst red warning"; the risk location field is taken from the average spatial coordinates of the strong interactive cluster; the risk time window field is the current window's end time plus twice the sliding step size; and the confidence level field is calculated based on the abnormal deviation of the parameters. Upon receiving the command, the safety system activates the emergency protocol: automatically triggering the evacuation command for personnel on the work face, starting the support system pressurization program, recording the decision log, and uploading it to the cloud audit platform.

[0113] The false alarm prevention mechanism includes: automatically downgrading the warning level when the fractal dimension rises in the subsequent two windows after a warning is triggered; and initiating a system-wide diagnostic calibration when more than three warnings are triggered in a single day. Historical warning records are stored in conjunction with geological structure maps, supporting spatial location clustering analysis. The system generates a monthly warning effectiveness report, statistically analyzing false alarm and missed alarm rates. When the false alarm rate exceeds 5%, the critical rate threshold is automatically adjusted to 3%.

[0114] Example 2: Figure 2 A schematic diagram of a microseismic signal pattern recognition system for rockburst risk in surrounding rock is provided according to the present invention. The microseismic signal pattern recognition system for rockburst risk in surrounding rock includes the following modules:

[0115] The microseismic acquisition module is used to acquire microseismic event signals generated by the surrounding rock within the target area in real time.

[0116] The fractal calculation module is used to extract microseismic event sequences based on a preset sliding time window and calculate the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window.

[0117] The cluster marking module is used to calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster.

[0118] The sequence storage module is used to update and store the sequence of changes in fractal dimension over time in real time;

[0119] The damage confirmation module is used to perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient.

[0120] The early warning generation module is used to detect the rate of decrease of the fractal dimension over time when there is a strong interactive cluster and the damage acceleration state is confirmed. When the rate of decrease exceeds the preset critical rate threshold, a high-risk rockburst command is generated.

[0121] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0122] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0123] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and inventive constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0124] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0125] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0126] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0127] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for identifying microseismic signal patterns in response to rockburst risk in surrounding rock, characterized in that, include: S1. Real-time acquisition of microseismic event signals generated by the surrounding rock within the target area; S2. Based on a preset sliding time window, extract the microseismic event sequence and calculate the fractal dimension of the spatiotemporal distribution of the microseismic events within the current sliding time window; S3. Calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster. S4. Update and store the sequence of changes in fractal dimension over time in real time; S5. Perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient. S6. When a strong interactive cluster exists and the damage acceleration state is confirmed, detect the rate of decrease of the fractal dimension over time. When the rate of decrease exceeds the preset critical rate threshold, generate a high-risk rockburst command.

2. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 1, characterized in that, Real-time acquisition of microseismic event signals generated by the surrounding rock within the target area, including: A microseismic sensor array is deployed on the surface of the surrounding rock in the target area, and the original microseismic waveform signal is acquired through the microseismic sensor array. The original micro-vibration waveform signal is denoised to eliminate equipment vibration noise and environmental electromagnetic interference noise, resulting in a denoised micro-vibration waveform signal. Microseismic events are identified based on the denoised microseismic waveform signals, and the occurrence time and spatial coordinates of each microseismic event are determined.

3. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 2, characterized in that, Based on a preset sliding time window, the sequence of microseismic events is extracted, and the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window is calculated, including: Based on the preset sliding time window and the occurrence time of the micro-seismic event, the micro-seismic events within the current sliding time window are extracted to form the micro-seismic event sequence of the current sliding time window; Extract the spatial coordinates of each microseismic event in the current sliding time window's microseismic event sequence; A spatial point set of microseismic events is constructed based on the extracted spatial location coordinates; The fractal dimension of the spatial point set of microseismic events was calculated using the box counting method.

4. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 3, characterized in that, The fractal dimension of a spatial point set of microseismic events calculated using the box counting method specifically includes: dividing the target area into multiple grid cells, counting the number of grid cells containing at least one microseismic event, and calculating the fractal dimension based on the relationship between the number of grid cells and the size of the grid cells.

5. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 3, characterized in that, The energy gradient between adjacent microseismic events is calculated by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds an energy gradient threshold, it is marked as a strong interaction cluster, including: Obtain the spatial coordinates and energy values ​​of all microseismic events within the current sliding time window; Based on the spatial location coordinates of microseismic events, a spatial density field is generated using a kernel density estimation algorithm. For any two microseismic events whose spatial distance is less than the preset neighborhood radius, calculate the energy gradient between the pair of microseismic events; When the energy gradient between the pair of microseismic events exceeds the energy gradient threshold, the pair of microseismic events is marked as a strong interaction event. A set of microseismic events that are spatially interconnected and contain at least one strong interactive event is labeled as a strong interactive cluster.

6. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 5, characterized in that, The energy gradient is the absolute value of the energy difference between two microseismic events divided by the spatial distance between the two microseismic events.

7. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 5, characterized in that, Real-time updates and storage of the fractal dimension change sequence over time, including: After each sliding time window update, obtain the fractal dimension value calculated for the current sliding time window; The data record is formed by associating the end time of the current sliding time window with the corresponding fractal dimension value. The data records are appended to the storage structure of the fractal dimension change sequence over time; The updated fractal dimension is persistently stored over time using a time series database.

8. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 7, characterized in that, Bandpass filtering of microseismic event signals within a sliding time window is performed on a preset characteristic frequency band. The energy proportion entropy value of the filtered signal is calculated. When the ratio of the rate of decrease of the energy proportion entropy value to the rate of decrease of the fractal dimension exceeds a preset critical coefficient, the accelerated damage state is confirmed, including: The microseismic event signal within the current sliding time window is subjected to bandpass filtering of a preset characteristic frequency band to obtain the filtered microseismic signal; The filtered microseismic signal is divided into multiple sub-time periods according to time, and the signal energy value of each sub-time period is calculated. Calculate the energy percentage entropy value based on the signal energy values ​​of all sub-time periods; Get the rate of decrease of the fractal dimension corresponding to the current sliding time window; Calculate the rate of decrease of the entropy value of the energy percentage in the current sliding time window; When the ratio of the rate of decrease of the energy percentage entropy to the rate of decrease of the fractal dimension exceeds a preset critical coefficient, the state of accelerated damage is confirmed.

9. The microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in claim 8, characterized in that, When a strong interactive cluster exists and the damage acceleration state is confirmed, the rate of decrease of the fractal dimension over time is detected. When the rate of decrease exceeds a preset critical rate threshold, a high-risk rockburst command is generated, including: Determine if a strong interaction cluster exists within the current sliding time window; Determine whether the accelerated damage state has been confirmed; When a strong interactive cluster exists and the damage acceleration state is confirmed, extract the latest consecutive preset number of fractal dimension values ​​from the fractal dimension change sequence over time. The rate of decrease of the fractal dimension over time is calculated based on the latest continuously preset number of fractal dimension values; When the rate of decrease of the fractal dimension over time exceeds a preset critical rate threshold, a high-risk rockburst warning command is generated and output to the safety monitoring system.

10. A microseismic signal pattern recognition system for rockburst risk in surrounding rock, used to implement the microseismic signal pattern recognition method for rockburst risk in surrounding rock as described in any one of claims 1-9, characterized in that, Includes the following modules: The microseismic acquisition module is used to acquire microseismic event signals generated by the surrounding rock within the target area in real time. The fractal calculation module is used to extract microseismic event sequences based on a preset sliding time window and calculate the fractal dimension of the spatiotemporal distribution of microseismic events within the current sliding time window. The cluster marking module is used to calculate the energy gradient between adjacent microseismic events by extracting the spatial density field of microseismic events within the current sliding time window. When the energy gradient exceeds the energy gradient threshold, it is marked as a strong interaction cluster. The sequence storage module is used to update and store the sequence of changes in fractal dimension over time in real time; The damage confirmation module is used to perform bandpass filtering on the microseismic event signal within the sliding time window using a preset characteristic frequency band, calculate the energy proportion entropy value of the filtered signal, and confirm the accelerated damage state when the ratio of the energy proportion entropy value decrease rate to the fractal dimension decrease rate exceeds a preset critical coefficient. The early warning generation module is used to detect the rate of decrease of the fractal dimension over time when there is a strong interactive cluster and the damage acceleration state is confirmed. When the rate of decrease exceeds the preset critical rate threshold, a high-risk rockburst command is generated.

Images