A tunnel rock burst grade prediction method based on microseismic monitoring and finite element inversion

Abstract

The present application relates to the technical field of tunnel engineering and rock burst disaster prediction, and particularly relates to a tunnel rock burst grade prediction method based on microseismic monitoring and finite element inversion, comprising: obtaining tunnel geological features to construct a geological feature vector, collecting microseismic time series data at each monitoring position, including event rate, energy release rate and spatial aggregation degree. Initial section division is carried out with geological mutation points as boundaries, and the final monitoring section is obtained by adjusting the microseismic time series difference. The microseismic time series data is classified into the corresponding section, spliced with the geological features to form a composite feature vector, input into a neural network prediction model, and the rock burst grade probability of each section is output. The model loss function adds finite element inversion stress mean square error constraint and Hoek-Brown criterion physical penalty constraint. Finally, the probability is projected to the corresponding position of the three-dimensional model by color mapping and displayed. The present application improves the accuracy of tunnel rock burst grade prediction by fusing multi-source data and introducing physical constraints.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering and rockburst disaster prediction technology, and in particular to a method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion. Background Technology

[0002] Predicting the rockburst severity in tunnels is a key challenge for safety control during deep-buried tunnel construction. Existing prediction methods mainly fall into two categories: one relies on geological sketches or single microseismic indices for empirical judgment, failing to fully utilize multi-source heterogeneous data and struggling to capture the time-varying evolution characteristics of surrounding rock fracturing during rockburst incubation; the other, while incorporating data-driven models such as neural networks, is prone to overfitting under small sample sizes and label imbalance conditions, and its prediction results lack physical and mechanical basis, often violating the actual stress state and strength criteria of the surrounding rock, leading to a high misjudgment rate. Furthermore, traditional methods are insufficient in characterizing the correlation between the energy storage capacity of the surrounding rock and its real-time fracturing response, lack a mechanism for effectively integrating geological static properties with microseismic dynamic responses, and the model outputs are often uninterpretable, failing to reflect the actual distribution of stress in the surrounding rock. Summary of the Invention

[0003] This invention provides a method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion. It addresses the problems in existing rockburst prediction methods, such as the failure to effectively integrate geological static characteristics and microseismic dynamic monitoring data, and the lack of physical and mechanical constraints in purely data-driven models, which lead to insufficient prediction accuracy and reliability.

[0004] The objective of this invention can be achieved through the following technical solutions:

[0005] The first aspect of this invention is to provide a method for predicting the rockburst intensity of tunnels based on microseismic monitoring and finite element inversion, comprising: Geological feature data and microseismic time series data in the tunnel are acquired; a geological feature vector is constructed using the geological feature data; the microseismic time series data includes the event rate, energy release rate, and spatial clustering of each segment within a preset time window; the segment is the result of division along the tunnel axis; the spatial clustering is obtained by calculating the covariance matrix of the source coordinates of each monitored microseismic event and taking the square root of the sum of the three eigenvalues; First, the tunnel is divided into initial sections using geological feature data to obtain an initial section set. Then, based on the differences in microseismic time series data between adjacent initial sections in the initial section set, the division results of the initial sections are adjusted to obtain several adjusted monitoring sections. The event rate, energy release rate, and spatial clustering within a time window of each monitoring segment are used to form the microseismic feature vector of each monitoring segment at each time point. The geological feature vector and the microseismic feature vector of each monitoring segment at each time point are concatenated to form the composite feature vector of each monitoring segment at each time point. The composite feature vector is input into a neural network prediction model for prediction to obtain the rockburst level probability of each monitoring segment. The total loss function in the neural network prediction model incorporates the mean square error constraint term of the finite element inversion stress distribution and the physical penalty constraint term of the Hoek-Brown strength criterion. The probability of rockburst levels in all monitored sections is input into a 3D visualization map for display.

[0006] Furthermore, the specific process for acquiring the microseismic time series data includes: The event rate is calculated by counting the total number of micro-seismic events in each segment within a preset time window; the energy release rate is calculated by counting the sum of the micro-seismic energy of all micro-seismic events in each segment within a preset time window. The process of obtaining the spatial clustering degree is as follows: Obtain the source coordinates of all microseismic events within each segment; form a symmetric covariance matrix based on the covariances of the source coordinates of all microseismic events in the three coordinate directions; solve the characteristic equation of the symmetric covariance matrix to obtain the three eigenvalues. , , The sum of the three factors, then the square root, gives the spatial density. ; The covariance symmetric matrix is ​​specifically represented as follows:

[0007] In the formula, This represents the covariance of the x-coordinate itself. This represents the covariance between the x-coordinate and y-coordinate directions. This represents the covariance between the x-coordinate and the z-coordinate. This represents the covariance between the y-coordinate and the z-coordinate. This represents the covariance along the y-coordinate direction itself. This represents the covariance of the z-coordinate itself. in, , , .

[0008] Furthermore, the process of obtaining the initial segment set includes: Along the tunnel mileage, arrange the geological feature data in the order of each cycle of advance, and check whether there are abrupt changes in lithology, fault markers, and structural surface density. If any of the lithology, fault markers, and structural surface density changes, mark it and record the marked position as the initial boundary. An initial segment is formed between adjacent initial boundaries. After all initial boundaries are marked, output the set of initial segments.

