Deep rock mass breaking dynamic prediction method based on microseismic monitoring

Abstract

The present application mainly relates to the technical field of geological monitoring, in order to predict the deep rock mass crushing index in real time and accurately, the present application provides a kind of deep rock mass crushing dynamic prediction method based on microseismic monitoring, its core is: collecting microseismic monitoring data, and integrating according to time sequence, constructs multidimensional feature input sequence;Rock mass crushing key time delay feature matrix is extracted from the multidimensional feature input sequence;The key time delay feature matrix is input into the prediction model containing the Transformer encoder to extract deep time sequence features, and the output high-dimensional features are mapped through the output layer full connection network, to obtain the future time rock mass crushing index prediction value, so that the prediction model can accurately capture the nonlinear time delay correlation between microseismic monitoring parameters and rock mass crushing index under complex geological conditions, and finally realize the significant improvement of rock mass crushing index prediction accuracy.

Description

Technical Field

[0001] This invention mainly relates to the field of geological monitoring technology, and in particular to a method for dynamic prediction of deep rock mass fracturing based on microseismic monitoring. Background Technology

[0002] Predicting deep rock mass fracturing is a core technical support for ensuring safe construction, optimizing support schemes, and improving construction efficiency in deep mining, tunnel excavation, and other engineering projects. Deep rock masses are affected by complex geological conditions (such as differences in ground stress and lithological variations), and there is a significant time delay effect between their microseismic activity (frequency, energy, focal depth, etc.) and the degree of fracturing. Accurately capturing this effect is key to improving prediction accuracy. However, the prediction accuracy and precision of rock mass fracturing indicators in existing technologies are limited, mainly because:

[0003] (1) The complexity and time lag of physical correlation. The process from microcrack initiation and convergence to macroscopic fracturing and instability of deep rock masses is a nonlinear energy accumulation and release process. There is a significant time lag effect, which is controlled by geological conditions, between microseismic monitoring parameters such as energy release rate and focal depth distribution and the final degree of fracturing. Under different lithologies and geostress states, the same microseismic activity pattern may correspond to completely different fracturing risks. Traditional statistical methods based on fixed time windows (such as moving average and exponential smoothing) or simple correlation analysis are difficult to characterize this dynamic and non-stationary correlation mapping.

[0004] (2) High-dimensional coupling and noise interference of data. There are complex coupling relationships among microseismic monitoring parameters. For example, high-frequency low-energy events are often related to surface spalling, while low-frequency high-energy events may be related to deep structural slippage. Moreover, field monitoring data are easily affected by electromechanical interference, signal attenuation, sensor drift, etc., resulting in the coexistence of "false anomalies" and "missing true values". Existing methods often perform independent analysis or simple linear combination of parameters, failing to effectively mine the deep damage precursor information in the multi-parameter co-evolution, and also failing to fully consider the correction effect of geological background on the correlation strength.

[0005] Therefore, there is an urgent need in this field for a method to predict rock mass fracture indicators that can deeply integrate geological knowledge and has a strong ability to extract temporal features. Summary of the Invention

[0006] The technical problem to be solved by this invention is to provide a method for dynamic prediction of deep rock mass fracturing based on microseismic monitoring, with the aim of improving the accuracy of prediction of deep rock mass fracturing indicators.

[0007] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0008] A method for dynamic prediction of deep rock mass fracturing based on microseismic monitoring, the method comprising:

[0009] Step 1: Collect multi-source microseismic monitoring parameters, integrate them according to time series, and construct a multi-dimensional feature input sequence;

[0010] Step 2: Adaptively determine the time window and step size based on the rock mass fracture grade, and divide the multi-dimensional feature input sequence into multi-scale sub-sequences;

[0011] Step 3: Extract the key time delay feature matrix of rock mass fracturing from the multi-scale subsequence;

[0012] Step 4: Based on the key time delay feature matrix, output the predicted value of future rock mass fracturing index through the prediction model containing the Transformer encoder.

[0013] Furthermore, the multi-source microseismic monitoring parameters include: frequency, energy, focal depth, and magnitude of the microseismic event; the rock mass fracturing indicators include fracture density, fracturing zone extent, and rock mass risk status.

[0014] Further, step 2 includes: calculating the BQ value based on the uniaxial saturated compressive strength of the rock and the rock mass integrity coefficient to determine the rock mass fragmentation level; determining the time window and step size according to the rock mass fragmentation level; and dividing the multidimensional feature input sequence into multi-scale sub-sequences within the time window according to the set step size. Specifically, when BQ > 550, the time window is 1h to 24h and the step size is 2h; when 250 ≤ BQ ≤ 550, the time window is 1h to 12h and the step size is 1h; and when BQ < 250, the time window is 0.5h to 6h and the step size is 0.5h.

[0015] Furthermore, step 3 includes:

[0016] Step 31: Calculate the Pearson correlation coefficient between each subsequence and the rock mass fracturing index;