[0009] Furthermore, the specific acquisition process of the aforementioned monitoring sections includes: The event rate, energy release rate, and spatial clustering degree within each initial segment are statistically analyzed. Based on these three parameters, the microseismic response difference between adjacent initial segments is calculated. When the microseismic response difference between adjacent initial segments is greater than or equal to a preset difference threshold, the microseismic index difference between the adjacent initial segments is determined to be significant. When the microseismic response difference between adjacent initial segments is less than the preset difference threshold, the microseismic index difference between the adjacent initial segments is determined to be insignificant. If the microseismic indices of adjacent initial segments differ significantly, the corresponding boundary is retained; if the microseismic indices of adjacent initial segments do not differ significantly, the corresponding boundary is deleted and the two initial segments are merged. The difference in microseismic response between adjacent initial sections is specifically expressed by the following formula:

[0010] In the formula, This represents the event rate of the first initial segment. This represents the event rate of the second initial segment. This represents the energy release rate of the first initial segment. This indicates the energy release rate of the second initial segment. This indicates the spatial clustering degree of the first initial segment. This indicates the spatial clustering degree of the second initial segment; This indicates the degree of difference in microseismic response between adjacent initial sections. Represents a linear normalization function. It is the absolute value symbol.

[0011] Furthermore, the neural network prediction model is an LSTM neural network model; The LSTM neural network model includes an LSTM layer, a fully connected layer, and an output layer; it also includes an independent single-neuron fully connected layer.

[0012] Furthermore, the total loss function consists of the cross-entropy loss function, the finite element stress calibration loss function, and the Hoek-Brown criterion constraint loss function; The total loss function is specifically expressed by the formula:

[0013] In the formula, This represents the cross-entropy loss value. This represents the finite element stress calibration loss value. This represents the constraint loss value of the Hoek-Brown criterion. This represents the first preset weighting coefficient. This represents the second preset weighting coefficient. This represents the total loss value; The finite element stress calibration loss function is specifically expressed as follows:

[0014] In the formula, Indicates the first Estimated surrounding rock stress for each sample Indicates the first The stress values ​​at the corresponding locations of each sample were obtained through finite element inversion calculation; Indicates the number of samples in the training batch; The Hoek-Brown criterion constraint loss function is specifically expressed as follows:

[0015] In the formula, Indicates the first Estimated surrounding rock stress for each sample represents the uniaxial compressive strength of the rock mass at the location corresponding to the kth sample, which is obtained by looking up the lithology table in the geological sketch or by calculating the Hoek-Brown parameter; N represents the number of samples in the training batch. This represents the function that takes the maximum value. Among them, a sample refers to the composite feature vector of a monitoring section at a single time point, whether a rockburst will occur in the monitoring section in the next day and the corresponding level label, and the stress value obtained by inversion calculation corresponding to the location of the monitoring section. The process of obtaining the estimated surrounding rock stress for each sample in the finite element stress calibration loss function and the Hoek-Brown criterion constraint loss function is as follows: the hidden state output by the LSTM layer in the LSTM neural network model is then passed through a single neuron fully connected layer to obtain the estimated surrounding rock stress for each sample.

[0016] Furthermore, the process of obtaining the stress value at each sample location in the finite element stress calibration loss function through finite element inversion calculation is as follows: Step 1: Obtain the tunnel's geometric parameters from the design drawings. These parameters include cross-sectional dimensions and burial depth. Based on these geometric parameters, obtain the geological parameters of each rock layer using geological sketches and borehole data. These geological parameters include the boundaries of each rock layer, fault locations, and dip angles. Input these geological parameters into finite element software to construct a three-dimensional geological model. Step 2: Obtain the physical and mechanical parameters of each rock layer through indoor rock mechanics tests; the physical and mechanical parameters include elastic modulus, Poisson's ratio, density, cohesion, and internal friction angle; assign the physical and mechanical parameters of each rock layer to the corresponding rock layer unit in the three-dimensional geological model; Step 3: Obtain the magnitude and direction of the principal stress at the tunnel burial depth by measuring the in-situ stress using the hydraulic fracturing method, and apply it as the initial stress field boundary condition to the outer boundary of the three-dimensional geological model. Step 4: Extract the stress components at the center of each monitoring section using a three-dimensional geological model. The stress components at the center of each monitoring section are then decomposed using the stress tensor eigenvalues ​​to obtain the maximum principal stress, which is used as the stress value calculated by finite element inversion at the center of each monitoring section. Step 5: Compare the stress value calculated in Step 4 with the actual stress data measured on site. If the deviation exceeds the allowable range, adjust the input parameters, return to Step 3 to reapply the boundary and solve again. Iterate repeatedly until the deviation between the calculated value and the measured value meets the requirements. At this point, the model is considered to represent the actual working condition. Output the stress value calculated by finite element inversion under this working condition for use in the loss function.

[0017] Furthermore, if the deviation in step 5 exceeds the allowable range, the input parameters are adjusted; the specific process includes: Step 51: The magnitude of structural stress, lateral compression coefficient, and elastic modulus reduction coefficient are the variables to be inverted, and a search range is set for each variable; Step 52: Using orthogonal experimental design or Latin hypercube sampling, generate multiple sets of parameter combinations within the search range set for each variable. Substitute each set of combinations into the finite element model for calculation to obtain the corresponding calculated stress value. Step 53: Construct a mapping model between the input parameter combination and the calculated stress values; Step 54: Within the search range set for each variable, an optimization algorithm is used to automatically search to minimize the error between the calculated stress value and the measured value; the objective function of the optimization algorithm is the error function.

[0018] In the formula, This represents the stress value obtained through finite element inversion calculation for each parameter combination. This represents the measured stress value. It is the absolute value symbol. Indicates stress value error; The optimization algorithm is either a genetic algorithm or a particle swarm optimization algorithm.

[0019] A second aspect of the present invention is to provide an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the aforementioned method for predicting tunnel rockburst levels based on microseismic monitoring and finite element inversion.

[0020] A third aspect of the present invention is to provide a computer-readable storage medium storing a computer program that, when executed by a processor, implements the method for predicting tunnel rockburst levels based on microseismic monitoring and finite element inversion.

[0021] Compared with the prior art, the beneficial effects of the present invention are: Geological feature data and microseismic time series data of the tunnel are acquired; a geological feature vector is constructed using the geological feature data; the microseismic time series data includes the event rate, energy release rate, and spatial clustering of each segment within a preset time window; the segments are the result of division along the tunnel axis; the spatial clustering is obtained by calculating the covariance matrix of the source coordinates of each monitored microseismic event and taking the square root of the sum of the three eigenvalues; by fusing geological and microseismic multi-source data, a feature basis that comprehensively reflects the static properties and dynamic response of the surrounding rock is provided for rockburst prediction; firstly, the tunnel is divided into initial segments using geological feature data to obtain an initial segment set; then, based on the differences in microseismic time series data between adjacent initial segments in the initial segment set, the initial segment division results are adjusted to obtain several adjusted monitoring segments; by defining boundaries through geological abrupt change points and verifying the adjustments with microseismic data, the consistency of the geomechanical properties within the segments and the uniformity of the microseismic response behavior are ensured, improving the targeting of subsequent predictions; The event rate, energy release rate, and spatial clustering of each monitoring segment within a time window are used to form a microseismic feature vector for each monitoring segment at each time point. The geological feature vector and the microseismic feature vector of each monitoring segment at each time point are concatenated to form a composite feature vector for each monitoring segment at each time point. This composite feature vector is then input into a neural network prediction model to obtain the rockburst probability for each monitoring segment. The total loss function of the neural network prediction model incorporates a mean square error constraint term for finite element inversion stress distribution and a physical penalty constraint term for the Hoek-Brown strength criterion. By introducing the mean square error constraint for finite element stress and the physical penalty constraint term for the Hoek-Brown criterion into the total loss function, the prediction results conform to mechanical laws, improving the prediction reliability under sparse data conditions. The rockburst probability of all monitoring segments is input into a 3D visualization map for display. The spatial distribution of rockburst risk is visualized through 3D model color mapping, facilitating intuitive location of high-risk segments by construction personnel. This approach solves the problem of insufficient accuracy and reliability in existing rockburst prediction methods, where static geological features and dynamic microseismic monitoring data are not effectively integrated, and purely data-driven models lack physical and mechanical constraints. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This invention provides a flowchart illustrating the steps of a method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion. Detailed Implementation

[0024] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0025] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0026] To address the problems existing in the background technology, a method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion was designed, which has important practical significance.

[0027] like Figure 1 As shown, the first aspect of this invention is to provide a method for predicting the rockburst severity of tunnels based on microseismic monitoring and finite element inversion, comprising the following steps: Step S1: Obtain geological feature data and microseismic time series data in the tunnel; construct a geological feature vector using the geological feature data; the microseismic time series data includes the event rate, energy release rate, and spatial clustering of each segment within a preset time window; the segment is the result of dividing along the tunnel axis; the spatial clustering is obtained by calculating the covariance matrix of the source coordinates of each monitored microseismic event and taking the square root of the sum of the three eigenvalues.

[0028] It should be noted that geological sketches record static features such as lithology, structural planes, and faults at different mileages, while microseismic signals reflect the dynamic changes in surrounding rock fracturing activity in real time. Both types of information characterize the conditions for rockburst formation and its evolution from different dimensions. To integrate this information for subsequent rockburst severity prediction, it is necessary to collect geological features and corresponding microseismic time-series data for each mileage segment as the foundational data source for model training and prediction.

[0029] Specifically, the process of acquiring geological feature data from the tunnel and constructing geological feature vectors is as follows: After each excavation cycle (each cycle advances 2-3m), the geological engineer goes to the working face to conduct a geological sketch, observes with the naked eye and measures with tools such as a compass and magnifying glass, and records the lithology, surrounding rock grade, structural surface density, groundwater status, and fault fracture zone conditions of that section; the lithology, surrounding rock grade, structural surface density, groundwater status, and fault fracture zone conditions of each section are combined into a 5-dimensional geological feature vector.

[0030] The numerical coding of lithology, surrounding rock grade, structural surface density, groundwater condition, and fault markers (fault fracture zone conditions) is as follows: Lithology: Assigned values ​​from table. Granite = 1, granodiorite = 2, slate interbedded with metamorphic sandstone = 3, diorite = 4; Rock mass classification: directly take the classification value. Class II = 2, Class III = 3, Class IV = 4, Class V = 5; Structural surface development degree: Values ​​are assigned based on joint density. <1 joint / m = 0.2, 1~3 joints / m = 0.5, 3~5 joints / m = 0.8, >5 joints / m = 1.0; Groundwater status: Dry = 0, Moist = 0.3, Drip = 0.6, Inflow = 1.0; Fault fracture zone: None = 0, Present = 1.