[0017] Step 32: Correct the Pearson correlation coefficient based on the set geological correction coefficient to obtain the corrected correlation degree;

[0018] Step 33: Calculate the mutual information value of any two microseismic monitoring parameters. Record the microseismic monitoring parameters with mutual information values ​​higher than the set value as strongly correlated parameter pairs. For strongly correlated parameter pairs, adjust the modified correlation degree based on the mutual information value and calculate the final correlation degree between the subsequence and the rock mass fracturing index. For microseismic monitoring parameters without strongly correlated parameters and with modified correlation degrees higher than the set value, adjust the modified correlation degree based on the set anomaly suppression factor and calculate the final correlation degree between the subsequence and the rock mass fracturing index.

[0019] Step 34: Assign weights to subsequences with a final correlation higher than a set threshold based on the Softmax function, and generate a weighted time delay feature matrix.

[0020] Furthermore, the calculation method for the geological correction factor mentioned in step 32 is as follows:

[0021] ;

[0022] in, This is a geological correction factor. For the maximum horizontal principal stress, For the minimum horizontal principal stress, For elastic modulus, Poisson's ratio, and It is a constant.

[0023] Furthermore, the method for calculating the corrected correlation degree in step 32 is as follows:

[0024] ;

[0025] in, For the first Microseismic monitoring parameters, For the first Rock mass fracturing index The Pearson correlation coefficient between the subsequence and the rock mass fracturing index. To correct the correlation.

[0026] Furthermore, step 33 includes:

[0027] Construct a mutual information matrix of microseismic monitoring parameters, and calculate the mutual information of the microseismic monitoring parameter pairs in the mutual information matrix:

[0028] ;

[0029] in, Microseismic monitoring parameters and mutual information value, and For its specific value, and For marginal probability density, For joint probability density;

[0030] For strongly correlated parameter pairs with mutual information values ​​greater than a set threshold, the final correlation between the subsequence and the rock mass fracturing index is calculated based on the mutual information value and the modified correlation degree.

[0031] ;

[0032] in, For the final correlation, This is the sum of mutual information of other microseismic monitoring parameters. and These are microseismic monitoring parameters. and Rock mass fracture index The corrected correlation degree;

[0033] For microseismic monitoring parameters with no strong correlation and whose correction correlation is higher than the set value, the anomaly suppression factor is used as the basis. Calculate the final correlation degree: .

[0034] Furthermore, step 34 uses the Softmax function to evaluate all final correlations higher than a set threshold. Assign weights to each subsequence, and assign weights to the t-th subsequence. for:

[0035] ;

[0036] in, For subsequence index, For subsequence The weighting of the allocation.

[0037] Furthermore, the method also includes: performing microseismic event classification and early warning based on the fracture density, fracture zone range, and rock mass risk status predicted by the trained prediction model.

[0038] Furthermore, the method also includes:

[0039] The trained prediction model serves as the teacher model, and the lightweight MobileViT architecture serves as the student model, employing a hierarchical structure. The loss function constrains the differences in output feature maps between the student model and the teacher model at the local feature learning layer and the global feature learning layer, respectively, in order to preserve key temporal features during model compression;

[0040] The network layers of the student model are divided into critical layers that are sensitive to accuracy and non-critical layers that have high tolerance for accuracy. Different integer quantization bits are used for critical layers and non-critical layers. The quantization bits of non-critical layers are dynamically adjusted according to the standard deviation of the input microseismic monitoring data.

[0041] The lightweight student model is optimized for edge devices by replacing the matrix multiplication operator with the Winograd convolution operator.

[0042] The beneficial effects of this invention are:

[0043] (1) This invention divides the microseismic monitoring parameters into different subsequences by dynamically classifying the rock mass fracturing level and adaptively configuring different time window ranges and step sizes for each level, so that the prediction model can accurately match the time delay characteristics of the response of microseismic and rock mass fracturing index under different rock mass conditions during the learning process.

[0044] (2) In the process of screening key time delay feature matrices, this invention first corrects the correlation coefficient between each subsequence and rock mass fracturing index by setting a geological correction coefficient based on the geological characteristics of high geostress, and obtains the corrected correlation degree; then, for strongly correlated microseismic monitoring parameter pairs, the mutual information value is used to correct the corrected correlation degree; for microseismic monitoring parameters without strong correlation but with a corrected correlation degree higher than the set value, the set anomaly suppression factor is used to correct them, and the final correlation degree between the subsequence and rock mass fracturing index is obtained. Finally, the Softmax function is used to assign weights to the subsequences with a final correlation degree higher than the set threshold, and a weighted time delay feature matrix is ​​generated. This matrix is ​​used as the input data of the prediction model, so that the prediction model can accurately capture the nonlinear time delay correlation between microseismic monitoring parameters and rock mass fracturing index under complex geological conditions during the learning and prediction process, and finally achieve a significant improvement in the prediction accuracy of rock mass fracturing index. Attached Figure Description

[0045] Figure 1 This is a flowchart of a method for dynamic prediction of deep rock mass fracturing based on microseismic monitoring in a specific embodiment of the present invention. Detailed Implementation

[0046] The core of the technical solution adopted by this invention to solve the above-mentioned technical problems is: collecting multi-source microseismic monitoring data and integrating them according to time series to construct a multi-dimensional feature input sequence; extracting the key time delay feature matrix of rock mass fracturing from the multi-dimensional feature input sequence; inputting the key time delay feature matrix into a prediction model containing a Transformer encoder for deep temporal feature extraction, and mapping the output high-dimensional features through a fully connected network of the output layer to obtain the predicted value of rock mass fracturing index for future time.

[0047] like Figure 1 As shown, the method for dynamic prediction of deep rock mass fracturing based on microseismic monitoring according to the present invention includes the following steps:

[0048] Step S1: Data Acquisition

[0049] (1) Microseismic monitoring parameter acquisition

[0050] High-precision microseismic monitoring equipment was deployed in the deep rock mass monitoring area according to the principle of "uniform coverage + intensified coverage in key areas". Data on the frequency, energy, focal depth, and magnitude of microseismic events were continuously collected at 15-30 minute intervals to ensure that the collected microseismic monitoring data covered different stages of rock mass activity. Multi-source microseismic monitoring parameters were integrated according to time series to construct a multi-dimensional feature input sequence.

[0051] (2) Rock mass physical parameter collection

[0052] The physical parameters of the rock mass include: uniaxial compressive strength of the rock. and rock mass integrity coefficient The uniaxial compressive strength of the rock was obtained by core sampling through on-site drilling, processing the cores into standard specimens, and then testing the uniaxial compressive strength in the laboratory using a pressure testing machine. The rock mass integrity factor was determined using an acoustic wave testing method, where acoustic wave transmitting and receiving devices were placed within the rock mass to measure the elastic longitudinal wave velocity. Simultaneously, elastic longitudinal wave velocity (pr) tests were conducted on rock samples taken from the same location, based on the formula... The rock mass integrity coefficient is calculated, with a value between 0 and 1. A value closer to 1 indicates better rock mass integrity, while a value further away indicates a higher degree of fragmentation. Alternatively, if elastic wave detection values ​​are unavailable, the number of volumetric joints in the rock mass can be statistically analyzed. Refer to relevant experience tables to determine the corresponding value.

[0053] The basic quality index BQ value of the rock mass is determined according to the People's Republic of China National Standard "Classification Standard for Engineering Rock Mass" (GB / T50218-2014), and the calculation formula is as follows: Calculations must adhere to specific conditions: when hour, Values ;when hour, Values The calculated BQ value can be used to preliminarily assess the basic quality of the rock mass and reflect its degree of fragmentation. For example, when the BQ value is greater than 550, the rock mass is hard and intact with a low degree of fragmentation; when the BQ value is less than 250, the rock mass is relatively soft and fragmented with a high degree of fragmentation.

[0054] (3) Geological parameter collection

[0055] Collect geostress and lithological data for the monitoring area, including the maximum horizontal principal stress. and minimum horizontal principal stress The lithological parameters include Poisson's ratio. and elastic modulus It is used for optimizing the geological adaptability of the model.

[0056] (4) Collection of rock mass fracturing index

[0057] Aligning with the acquisition time of microseismic monitoring data, rock mass fracturing indicators are obtained, including fracture density, fracturing zone extent, and rock mass risk status.

[0058] Step S2: Preprocess the data collected in step S1

[0059] Data cleaning: Data cleaning is a crucial step in ensuring data quality, aiming to remove interfering data, repair missing values, and improve data integrity and usability. First, it's necessary to identify and remove data points with obvious errors or abnormal fluctuations. These outliers are usually caused by factors such as monitoring equipment malfunctions or signal transmission interference, such as sudden changes in values ​​due to momentary sensor failure or extreme data exceeding reasonable ranges. For missing data, appropriate interpolation methods are selected based on the data characteristics for imputation, specifically including:

[0060] a) Linear interpolation: This method is suitable for situations where the data trend is relatively stable. It estimates missing values ​​by connecting adjacent known data points and using linear equations. This method is simple and efficient and can better maintain the overall trend of the data.

[0061] b) K-Nearest Neighbor Interpolation (KNN): By finding the K known data points that are most similar to the missing data points in terms of features, the missing value is calculated by weighting the data values ​​of these nearest neighbors. It is suitable for scenarios where the data distribution has local similarity and can effectively reflect the local features of the data.

[0062] Data standardization: Since actual collected data often have different dimensions and significantly varying numerical ranges, directly using the raw data for analysis may lead to model bias towards features with larger numerical values, affecting the accuracy and reliability of the analysis results. Therefore, Z-score standardization is necessary to perform dimensionless processing of the data, unifying all data to the same scale. Standardized data has a mean of 0 and a standard deviation of 1, eliminating the influence of dimensions and ensuring that different feature data are given equal importance. This facilitates subsequent data comparison, feature extraction, and model training, improving the convergence speed and prediction accuracy of the analysis model.