[0031] It should be noted that the values ​​of the five-dimensional geological feature vectors remain unchanged.

[0032] The specific process for obtaining microseismic time series data in the tunnel is as follows: Several monitoring locations are determined at predetermined, evenly spaced intervals. Accelerometers are installed on the tunnel wall at each location. These sensors convert the vibration signals generated by the surrounding rock fracturing into electrical signals. These signals are then converted into digital waveforms by a digital signal acquisition instrument at a sampling rate of 10,000–20,000 Hz. After continuous data acquisition, three parameters—event rate, energy release rate, and spatial clustering—are statistically analyzed for each segment within a 2-hour time window, forming multidimensional time-series data that evolves over time. The length of the time window is not specifically limited and can be set by the implementer according to specific circumstances.

[0033] In this embodiment, the preset spacing is 5m. However, the preset spacing is not specifically limited in this embodiment, and the implementer can determine it according to the specific situation.

[0034] The event rate is calculated by counting the total number of microseismic events in the section within a preset time window (2 hours). The energy release rate is calculated by counting the sum of the microseismic energy of all microseismic events in the section within a preset time window (2 hours); for each event, the original energy value at the source is calculated by back-calculating the amplitude received by multiple sensors.

[0035] The process for obtaining the spatial clustering degree of each segment is as follows: Obtain the source coordinates of all microseismic events within each segment; based on the covariance of the source coordinates of all microseismic events in the three coordinate directions, form a 3×3 covariance symmetric matrix; solve the characteristic equation of this matrix to obtain the three eigenvalues. , , The sum of the three factors, then the square root, gives the spatial density. Specifically, this can be expressed by the formula: The smaller the value, the closer the points are, and the more concentrated they are; the larger the value, the more dispersed the points are.

[0036] The covariance symmetric matrix is ​​specifically represented as follows: .

[0037] In the formula, This represents the covariance of the x-coordinate itself. This represents the covariance between the x-coordinate and y-coordinate directions. This represents the covariance between the x-coordinate and the z-coordinate. This represents the covariance between the y-coordinate and the z-coordinate. This represents the covariance along the y-coordinate direction itself. This represents the covariance of the z-coordinate itself.

[0038] in, , , The formula for calculating covariance is as follows (one is given, others are similar):

[0039] In the formula, This represents the value in the x-coordinate direction of the i-th earthquake source. This represents the value in the y-coordinate direction of the i-th seismic source. This represents the average value of the x-coordinate values ​​for all earthquake source coordinates. This represents the average value of the y-coordinate values ​​for all earthquake source coordinates. This indicates the number of coordinates of all earthquake sources.

[0040] Step S2: First, the tunnel is divided into initial sections using geological feature data to obtain an initial section set; then, based on the differences in microseismic time series data between adjacent initial sections in the initial section set, the division results of the initial sections are adjusted to obtain several adjusted monitoring sections.

[0041] It should be noted that geological characteristics (lithology, faults, structural planes) determine the mechanical properties and energy storage capacity of the surrounding rock, forming the material basis for rockbursts. These characteristics exhibit discrete and abrupt changes in space, making them suitable for defining hard boundaries as the initial framework. Microseismic data, on the other hand, reflects the real-time fracture response of the surrounding rock under excavation disturbance, verifying whether there are genuine differences in microseismic behavior on both sides of the boundary. Therefore, to obtain reasonable segmentation results that conform to geomechanical laws and are verified by actual microseismic responses, we first use geological characteristics for initial segmentation, then use microseismic data to verify the validity of the boundaries and perform merging or densification adjustments.

[0042] Specifically, the process of obtaining the initial segment set includes: Along the tunnel mileage, geological feature data are arranged in order of each cycle advance (2~3m), and the three key indicators (lithology number, fault mark, and structural surface density) are checked one by one to see if there is a sudden change. If the lithology number of a location changes, or the fault mark changes from 0 to 1, or the structural surface density changes across grades (marking is performed as long as one of the three changes), then the location is marked as the initial boundary. An initial segment is formed between adjacent boundaries. After all boundaries are marked, the set of initial segments is output.

[0043] Then, based on the differences in microseismic time series data between adjacent initial segments in the initial segment set, the initial segment division results are adjusted to obtain several adjusted monitoring segments; the specific process is as follows: The event rate, energy release rate, and spatial clustering degree within each initial segment are statistically analyzed. Based on the event rate, energy release rate, and spatial clustering degree of adjacent initial segments, the microseismic response difference between adjacent initial segments is calculated. When the microseismic response difference between adjacent initial segments is greater than or equal to a preset difference threshold, the microseismic index difference between adjacent initial segments is determined to be significant. When the microseismic response difference between adjacent initial segments is less than the preset difference threshold, the microseismic index difference between adjacent initial segments is determined to be insignificant.

[0044] In this embodiment, the preset difference threshold is set to 0.5. However, the preset difference threshold is not specifically limited in this embodiment, and the implementer can determine it according to the specific circumstances.