[0063] Step S3: Extract multi-scale subsequences from microseismic monitoring parameters

[0064] Based on the mapping rule between rock mass fracturing grade and time window parameters, the window range and step size under different fracturing grades are defined to ensure that the window matches the rock mass response characteristics. Multiple subsequences are extracted from the multidimensional feature sequence, specifically including:

[0065] When the BQ value is greater than 550, it indicates that the rock mass is hard and intact, with strong lag in rock mass fracturing. It is necessary to cover long-period delays (such as significant fracturing triggered by microseismic energy accumulation after 24 hours) and divide the time window of the multidimensional feature input sequence. and step length The values ​​are 1h~24h and 2h respectively; when 250≤BQ≤550, it indicates moderately fractured rock. The fracture response has both hysteresis and timeliness, balancing the comprehensiveness of features and computational efficiency, and dividing the multidimensional feature input sequence into time windows. and step length The time windows are 1h~12h and 1h, respectively. When BQ < 250, it indicates that the rock mass is relatively soft and fractured, with a rapid fracture response. It is necessary to accurately capture short-period sudden signals (such as local fracture triggered by high-frequency micro-vibrations within 0.5h) and divide the multi-dimensional feature input sequence into time windows. and step length The durations are 0.5h to 6h and 0.5h, respectively.