[0045] If the microseismic indices of adjacent initial sections differ significantly (indicating that the initial division does indeed have a distinguishing effect on microseismic response), the corresponding boundary is retained; if the microseismic indices of adjacent initial sections do not differ significantly (indicating that the microseismic behavior of the surrounding rock on both sides of the boundary is not substantially different), the corresponding boundary is deleted and the two initial sections are merged.

[0046] The microseismic response difference between adjacent initial sections is calculated based on three parameters: event rate, energy release rate, and spatial clustering. The specific formula for expressing the microseismic response difference between adjacent initial sections is as follows:

[0047] In the formula, This represents the event rate of the first initial segment. This represents the event rate of the second initial segment. This represents the energy release rate of the first initial segment. This indicates the energy release rate of the second initial segment. This indicates the spatial clustering degree of the first initial segment. This indicates the spatial clustering degree of the second initial segment; This indicates the degree of difference in microseismic response between adjacent initial sections. Represents a linear normalization function. It is the absolute value symbol.

[0048] in, This indicates the degree of difference in event rates between two adjacent initial segments. This indicates the degree of difference in energy release rate between two adjacent initial sections. It indicates the degree of difference in spatial clustering between two adjacent initial segments.

[0049] At this point, several monitoring sections have been adjusted.

[0050] Step S3: The event rate, energy release rate, and spatial clustering within a time window of each monitoring segment are used to form the microseismic feature vector of each monitoring segment at each time point; the geological feature vector and the microseismic feature vector of each monitoring segment at each time point are concatenated to form the composite feature vector of each monitoring segment at each time point; the composite feature vector is input into the neural network prediction model for prediction to obtain the rockburst level probability of each monitoring segment; the total loss function in the neural network prediction model introduces the mean square error constraint term of the finite element inversion stress distribution and the physical penalty constraint term of the Hoek-Brown strength criterion.

[0051] It should be noted that because single microseismic time series data can only reflect the degree of rock fracturing activity and lacks characterization of the rock's energy storage capacity and mechanical properties, it is necessary to concatenate geological features and microseismic features into a composite feature vector as input. This allows the model to simultaneously acquire information on both the "material basis" and "activity signal" of rockbursts. Furthermore, due to the limited training data and the small sample size of rockburst events, relying solely on data-driven approaches can easily lead to predictions deviating from physical laws. Therefore, to ensure the model output is mechanically reasonable, finite element inversion stress distribution and the Hoek-Brown strength criterion are introduced as two physical constraints to modify the loss function of the neural network model. This allows the model to both fit the data and adhere to mechanical laws during training, ultimately resulting in an improved loss function.

[0052] Specifically, the event rate, energy release rate, and spatial clustering within a time window of each monitoring segment are used to form the microseismic feature vector of each monitoring segment at each time point; the geological feature vector and the microseismic feature vector of each monitoring segment at each time point are concatenated to form the composite feature vector of each monitoring segment at each time point; the composite feature vector is input into a neural network prediction model for prediction to obtain the rockburst level probability of each monitoring segment.

[0053] In this embodiment, the neural network prediction model is an LSTM (Long Short-Term Memory) neural network model. The LSTM neural network model includes an LSTM layer, a fully connected layer, and an output layer; it also includes a separate single-neuron fully connected layer: the hidden state output by the LSTM layer is passed to the single-neuron fully connected layer (1 unit, linear activation), which outputs the estimated surrounding rock stress. Scalar values ​​are used for constraint calculations.

[0054] The total loss function in the LSTM neural network consists of the cross-entropy loss function, the finite element stress calibration loss function, and the Hoek-Brown criterion constraint loss function.

[0055] The total loss function is specifically expressed by the formula:

[0056] In the formula, This represents the cross-entropy loss value. This represents the finite element stress calibration loss value. This represents the constraint loss value of the Hoek-Brown criterion. This represents the first preset weighting coefficient. This represents the second preset weighting coefficient. This represents the total loss value. In this embodiment, the first preset weighting coefficient... Set to 0.1, the second preset weighting coefficient. The value is set to 0.05; where, in this embodiment, the first preset weighting coefficient is... Second preset weight coefficient No specific restrictions are imposed; implementers can decide based on the specific circumstances.

[0057] The cross-entropy loss function is a well-known technique and will not be elaborated on here.

[0058] The finite element stress calibration loss function is specifically expressed as follows:

[0059] In the formula, Indicates the first Estimated surrounding rock stress for each sample Indicates the first The stress values ​​at the corresponding locations of each sample were obtained through finite element inversion calculation; This indicates the number of samples in the training batch.

[0060] The formula represents the mean square error between the stress predicted by the calculation model and the stress calculated by the finite element method. The smaller the error, the closer the stress estimate output by the model is to the mechanical calculation result.

[0061] The Hoek-Brown criterion constraint loss function is specifically expressed as follows:

[0062] In the formula, Indicates the first Estimated surrounding rock stress for each sample The uniaxial compressive strength (scalar, unit MPa) of the rock mass at the location corresponding to the k-th sample is obtained by looking up the lithology table in the geological sketch or by calculating the Hoek-Brown parameter; N represents the number of samples in the training batch. This represents the function that takes the maximum value.