[0066] According to the set step size Within a time window, a sliding truncation of the microseismic monitoring parameter time series is performed to generate a set of subsequences. This is based on microseismic energy. For example, according to the set step size Obtain delay time ( The corresponding subsequence is ,in, As the starting index, Original sequence length, The values ​​are positive integers, and the final result is a time delay sequence matrix consisting of "microseismic monitoring parameters × window steps × subsequence length".

[0067] Compared to traditional fixed windows, in the hard and intact rock scenario, subsequence division based on rock mass fracture level improves the capture rate of long-period delayed rock mass fracture features by 40%, and reduces the redundancy of short-period signals by 60% compared to the soft fractured rock scenario, providing accurate subsequence data for subsequent correlation analysis.

[0068] Step S4: Filter out the key time delay feature matrix from the subsequences

[0069] The purpose of this step is to accurately capture the nonlinear time-delay correlation between microseismic monitoring parameters and rock mass fracturing indicators under complex geological conditions, and use it as input to the prediction model. The specific process is as follows:

[0070] Step S41: Calculate the Pearson correlation coefficient between each subsequence and the rock mass fracturing index.

[0071] The basic correlation between each microseismic monitoring parameter subsequence (such as the "high-frequency subsequence delayed by 2 hours" and the "low-energy subsequence delayed by 4 hours") and the fracture index was calculated based on the Pearson correlation coefficient formula. The original correlation strength between "single microseismic monitoring parameter and fracture index" is obtained. The calculation formula is as follows:

[0072]

[0073] in , The first Microseismic monitoring parameters and rock mass fracturing indicators and These are microseismic monitoring parameters. and rock mass fracturing index The A specific value, and These represent the mean values ​​of microseismic monitoring parameters and rock mass fracturing indices in the subsequences, respectively. The number of elements in the subsequence.

[0074] Step S42: Calculate the corrected correlation degree

[0075] The geological correction coefficient Based on the maximum horizontal principal stress collected in step S1 Minimum horizontal principal stress Poisson's ratio and elastic modulus To determine this, the calculation formula is as follows:

[0076] ;

[0077] in, and It is a constant.

[0078] Corrected correlation based on geological correction coefficient for:

[0079] ;

[0080] in, For the first Microseismic monitoring parameters, For the first Individual rock mass fracturing indicators.

[0081] Step S43: Calculate the final correlation degree

[0082] Because microseismic activity exhibits multi-parameter coupling characteristics—for example, "high frequency, low energy" corresponds to localized minor rupture, while "deep source, high energy" corresponds to widespread rupture—traditional single-parameter correlation methods are prone to producing "spurious correlations" (high correlation due to equipment interference but without practical significance) or "correlation omissions" (low correlation of a single parameter but strong correlation after multi-parameter synergy). By identifying parameter coupling relationships through mutual information analysis and integrating synergistic correlation characteristics, we can achieve "removing false correlations and strengthening aggregation," specifically as follows:

[0083] (1) Construct a mutual information matrix of microseismic monitoring parameters, and calculate the mutual information to quantify the dependency between two microseismic monitoring parameters (value ≥ 0, the larger the value, the stronger the dependency). The calculation formula is as follows:

[0084] ;

[0085] in, Microseismic monitoring parameters and mutual information value, and For its specific value, and These are the marginal probability densities, Let be the joint probability density.

[0086] filter Strongly correlated microseismic monitoring parameter pairs, among which For the set mutual information threshold, an example value in this embodiment is taken as follows: .

[0087] (2) For strongly correlated parameter pairs with mutual information values ​​greater than a set threshold, calculate the final correlation between the subsequence and the rock mass fracturing index based on the mutual information value and the corrected correlation degree;

[0088] ;

[0089] in, For the final correlation, This is the sum of mutual information of other microseismic monitoring parameters. and These are microseismic monitoring parameters. and Rock mass fracture index The corrected correlation.

[0090] (3) Introduce abnormal inhibitory factors ( For microseismic monitoring parameters with no strong correlation and whose correction correlation is higher than the set value, the final correlation is calculated based on the set anomaly suppression factor. This breaks the false correlation with the brokenness indicator.

[0091] After the key delay feature screening process (1)-(3) above, the accuracy rate is improved by 35% compared with single parameter correlation, and the false correlation signal elimination rate is improved by 50%, providing a reliable correlation basis for geological correction.

[0092] Step S43: Weight Normalization

[0093] Based on the Softmax function, the final correlation degree is higher than a set threshold. The subsequences with a weight of 0.5 (the optimal value for engineering verification) are assigned weights to generate a weighted time delay feature matrix.

[0094] Softmax weighting is used to determine the final correlation. As input, the weights of each subsequence are calculated using the softmax function, with the following formula: ,in, The weight is the number of subsequences with a final correlation higher than a set threshold. And the sum is 1, strongly correlated subsequences receive higher weights (e.g. The subsequences are 0.8, 0.6, and 0.5, with weights of approximately 0.48, 0.32, and 0.20, respectively.

[0095] Step S5: Train the prediction model containing the Transformer encoder based on the key time delay feature matrix.

[0096] The Transformer encoder layer employs a multi-layer stacked Transformer encoder architecture. Each encoder unit integrates a multi-head self-attention mechanism, a feedforward neural network (FFNN), residual connections, and layer normalization. The self-attention mechanism deeply mines the nonlinear coupling relationships between multiple parameters by calculating the attention weight matrix between sequence elements. The FFNN module performs a linear transformation on the multi-head attention output (512→2048→512) and enhances the nonlinear feature extraction capability through the ReLU activation function. The output of each sub-layer (multi-head attention, FFNN) performs a "residual connection + layer normalization" operation to avoid gradient vanishing and accelerate training convergence. The encoder layer is stacked with 2-3 layers. The first layer focuses on local time delay features, and the second layer captures global correlations across delays, enhancing the model's ability to represent complex temporal relationships.

[0097] A weighted time-delay feature matrix is ​​used as input data, and the corresponding rock mass fracturing index is used as the ground truth to build the model training dataset. The data is divided into a 70% training set, 15% validation set, and 15% test set, using a stratified sampling strategy to ensure similar data distribution across sets and to guarantee that the dataset covers data samples from different geological conditions (e.g., burial depth, lithology, and geostress state) and rock mass fracturing stages (intact, slightly fracturing, severely fracturing). Before partitioning, the training dataset needs to be standardized to unify the units and avoid the impact of data scale differences on model training performance. After partitioning, the training, validation, and test sets are stored separately for later retrieval.

[0098] Loss Function and Optimizer: The Mean Squared Error (MSE) loss function is used to measure the difference between the predicted and actual values. MSE, by calculating the mean of the squared errors, can amplify the impact of larger errors, prompting the model to pay more attention to samples with large prediction biases. The Adam optimizer is selected to tune the model parameters. This optimizer combines the advantages of Adagrad in handling sparse gradients and RMSProp in handling non-stationary objectives, with an initial learning rate set to 0.001. During training, the Adam optimizer dynamically adjusts the learning rate of each parameter based on the first and second moment estimates of the gradient, adapting to the update needs of different parameters.

[0099] Step S6: Based on the fracture density and fracture zone range predicted by the trained prediction model, perform microseismic event classification and early warning.

[0100] Based on the risk level of rock mass fracturing indicators, a three-level early warning system is established, with corresponding differentiated response mechanisms to avoid "over-warning (wasting resources)" or "under-warning (missing risks)," ensuring precise and efficient risk management. The early warning level classification criteria (based on predicted fracturing indicators for the next 24 hours) are as follows:

[0101] Blue alert: Fractured density D < 0.1, fractured zone range S < 5, rock mass risk status is stable (no risk), the system automatically generates a "Routine Monitoring Report" which includes microseismic data trends and fracture index prediction values, and pushes it to the construction management terminal. No adjustment to the construction plan is required, only monitoring at the original frequency is needed.

[0102] Yellow Alert: Crack density 0.1≤D<0.3, fractured zone range 5≤S<20, slight risk (local fracture possible). 1. Automatically trigger the generation of "local reinforcement support plan": based on the predicted fractured zone location (combined with microseismic source depth positioning), recommend support parameters (such as increasing anchor density to 1 bolt / m², increasing shotcrete thickness by 50mm, and anchor cable length ≥4m); 2. Push the plan to the construction team terminal, indicating the implementation scope of the plan (3m around the fractured zone), the impact on the construction period (increase of 0.5-1 hour / cycle), and the material usage (anchor bolts increased by 20%); 3. Update the predicted value every 2 hours to track risk changes.

[0103] Red Alert: Crack density D≥0.3, fractured zone range S≥5, severe risk (large-scale fracture possible). 1. Immediately trigger audible and visual alarms (on-site equipment horn + flashing red light), and simultaneously push emergency warning information (including risk area location and estimated fracture time) to the management personnel's mobile APP; 2. Automatically generate a "Personnel Evacuation + Work Stoppage Plan": specifying the evacuation radius (≥50m, calculated based on the fractured zone range and ground stress distribution), evacuation route (avoiding high-stress areas), and work stoppage duration (≥4 hours, until the risk decreases); 3. Link with the engineering monitoring platform: retrieve camera footage from the risk area in real time to confirm personnel evacuation, and update the predicted value every 30 minutes during the work stoppage until the warning level drops to yellow or below.

[0104] Preferably, this invention also performs accuracy-balanced optimization on the trained prediction model. Through "hierarchical knowledge distillation, dynamic quantization adaptation, and edge device hardware optimization," it significantly reduces model size and improves inference speed while strictly controlling accuracy loss, ensuring that the model can run stably on edge devices. The specific implementation is as follows:

[0105] a) Layered knowledge distillation and compression

[0106] Traditional knowledge distillation only transmits the probability distribution of the model's final output, easily losing key features from intermediate layers (such as local / global correlation features of the Transformer encoder), leading to a significant decrease in the accuracy of the student model. This mechanism, through "hierarchical feature transmission," transmits the features of each layer of the teacher model to the corresponding layers of the student model, ensuring that both local delayed features and global correlation features are accurately learned, maximizing accuracy preservation while compressing parameters.

[0107] Teacher model setting: The original Transformer model after training is used as the teacher model. The output feature maps of each encoder layer are saved. The first layer is a local feature map, focusing on the correlation of local time delay in microseismic events. The second and third layers are global feature maps, focusing on the global parameter coupling across the delay step.

[0108] Student Model Construction: The student model is built based on the MobileViT architecture, with corresponding local feature learning layers and global feature learning layers. MobileViT's lightweight backbone network (including attention and convolutional fusion structure) can maintain feature extraction capabilities with a small number of parameters.

[0109] Layer loss constraint: adopt The loss function constrains the difference between the corresponding layer feature maps of the student model and the teacher model, and its formula is as follows: ,in For the feature map of the i-th layer of the student model, For the feature map of the i-th layer of the teacher model, The number of elements in the feature map is denoted by a constraint loss threshold of <0.03. Hierarchical distillation weights are also set: local feature distillation weight 0.4 (to ensure the propagation of locally delayed features), global feature distillation weight 0.6 (to prioritize the propagation of globally related features), and the total distillation loss is the weighted sum of the losses from the two layers.

[0110] Compression effect: After stratified distillation, the student model parameters were compressed from 14.8MB to 2.8MB, the inference speed was improved by more than 4 times compared with the original model, and the increase in RMSE of the test set was strictly controlled within 1%, achieving a significant reduction in data processing volume.

[0111] b) Dynamic quantization adaptation

[0112] Traditional fixed-bit quantization (such as uniform 16-bit) does not differentiate the impact of each model layer on accuracy. Critical layers (such as multi-head attention layers) are sensitive to accuracy, and low bit quantization can lead to large deviations in attention weight calculations. Non-critical layers (such as feedforward neural network layers) have a high tolerance for accuracy and can use even lower bit quantization to further compress the size. This mechanism, through "layered differentiated quantization + dynamic threshold adjustment," maximizes size compression while avoiding the loss of accuracy in critical layers.

[0113] Model layer classification: The student model layer is divided into "critical layer" and "non-critical layer". The feedforward neural network (FFNN) layer is a non-critical layer (only performs linear transformations and has high accuracy tolerance).

[0114] Differentiated quantization strategy: Key layer (multi-head attention): Use 16-bit integer quantization (int16) to preserve the precision of weight calculation and avoid distortion of attention correlation;

[0115] Non-critical layer (FFNN): 8-bit integer quantization is used by default, which greatly compresses parameters while keeping the accuracy loss controllable (<0.5%).

[0116] Dynamic threshold adjustment: Introducing the standard deviation of input data The formula used as the basis for dynamically adjusting the quantization threshold is: ,when When data fluctuations are small and feature stability is high, non-critical layers can be further reduced to 4-bit quantization (int4), at which point the parameter volume is reduced by another 50%, and the RMSE increase is <0.5%; when When data fluctuates significantly, non-critical layers maintain 8-bit quantization to ensure that features are not lost.

[0117] Quantization tools and results: Quantization was achieved using PyTorch's torch.quantization tool. After quantization, the model parameters were further compressed from 2.8MB to less than 2.2MB, which is more than 85% of the original model parameters, and the total accuracy loss was less than 2% (including distillation and quantization).

[0118] c) Edge device hardware adaptation