[0063] The formula means that when the estimated stress of the surrounding rock exceeds the uniaxial compressive strength of the rock mass, it indicates that the prediction result does not meet the strength criterion, and a secondary penalty is applied, with the penalty increasing as the stress exceeds the uniaxial compressive strength.

[0064] Among them, a sample refers to the composite feature vector of a monitoring section at a single time point, whether a rockburst will occur in the monitoring section within the next day (24 hours) and the corresponding level label, and the stress value obtained by inversion calculation corresponding to the location of the monitoring section.

[0065] The process for obtaining the estimated surrounding rock stress for each sample in the finite element stress calibration loss function and the Hoek-Brown criterion constraint loss function is as follows: The hidden states output by the LSTM layer in the LSTM neural network model are then passed through a single neuron fully connected layer to obtain the estimated surrounding rock stress for each sample.

[0066] The process of obtaining the stress value at each sample location in the finite element stress calibration loss function through finite element inversion calculation is as follows: Step 1: Obtain the tunnel's geometric parameters from the design drawings, including cross-sectional dimensions and burial depth. Based on these geometric parameters, obtain the geological parameters of each rock stratum using geological sketches and borehole data. These parameters include the boundaries of each rock stratum, fault locations, and dip angles. Input the geological parameters into finite element software to construct a three-dimensional geological model, including stratigraphic interfaces and fault structures. The rock mass material in the three-dimensional geological model uses a Mohr-Coulomb or Hoek-Brown constitutive model.

[0067] Step 2: Obtain the physical and mechanical parameters of each rock layer through indoor rock mechanics tests (cored uniaxial / triaxial compression tests); the physical and mechanical parameters include elastic modulus, Poisson's ratio, density, cohesion, and internal friction angle; assign the physical and mechanical parameters of each rock layer to the corresponding rock layer unit in the three-dimensional geological model; Step 3: Obtain the magnitude and direction of the principal stress at the tunnel burial depth by measuring the in-situ stress using the hydraulic fracturing method, and apply it as the initial stress field boundary condition to the outer boundary of the three-dimensional geological model (the vertical boundary is subjected to horizontal stress, the bottom is fixed, and the upper surface is subjected to load according to the self-weight of the overlying rock layer).

[0068] Step 4: Extract the stress components at the center of each monitoring section using a 3D geological model. Then, decompose each stress component at the center of each monitoring section using the stress tensor eigenvalue, obtaining the maximum principal stress as the stress value calculated by finite element inversion at the center of each monitoring section. The process of decomposing the stress tensor eigenvalue is a well-known technique and will not be described in detail here.

[0069] Step 5: Compare the stress value calculated in Step 4 with the actual stress data measured in the field (such as the stress state corresponding to the concentrated area of ​​microseismic energy release or the measured value of the borehole stress gauge). If the deviation exceeds the allowable range, adjust the input parameters (such as the magnitude of structural stress, lateral pressure coefficient, elastic modulus reduction factor, etc.), return to Step 3, reapply the boundary, and solve again. Iterate repeatedly until the deviation between the calculated value and the measured value meets the requirements (such as...). deviation If the deviation between the calculated and measured values ​​is ≤5%, the model is considered to represent the actual working conditions, and the stress value obtained by finite element inversion under these conditions is output for use in the loss function. In this embodiment, the 5% deviation between the calculated and measured values ​​is not specifically limited and can be determined by the implementer according to the specific circumstances.

[0070] If the deviation exceeds the allowable range in step 5, the input parameters are adjusted; the specific process includes: Step 51: The magnitude of structural stress, lateral compression coefficient, elastic modulus reduction coefficient, etc. are the variables to be inverted, and a search range is set for each variable (e.g., the range of structural stress ±20%). Step 52: Using orthogonal experimental design or Latin hypercube sampling, generate multiple parameter combinations (e.g., 50 combinations) within the search range set for each variable. Substitute each combination into the finite element model for calculation to obtain the corresponding calculated stress value. ; Step 53: Construct a system based on the combination of input parameters and Mapping relationship model between them; Step 54: Within the search range set for each variable, automatically search using an optimization algorithm (such as a genetic algorithm or particle swarm optimization algorithm) to achieve the desired result. Minimize the error with the measured value; the objective function of the optimization algorithm is the error function:

[0071] In the formula, This represents the stress value obtained through finite element inversion calculation for each parameter combination. This represents the measured stress value. It is the absolute value symbol. This indicates the stress value error.

[0072] It's important to note that the search process is iterative. Each time, the parameters are adjusted based on the current error direction to gradually approach the optimal solution, rather than through random attempts. The iteration stops when the error is less than a set threshold (e.g., 5%), and the set of parameters is output as the final "calibrated parameters".

[0073] Step S4: Input the rockburst probability of all monitored sections into the 3D display map for display.

[0074] It should be noted that because the rockburst probability of each section output by the model is an independent and dispersed value, it is difficult to intuitively present the overall risk distribution and relative level of each location. Therefore, in order to facilitate on-site personnel to quickly locate high-risk locations and formulate targeted construction measures, the predicted probability of each monitoring section is projected onto the tunnel 3D space display map in the form of color gradient or bar chart according to its mileage location, so as to realize the visualization of risk distribution.

[0075] Specifically, the mileage markers of each monitoring section and the corresponding rockburst probability levels (minor / moderate / severe) are used as input data; a three-dimensional model of the tunnel is established, and the location of each monitoring section is marked along the tunnel axis; the predicted probability value of each section is projected onto the corresponding location of the three-dimensional model in the form of color mapping (e.g., green = low risk, yellow = medium risk, red = high risk) or bar chart; and a three-dimensional risk distribution map is output for on-site personnel to view.

[0076] Specifically, when mapping the predicted probability value of each segment with color, the process is as follows: When the probability of a rockburst is less than the first-level threshold, the corresponding monitoring section is classified as no rockburst or minor rockburst, and the corresponding color is green; when the probability of a rockburst is less than or equal to the second-level threshold and greater than or equal to the first-level threshold, the corresponding monitoring section is classified as moderate rockburst, and the corresponding color is yellow; when the probability of a rockburst is greater than the second-level threshold, the corresponding monitoring section is classified as severe rockburst, and the corresponding color is red.

[0077] In this embodiment, the first grading threshold is 0.4 and the second grading threshold is 0.7. In this embodiment, the first grading threshold and the second grading threshold are not specifically limited, and the implementer can determine them according to the specific circumstances.

[0078] This concludes the embodiment.

[0079] A second aspect of the present invention is to provide an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement a method for predicting tunnel rockburst levels based on microseismic monitoring and finite element inversion.

[0080] A third aspect of the present invention is to provide a computer-readable storage medium storing a computer program that, when executed by a processor, implements a method for predicting tunnel rockburst levels based on microseismic monitoring and finite element inversion.

[0081] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, optical storage, etc.) containing computer-usable program code.

[0082] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, systems, and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0083] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0084] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0085] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the protection scope of the present invention.

Claims

1. A method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion, characterized in that, include: Geological feature data and microseismic time series data in the tunnel are acquired; a geological feature vector is constructed using the geological feature data; the microseismic time series data includes the event rate, energy release rate, and spatial clustering of each segment within a preset time window; The section is the result of dividing along the tunnel axis; The spatial clustering degree is obtained by calculating the covariance matrix of the source coordinates of each monitored microseismic event and taking the square root of the sum of the three eigenvalues. First, the tunnel is divided into initial sections using geological feature data to obtain an initial section set. Then, based on the differences in microseismic time series data between adjacent initial sections in the initial section set, the division results of the initial sections are adjusted to obtain several adjusted monitoring sections. The event rate, energy release rate, and spatial clustering within a time window of each monitoring segment are used to form the microseismic feature vector of each monitoring segment at each time point. The geological feature vector and the microseismic feature vector of each monitoring segment at each time point are concatenated to form the composite feature vector of each monitoring segment at each time point. The composite feature vector is input into a neural network prediction model for prediction to obtain the rockburst level probability of each monitoring segment. The total loss function in the neural network prediction model incorporates the mean square error constraint term of the finite element inversion stress distribution and the physical penalty constraint term of the Hoek-Brown strength criterion. The probability of rockburst levels in all monitored sections is input into a 3D visualization map for display.

2. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 1, characterized in that, The specific process for acquiring the microseismic time series data includes: The event rate is calculated by counting the total number of micro-seismic events in each segment within a preset time window; the energy release rate is calculated by counting the sum of the micro-seismic energy of all micro-seismic events in each segment within a preset time window. The process of obtaining the spatial clustering degree is as follows: Obtain the source coordinates of all microseismic events within each segment; form a symmetric covariance matrix based on the covariances of the source coordinates of all microseismic events in the three coordinate directions; solve the characteristic equation of the symmetric covariance matrix to obtain the three eigenvalues. , , The sum of the three factors, then the square root, gives the spatial density. ; The covariance symmetric matrix is ​​specifically represented as follows: In the formula, This represents the covariance of the x-coordinate itself. This represents the covariance between the x-coordinate and y-coordinate directions. This represents the covariance between the x-coordinate and the z-coordinate. This represents the covariance between the y-coordinate and the z-coordinate. This represents the covariance along the y-coordinate direction itself. This represents the covariance of the z-coordinate itself. in, , , .

3. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 1, characterized in that, The process of obtaining the initial segment set includes: Along the tunnel mileage, arrange the geological feature data in the order of each cycle of advance, and check whether there are abrupt changes in lithology, fault markers, and structural surface density. If any of the lithology, fault markers, and structural surface density changes, mark it and record the marked position as the initial boundary. An initial segment is formed between adjacent initial boundaries. After all initial boundaries are marked, output the set of initial segments.

4. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 1, characterized in that, The specific acquisition process for the aforementioned monitoring sections includes: The event rate, energy release rate, and spatial clustering degree within each initial segment are statistically analyzed. Based on these three parameters, the microseismic response difference between adjacent initial segments is calculated. When the microseismic response difference between adjacent initial segments is greater than or equal to a preset difference threshold, the microseismic index difference between the adjacent initial segments is determined to be significant. When the microseismic response difference between adjacent initial segments is less than the preset difference threshold, the microseismic index difference between the adjacent initial segments is determined to be insignificant. If the microseismic indices of adjacent initial segments differ significantly, the corresponding boundary is retained; if the microseismic indices of adjacent initial segments do not differ significantly, the corresponding boundary is deleted and the two initial segments are merged. The difference in microseismic response between adjacent initial sections is specifically expressed by the following formula: In the formula, This represents the event rate of the first initial segment. This represents the event rate of the second initial segment. This represents the energy release rate of the first initial segment. This indicates the energy release rate of the second initial segment. This indicates the spatial clustering degree of the first initial segment. This indicates the spatial clustering degree of the second initial segment; This indicates the degree of difference in microseismic response between adjacent initial sections. Represents a linear normalization function. It is the absolute value symbol.

5. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 1, characterized in that, The neural network prediction model is an LSTM neural network model; The LSTM neural network model includes an LSTM layer, a fully connected layer, and an output layer; it also includes an independent single-neuron fully connected layer.

6. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 5, characterized in that, The total loss function consists of the cross-entropy loss function, the finite element stress calibration loss function, and the Hoek-Brown criterion constraint loss function. The total loss function is specifically expressed by the formula: In the formula, This represents the cross-entropy loss value. This represents the finite element stress calibration loss value. This represents the constraint loss value of the Hoek-Brown criterion. This represents the first preset weighting coefficient. This represents the second preset weighting coefficient. This represents the total loss value; The finite element stress calibration loss function is specifically expressed as follows: In the formula, Indicates the first Estimated surrounding rock stress for each sample Indicates the first The stress values ​​at the corresponding locations of each sample were obtained through finite element inversion calculation; Indicates the number of samples in the training batch; The Hoek-Brown criterion constraint loss function is specifically expressed as follows: In the formula, Indicates the first Estimated surrounding rock stress for each sample This represents the uniaxial compressive strength of the rock mass at the location corresponding to the kth sample, which is obtained by looking up the lithology table in the geological sketch or by calculating the Hoek-Brown parameter. Indicates the number of samples in the training batch; This represents the function that takes the maximum value. Among them, a sample refers to the composite feature vector of a monitoring section at a single time point, whether a rockburst will occur in the monitoring section in the next day and the corresponding level label, and the stress value obtained by inversion calculation corresponding to the location of the monitoring section. The process of obtaining the estimated surrounding rock stress for each sample in the finite element stress calibration loss function and the Hoek-Brown criterion constraint loss function is as follows: the hidden state output by the LSTM layer in the LSTM neural network model is then passed through a single neuron fully connected layer to obtain the estimated surrounding rock stress for each sample.

7. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 6, characterized in that, The process of obtaining the stress value at each sample location in the finite element stress calibration loss function through finite element inversion calculation is as follows: Step 1: Obtain the tunnel's geometric parameters from the design drawings. These parameters include cross-sectional dimensions and burial depth. Based on these geometric parameters, obtain the geological parameters of each rock layer using geological sketches and borehole data. These geological parameters include the boundaries of each rock layer, fault locations, and dip angles. Input these geological parameters into finite element software to construct a three-dimensional geological model. Step 2: Obtain the physical and mechanical parameters of each rock layer through indoor rock mechanics tests; the physical and mechanical parameters include elastic modulus, Poisson's ratio, density, cohesion, and internal friction angle; assign the physical and mechanical parameters of each rock layer to the corresponding rock layer unit in the three-dimensional geological model; Step 3: Obtain the magnitude and direction of the principal stress at the tunnel burial depth by measuring the in-situ stress using the hydraulic fracturing method, and apply it as the initial stress field boundary condition to the outer boundary of the three-dimensional geological model. Step 4: Extract the stress components at the center of each monitoring section using a three-dimensional geological model. The stress components at the center of each monitoring section are then decomposed using the stress tensor eigenvalues ​​to obtain the maximum principal stress, which is used as the stress value calculated by finite element inversion at the center of each monitoring section. Step 5: Compare the stress value calculated in Step 4 with the actual stress data measured on site. If the deviation exceeds the allowable range, adjust the input parameters, return to Step 3 to reapply the boundary and solve again. Iterate repeatedly until the deviation between the calculated value and the measured value meets the requirements. At this point, the model is considered to represent the actual working condition. Output the stress value calculated by finite element inversion under this working condition for use in the loss function.

8. The method for predicting tunnel rockburst severity based on microseismic monitoring and finite element inversion according to claim 7, characterized in that, If the deviation in step 5 exceeds the allowable range, the input parameters are adjusted; the specific process includes: Step 51: The magnitude of structural stress, lateral compression coefficient, and elastic modulus reduction coefficient are the variables to be inverted, and a search range is set for each variable; Step 52: Using orthogonal experimental design or Latin hypercube sampling, generate multiple sets of parameter combinations within the search range set for each variable. Substitute each set of combinations into the finite element model for calculation to obtain the corresponding calculated stress value. Step 53: Construct a mapping model between the input parameter combination and the calculated stress values; Step 54: Within the search range set for each variable, an optimization algorithm is used to automatically search to minimize the error between the calculated stress value and the measured value; the objective function of the optimization algorithm is the error function. In the formula, This represents the stress value obtained through finite element inversion calculation for each parameter combination. This represents the measured stress value. It is the absolute value symbol. Indicates stress value error; The optimization algorithm is either a genetic algorithm or a particle swarm optimization algorithm.

9. An electronic device, characterized in that, The method includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the tunnel rockburst level prediction method based on microseismic monitoring and finite element inversion as described in any one of claims 1-8.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the tunnel rockburst level prediction method based on microseismic monitoring and finite element inversion as described in any one of claims 1-8.

Images