[0119] Edge devices used in deep engineering typically suffer from low CPU clock speeds (mostly 1.5-2GHz), small memory (2-4GB), and susceptibility to electromagnetic interference, high temperatures, and humidity. Traditional lightweight design focuses only on parameter compression, neglecting hardware characteristics and environmental interference, leading to model stuttering or interruptions. This mechanism adapts to edge devices in terms of computation, memory, and stability through "operator optimization, memory reuse, and anti-interference encapsulation," ensuring reliable model operation. The specific implementation is as follows:

[0120] Computation operator optimization: The traditional matrix multiplication operator in Transformer is replaced with the Winograd convolution operator. The Winograd operator transforms matrix multiplication into more efficient convolution calculation through mathematical transformation, which can reduce the amount of computation by 30% (e.g., the number of calculations for 512×512 matrix multiplication is reduced from about 2.6 million to 1.8 million), and can significantly improve inference speed on low-frequency CPUs.

[0121] Memory reuse optimization: A "dynamic memory allocation" strategy is adopted to allocate memory blocks in a loop for intermediate feature maps of the model (such as quantized feature matrices). The same memory block is used to store intermediate features of different layers in sequence (such as storing the output of the attention layer first, and then covering the output of the FFNN layer). There is no need to allocate memory separately for each layer, reducing memory usage by 45% (from 1.2GB to 0.66GB), which is suitable for the small memory characteristics of edge devices.

[0122] Anti-interference encapsulation: An electromagnetic interference detection module is embedded in the model inference code to monitor the fluctuation of the device power supply voltage in real time. When a voltage fluctuation > ±10% is detected, the cached feature map recalculation mechanism is immediately activated: the key feature map before the fluctuation (such as the quantized feature) is cached locally. If the fluctuation causes the current calculation to be interrupted, the cached feature is directly called to recalculate, avoiding inference failure. At the same time, the model file is encapsulated in a moisture-resistant and corrosion-resistant code (such as adding data check bits to prevent file storage errors).

[0123] Adaptation effect: The optimized model can run stably on edge devices with a main frequency of 1.5GHz and 2GB of memory, with a single sample inference time of ≤0.06 seconds. It can run continuously for 72 hours without interruption in harsh environments with temperatures of -20~60℃, relative humidity of 85%, and electromagnetic interference intensity of ≤10V / m, meeting the real-time and stability requirements of deep engineering.

[0124] Matrix multiplication is replaced with the Winograd convolution operator; dynamic memory reuse is adopted; electromagnetic interference detection is added, and cached feature map recalculation is enabled when voltage fluctuation is > ±10% to ensure stable operation in harsh environments.

[0125] As a preferred approach, considering that deep engineering geological conditions dynamically change with the progress of construction (such as stress release and lithological stratification changes), fixed models are prone to increasing prediction bias. This invention also employs a closed loop of "post-construction verification data - bias assessment - incremental training" to adaptively iterate the trained prediction model, as specifically implemented below;

[0126] Verification of data acquisition: Within 72 hours after construction, the actual fracture density D was obtained using borehole imaging technology. 实际 Ground-penetrating radar acquires the actual fractured zone range S 实际 This creates a validation dataset that pairs predicted values ​​with actual values.

[0127] Deviation rate calculation and evaluation: The deviation rate is used to quantify the prediction accuracy. The formula is as follows: when the deviation rate is >10% (the engineering acceptable threshold), the model is determined to need iteration; if the deviation rate is ≤10%, only the validation data is saved to the database and iteration is not triggered.

[0128] Incremental training execution: The newly collected validation data is added to the original training set and validation set at a ratio of 7:3 (stratified sampling to ensure that the distribution of the new data is consistent with the original data).

[0129] Training strategy: The "freeze-update" mode is adopted, freezing the encoder layer and lightweight module of the model (to avoid degradation of core feature extraction capabilities), and only updating the parameters of the fully connected layer of the output layer (to reduce the amount of computation).

[0130] The iteration is completed within 1 hour in a GPU environment (full training takes 20-30 hours). After each iteration, the RMSE of the new model is automatically calculated. The old model is overwritten only when the RMSE is less than that of the original model. The iteration log is saved: the "data source (e.g., a tunnel K12+300 section), deviation cause analysis (e.g., changes in ground stress leading to changes in the microseismic-fracture correlation), and parameter adjustment records (e.g., updated values ​​of the weights of the fully connected layer)" of each iteration are automatically recorded to form a traceable optimization archive, which is convenient for subsequent analysis of the model performance change pattern.

Claims

1. A deep rock mass breaking dynamic prediction method based on microseismic monitoring, characterized in that, The method comprises: Step 1: Collecting multi-source microseismic monitoring parameters, integrating in time sequence, and constructing multi-dimensional feature input sequence; Step 2: According to the rock mass fragmentation level, the time window and the step are adaptively determined, and the multi-dimensional feature input sequence is divided into multi-scale subsequences; Step 3: Extracting the rock mass fragmentation key time delay feature matrix from the multi-scale subsequence; specifically comprising: Step 31: Calculate the Pearson correlation coefficient of each subsequence and the rock mass fragmentation index; Step 32: Based on the set geological correction coefficient, the Pearson correlation coefficient is corrected to obtain the corrected correlation degree; Step 33: Calculate the mutual information value of any two microseismic monitoring parameters, and record the microseismic monitoring parameters with mutual information value higher than the set value as the strong correlation parameter pair; for the strong correlation parameter pair, the final correlation degree of the subsequence and the rock mass fragmentation index is calculated based on the mutual information value and the corrected correlation degree; for the microseismic monitoring parameters without strong correlation parameters and with corrected correlation degree higher than the set value, the final correlation degree of the subsequence and the rock mass fragmentation index is calculated based on the set abnormality suppression factor; Step 34: Based on the Softmax function, the weight of the subsequence with final correlation degree higher than the set threshold is allocated to generate a weighted time delay feature matrix; Step 4: Based on the key time delay feature matrix, a prediction model containing a Transformer encoder is used to output the prediction value of the future rock mass fragmentation index.

2. The deep rock mass breakage dynamic prediction method based on microseismic monitoring according to claim 1, characterized in that, The multi-source microseismic monitoring parameters include the frequency, energy, focal depth and magnitude of the microseismic event; and the rock mass fragmentation index includes the fracture density, fragmentation zone range and rock mass risk state.

3. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 1, characterized in that, Step 2 comprises: determining the rock mass fragmentation level based on the uniaxial saturated compressive strength of rock and the rock mass integrity coefficient, determining the time window and the step according to the rock mass fragmentation level, and dividing the multi-dimensional feature input sequence into multi-scale subsequences according to the set step within the time window, wherein when BQ>550, the time window is 1h~24h and the step is 2h; when 250≤BQ≤550, the time window is 1h~12h and the step is 1h; and when BQ<250, the time window is 0.5h~6h and the step is 0.5h.

4. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 1, characterized in that, The calculation method of the geological correction coefficient in step 32 is: ； wherein, is a geological correction factor, is the maximum horizontal principal stress, is the minimum horizontal principal stress, is the elastic modulus, is the Poisson's ratio, and is a constant.

5. The deep rock mass breakage dynamic prediction method based on microseismic monitoring of claim 4, characterized in that, The calculation method of the corrected correlation degree in step 32 is: ； wherein, is the first microseismic monitoring parameter, is the second microseismic monitoring parameter, is the third microseismic monitoring parameter, is the fourth microseismic monitoring parameter, is the Pearson correlation coefficient of the subsequence and the rock mass fragmentation index, is the corrected correlation degree.

6. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 5, characterized in that, Step 33 comprises: Constructing a microseismic monitoring parameter mutual information matrix, and calculating the mutual information of the microseismic monitoring parameter pair in the microseismic monitoring parameter mutual information matrix: ； wherein, is a microseismic monitoring parameter and is a mutual information value, and are specific values thereof, and are marginal probability densities, is a joint probability density; For the strong correlation parameter pair with mutual information value greater than the set threshold, the final correlation degree of the subsequence and the rock mass fragmentation index is calculated based on the mutual information value and the corrected correlation degree; ； wherein, is the final correlation degree, is the sub-sequence sequence number index, is the sum of other microseismic monitoring parameter mutual information, and are microseismic monitoring parameters and are the modified correlation degrees of the microseismic monitoring parameters and the rock mass fragmentation index. For microseismic monitoring parameters without strong correlation parameters and with a correction correlation degree higher than a set value, based on the set anomaly suppression factor Calculate the final correlation degree: .

7. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 6, characterized in that, Step 34 uses a Softmax function to assign weights to all sub-sequences whose final correlation is above a set threshold, the weight for the tth sub-sequence is: where N is the number of sub-sequences above the threshold. where N is the number of sub-sequences above the threshold. ； wherein, is a sub-sequence index, is a sub-sequence allocation weight.

8. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 2, characterized in that, The method further comprises: based on the prediction output of the fracture density, the fragmentation zone range and the rock mass risk state of the trained prediction model, a microseismic event grading early warning is performed.

9. The microseismic monitoring based deep rock mass breakage dynamic prediction method according to claim 1, characterized in that, The method further comprises: The trained prediction model is used as a teacher model, a lightweight MobileViT architecture is used as a student model, and a hierarchical Loss functions are used to constrain the output feature maps of the student model and the teacher model in the local feature learning layer and the global feature learning layer, respectively, so as to retain key timing features in model compression. The network layers of the student model are divided into key layers sensitive to accuracy and non-key layers with high accuracy tolerance, different integer quantization bit numbers are used for the key layers and the non-key layers, and the quantization bit number of the non-key layers is dynamically adjusted according to the standard deviation of the input microseismic monitoring data; The lightweight child model is adapted and optimized for edge devices, and the matrix multiplication operator is replaced by a Winograd convolution operator.

Images