A method for monitoring fissure field evolution of slicing mining of super-thick coal seam under rock burst

Abstract

The application discloses a kind of special thick coal seam layered mining fissure field evolution monitoring methods under rock burst, method includes: using stress correction coefficient to dynamic adjustment wave velocity-fracture opening, obtain dynamic fracture opening;Combining dynamic fracture opening, real-time stress and equivalent permeability, reconstruct to obtain time series crack state feature sequence;Mining disaster risk prediction model is constructed and trained, five-dimensional crack-stress tensor and time series crack state feature sequence are input into mining disaster risk prediction model, after model coupling deduction, output four-dimensional risk probability distribution tensor and physical credibility index of future time;Track the volume change sequence of potential disaster body, calculate to obtain bite intensity index and disaster crack predicted breakthrough time, finally carry out disaster risk grade division, finally output disaster risk grade, leading disaster type and disaster crack predicted breakthrough time;Intelligent monitoring is realized and the accuracy of prediction is improved.

Description

Technical Field

[0001] This invention relates to the field of coal mine safety technology, and relates to, but is not limited to, a method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst. Background Technology

[0002] In the process of coal resource mining, the efficient development of extra-thick coal seams (usually referring to those with a thickness exceeding 8 meters) is of great significance to ensuring the country's energy supply, and layered mining technology has become the mainstream process for such coal seams. However, with the increase in mining depth and intensity, layered mining of extra-thick coal seams faces the coupled risks of "rockburst" and "spontaneous combustion of coal seams"; the complex structure of the recycled roof formed by upper layered mining can easily induce rockbursts when the lower layered mining is carried out; at the same time, the coal pillars, roof coal fracture zones, and vertical fissures generated by roof fractures left by layered mining form air leakage and oxygen supply channels, which can easily cause spontaneous combustion of coal seams in coal seams with strong natural tendency to spread, as oxygen seeps into the deep goaf.

[0003] Existing fiber-optic microseismic monitoring systems can only invert the type, intensity, and location of microseismic events, but cannot distinguish between rockbursts and coal seam spontaneous combustion hazards, nor can they respond in real time to the dynamic evolution of regenerated roof fractures caused by layered mining. Therefore, it is necessary to achieve accurate inversion of the dynamic evolution of the fracture field and disaster risk classification through multi-source data fusion and physical model correction. Summary of the Invention

[0004] In view of this, embodiments of the present invention provide a method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst, aiming to overcome the shortcomings of the prior art. The prior art cannot respond in real time to the dynamic evolution of the regenerated roof stress field and the fracture failure activation process caused by the mining of the lower layer in layered mining, and cannot distinguish the disaster type, resulting in inaccurate prevention and control strategies.

[0005] The specific technical solutions of this invention are as follows: This invention provides a method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst, comprising: Acquire microseismic monitoring point data and stress data, process the microseismic monitoring point data and stress data to obtain the five-dimensional crack-stress tensor; Based on real-time microseismic wave velocity data and real-time fiber optic strain data, a stress correction coefficient is introduced; the wave velocity-fracture aperture is dynamically adjusted using the stress correction coefficient to obtain the dynamic fracture aperture. Acquire fiber stress data, calculate the fiber stress data to obtain real-time stress; combine the dynamic fracture aperture, real-time stress and equivalent permeability to reconstruct the time-series fracture state characteristic sequence. A mining disaster risk prediction model is constructed and trained. The five-dimensional fracture-stress tensor and the temporal fracture state feature sequence are input into the mining disaster risk prediction model. After coupling and extrapolation by the model, the four-dimensional risk probability distribution tensor and physical reliability index at future time are output. The mining disaster risk prediction model sets up a static background encoding branch and a dynamic damage encoding branch. The static background encoding branch extracts the spatial features of the five-dimensional fracture-stress tensor, and the dynamic damage encoding branch extracts the temporal evolution features of the temporal fracture state feature sequence. After the two features are encoded by the physical information embedding layer to encode the fracture geometry-seepage constitutive relationship, the data-driven features and physical constraint features are adaptively fused through a physical consistency gating network, and then the prediction result is output after correction by a physical correction head. The four-dimensional risk probability distribution tensor and the physical reliability index are subjected to multi-channel separation processing to extract the potential disaster bodies in each channel. The volume change sequence of the potential disaster bodies is tracked to calculate the interlocking force index and the expected penetration time of the disaster crack. The disaster risk level is classified according to the interlocking force index, the expected penetration time of the disaster crack and the risk probability of each channel. Finally, the disaster risk level, the dominant disaster type and the expected penetration time of the disaster crack are output.

[0006] In some embodiments, the multi-channel separation processing of the four-dimensional risk probability distribution tensor and the physical credibility index to extract potential hazard bodies for each channel includes: By combining the four-dimensional risk probability distribution tensor with the physical credibility index, multi-channel separation processing is performed. A connected component analysis algorithm is used to extract continuous regions in each channel whose probability exceeds a preset threshold, and these continuous regions are identified as potential disaster bodies.

[0007] In some embodiments, the step of classifying disaster risk levels based on the interlocking force index, the estimated breakthrough time of the disaster fissure, and the risk probability of each channel, and finally outputting the disaster risk level, the dominant disaster type, and the estimated breakthrough time of the disaster fissure, includes: When the impact-spontaneous combustion coupling probability is greater than the first probability threshold and the biting force index is greater than the first force threshold, it is classified as a level one disaster risk; When the impact-spontaneous combustion coupling probability is less than the second probability threshold, the rockburst probability is greater than the first probability threshold, or the coal seam spontaneous combustion probability is greater than the first probability threshold, and the interlocking force index is less than the first force threshold, it is classified as a level two disaster risk. When the safety probability, rockburst probability, coal seam spontaneous combustion probability, and rockburst-spontaneous combustion coupling probability are all greater than the second probability threshold and less than the first probability threshold, or when the biting force index is less than the second force threshold, the risk is classified as Level III disaster risk.

[0008] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following: In this embodiment of the invention, firstly, a background tensor containing parameters such as fracture density and stress is constructed using microseismic and fiber optic data. The fracture aperture and permeability are dynamically inverted through a stress correction model to form a temporal feature sequence. A dual-stream encoding-fusion-decoding deep learning architecture is designed to extract the spatiotemporal features of static geological background and dynamic damage evolution, respectively, and output a four-dimensional disaster probability distribution. The "biting force" index is introduced to quantify the spatial correlation between upper and lower layered fractures. Combined with probability thresholds, a three-level risk classification rule is established to realize intelligent monitoring of the entire process from disaster identification to risk management.

[0009] By integrating multi-source monitoring data with a physical inversion model, a dynamic inversion mechanism for fracture characteristics based on stress correction and a dual-stream encoding-fusion-decoding architecture embedding physical information were constructed. By introducing physical information neural network constraints, prior knowledge of rock mechanics such as mass conservation, momentum conservation, and the cubic law were explicitly embedded into the deep learning model, solving the problems of poor physical consistency and weak extrapolation ability in purely data-driven models. Through the design of a physical consistency gating network and an adaptive weight scheduling strategy, a dynamic balance between data fitting and physical constraints was achieved, avoiding training convergence difficulties caused by premature introduction of physical priors. Through a physical correction head and post-processing mechanism, the prediction results were ensured to strictly obey the laws of continuum mechanics. Attached Figure Description

[0010] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments 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, wherein: Figure 1 A schematic flowchart illustrating the method for monitoring the fracture field evolution in layered mining of extra-thick coal seams under rockburst, as provided in an embodiment of the present invention. Figure 2 This is a model architecture diagram of the method for monitoring the fracture field evolution in layered mining of extra-thick coal seams under rockburst, provided in an embodiment of the present invention. Detailed Implementation

[0011] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. The following embodiments are used to illustrate the present invention, but are not intended to limit the scope of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0012] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.

[0013] It should be noted that the terms "first, second, and third" used in the embodiments of the present invention are only used to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, and third" can be interchanged in a specific order or sequence where permitted, so that the embodiments of the present invention described herein can be implemented in an order other than that illustrated or described herein.

[0014] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which these embodiments of the invention pertain. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.

[0015] Figure 1 This is a flowchart illustrating a method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst, as provided in an embodiment of the present invention. Figure 1 As shown, the method includes at least the following steps: Step S110: Process the microseismic monitoring point data and stress data to obtain the five-dimensional fracture-stress tensor.

[0016] Historical monitoring data was collected during the upper-layer mining of Coal Mine A, covering three working faces and spanning approximately 90 days. Through underground total station coordinate reference network and fiber optic grating spatial positioning calibration, the three-dimensional coordinates of microseismic events and the coordinates of fiber optic strain gauge points were unified into the same mine coordinate system. Nearest neighbor interpolation and inverse distance weighted fusion were used to achieve spatial registration of multi-source data to construct a static five-dimensional background fracture-stress tensor. To address the issues of missing or unevenly sampled historical data in the upper layers, a spatiotemporal interpolation algorithm based on the physical constraints of fracture evolution was used to reasonably fill in the missing areas. Microseismic monitoring and distributed fiber optic sensors were used to acquire microseismic monitoring point data and stress data. For microseismic monitoring point data and stress data, an inverse distance weighted interpolation algorithm is used to interpolate and map discrete microseismic monitoring points to a three-dimensional grid space, constructing a continuous physical field throughout the space. Five key physical parameters—fracture density, fracture aperture, equivalent permeability, historical stress concentration factor, and fracture attitude tensor—are calculated and stacked in three-dimensional space to construct a five-dimensional tensor. Simultaneously, a fracture aperture-permeability physical relationship matrix is ​​constructed based on the cubic law, providing a priori constraints for subsequent physical information embedding layers.

[0017] During layered mining in the mine, historical data from the entire mining cycle is reviewed. This historical data refers to data from a period prior to the current moment. The three-dimensional spatial coordinates (X, Y, Z), radiation energy E, and seismic moment of each microseismic event are extracted. Simultaneously, high-frequency fiber optic stress monitoring data is acquired through a distributed fiber optic sensor network deployed at key layers of the roof. The monitoring area is spatially divided into... A three-dimensional voxel mesh is constructed, and an inverse distance-weighted interpolation algorithm is used, with the reciprocal of the distance as the weight, to map and assign discrete point-like microseismic data and linear fiber optic data to each mesh voxel. .

[0018] First, calculate the crack density using the following formula: in, Crack density, For the frequency of microseismic events, Voxel volume; The formula for calculating crack opening is as follows: in, For crack opening, For seismic moment, is the rock mass shear modulus (a preset value in the mining area geomechanical parameter database), and A is the fracture surface area; Using an equivalent area strategy, the equivalent radius r of the dominant rupture surface is obtained through micro-moment tensor inversion, in order to... Calculate the equivalent fracture surface area and introduce the fracture surface shape factor. (Usually taken as 1.2~2.5, determined by laboratory acoustic emission calibration) is corrected, i.e., A_eff=ξ·A, where ξ is the fracture surface shape factor; Substitute again ;in, The corrected fracture surface; It should be noted that during the dynamic inversion stage, fiber optic strain monitoring data is used as a hard constraint to perform closed-loop correction on the initial value of δ.

[0019] based on and The equivalent value of permeability is derived using the following formula: in, This is the equivalent value of penetration rate. Porosity is C, and shape factor is C. For the degree of curvature, The crack opening; In regions with relatively uniform fracture distribution and good connectivity, the method can be further simplified to: ; Historical stress concentration factor ,in, The peak stress within the grid. This represents the regional background stress. It should be noted that the fracture orientation tensor The strike, dip, and dip angle of the fracture surface are extracted by calculating the moment tensor inversion results of microseismic events and using eigenvalue decomposition.

[0020] The formula for calculating the rock mass splitting index is as follows: in, This represents the rock mass splitting index. This represents the cumulative increment of energy released during micro-seismic events within the current window. This is equivalent to the average energy release for the same period in history. The local b value of the current grid. The background b value for the region; It should be noted that the b-value represents the slope of the magnitude-frequency distribution line, which is a core parameter in the field of rock mechanics.

[0021] The five physical parameters calculated above are stacked in three-dimensional space to construct the background crack-stress five-dimensional tensor of the regenerated roof. .

[0022] Step S120: Dynamically adjust the wave velocity-fracture aperture using the stress correction coefficient to obtain the dynamic fracture aperture; combine the dynamic fracture aperture, real-time stress, and equivalent permeability to reconstruct the time-series fracture state characteristic sequence.

[0023] By integrating real-time microseismic wave velocity data and real-time fiber optic strain data, a stress correction coefficient is introduced to dynamically weight and adjust the classical empirical formula for wave velocity-fracture aperture, thereby inverting and calculating the real-time dynamic fracture aperture. The equivalent permeability is calculated using the cubic law, and the real-time stress state is solved based on fiber optic data. Finally, the dynamic fracture aperture, equivalent permeability, and real-time stress state are combined into a time-series fracture state characteristic sequence that evolves over time, achieving dynamic physical field reconstruction of rock mass damage evolution. Simultaneously, mass conservation residuals and momentum conservation residuals are calculated based on the inversion results, and a physical consistency monitoring signal set is constructed for physical information constraints during the training phase.

[0024] It should be noted that the microseismic monitoring system captures the elastic wave signal generated by rock mass fracture in real time, extracts the first arrival wave travel time and amplitude attenuation characteristics, and the distributed optical fiber sensor network (such as FBG grating) collects the grating wavelength offset at high frequency.

[0025] The formula for the stress correction factor is as follows: in, This is the stress correction factor; This is the empirical benchmark coefficient; The real-time principal stress is inverted from the grating wavelength offset; The initial background stress; This is the lithological stress sensitivity coefficient; It should be noted that the lithological stress sensitivity coefficient was calibrated using laboratory triaxial compression-seepage experiments, with different lithological ranges (sandstone 0.05~0.15, siltstone 0.15~0.30, mudstone and coal 0.30~0.60), to establish the lithological-stress sensitivity coefficient. Corresponding database. The fracture width is measured downhole using borehole cameras and then verified using least-squares inversion. The recommended calibration period is every 50-100m of drilling or when traversing zones of significant lithological change.

[0026] The formula for calculating dynamic fracture aperture is as follows: in, For dynamic fracture aperture, The initial fracture aperture, The theoretical wave velocity of the intact rock mass, The wave velocity is monitored in real time; the exponent "2" comes from the theory of equivalent wave velocity in fractured media. This is the stress correction factor.

[0027] Based on the cubic law, the dynamic fracture aperture is converted into equivalent permeability. : in, For equivalent penetration rate, For dynamic crack aperture; The formula for calculating real-time stress is as follows: in, For real-time stress, The initial stress; For stress increment, ,in, The elastic modulus of the rock mass. For real-time response.

[0028] Based on the physical field data of the three key dimensions obtained from the above inversion, the fracture type is determined. like and If the stress level is below the low stress threshold, it is classified as an open fracture. The low stress threshold is... ; like and When the stress exceeds the high stress threshold, it is determined to be a closed crack; the high stress threshold is... ; like or and When the combination of [condition] does not satisfy the above two hard classification thresholds, it is determined to be a transitional fracture; among which, The uniaxial compressive strength of the rock mass; It should be noted that, The crack density threshold is determined empirically; the uniaxial compressive strength of the rock mass is a rock mechanics parameter obtained from rock mechanics experiments.

[0029] If a fracture is identified as a transitional fracture, it is no longer forcibly classified into a single category. Instead, the continuous probability distribution of its classification into each fracture category is output through a Softmax layer in the mining disaster risk prediction model. ( The probability of a complete rock mass. For closed fracture probability, This represents the probability of an open-type crack. The probability of the disruption / instability region satisfies ).

[0030] It should be noted that, Characterizing the mechanical state of the crack is a soft classification based on physical parameter thresholds. Its purpose is to avoid information loss caused by hard classification and to provide data-driven fine-grained discrimination for deep learning modules; p represents disaster risk and is used for subsequent risk classification.

[0031] Real-time acquisition of microseismic and stress data at T time steps; reconstruction of the temporal fracture state characteristic sequence from the inversion results of each time step t. .

[0032] Step S130: Construct and train the mining disaster risk prediction model. Input the five-dimensional fracture-stress tensor and the temporal fracture state feature sequence into the mining disaster risk prediction model. After model coupling and deduction, output the four-dimensional risk probability distribution tensor and physical credibility index for future time moments.

[0033] It should be noted that the constructed background tensor and temporal feature sequence are input into the trained mining disaster risk prediction model. During inference on the edge device, the lightweight student network directly replaces the teacher network to perform forward propagation. The background tensor calls the pre-computation cache, and the temporal sequence is quickly encoded by 1D depthwise separable convolution. After being fused by the physical information embedding layer and gating, the 3D transposed convolution implemented by group convolution maps back to the spatial resolution, and finally outputs a four-dimensional probability distribution tensor. Simultaneously output physical consistency indicators. Quantify the physical reliability of the current prediction results, whereby, Indicates physical loss. This represents the physical loss threshold.

[0034] A spatiotemporal deep learning framework based on physical information embedding, employing a dual-stream encoding-fusion-decoding approach, employs a static background encoding branch to extract spatial features of a five-dimensional background tensor using a 3D-CNN, and a dynamic damage encoding branch to extract the temporal evolution of temporal feature sequences using a hybrid architecture of 1D-CNN and Transformer. The static background encoding branch extracts spatial features of a five-dimensional fracture-stress tensor, while the dynamic damage encoding branch extracts temporal evolution features of a temporal fracture state feature sequence. These two feature paths are encoded through a physical information embedding layer to represent the fracture geometry-seepage constitutive relationship. Then, a physical consistency gating network adaptively fuses data-driven features and physical constraint features. After physical correction, the data is mapped back to the original spatial resolution via a 3D transposed convolutional network, outputting a four-dimensional risk probability distribution tensor and a physical reliability index representing four future states: safety, rockburst, coal seam spontaneous combustion, and coupled disasters.

[0035] Based on a two-stream coding-fusion-decoding architecture, and designed for real-time inference with a teacher-student dual-network architecture, the teacher network is used for offline high-precision training, while the student network, trained through knowledge distillation, is used for lightweight real-time inference at the edge. Both networks share the same physical information embedding and fusion-decoding logic. The two-stream coding includes a static background coding branch and a dynamic impairment coding branch.

[0036] It should be noted that, in order to meet the requirements of real-time inference at the edge, a lightweight student network is trained using knowledge distillation, which reduces the number of parameters by about 78% while maintaining an accuracy of over 95%.

[0037] Specifically, the static background encoding branch employs a 3D convolutional neural network, using a five-dimensional background crack-stress tensor. As input, output a high-dimensional static feature map. The dynamic damage encoding branch employs a hybrid architecture of 1D convolutional neural network (1D-CNN) and Transformer to achieve... As input, output feature vector .

[0038] To meet the computational constraints of real-time inference at the edge, the student network architecture replaces the standard 3D-CNN with a static background encoding branch with a 3D depthwise separable convolution (3D-MobileNetV3 backbone), and simplifies the Transformer to a single-layer 4-head attention mechanism; the 3D transposed convolution at the decoding end is implemented using grouped convolution.

[0039] It should be noted that after knowledge distillation training, the student network reduced the number of parameters by about 78% and the single-step inference computation by about 85% while maintaining more than 95% of the prediction accuracy of the teacher network.

[0040] Before feature fusion, it is explicitly agreed that a physical information embedding layer will be introduced during the feature fusion stage. This is achieved through a learnable physical encoder. The static background fracture geometry parameters and the dynamically inverted fracture state parameters are jointly mapped to the physical embedding space.

[0041] in, This is a time-series feature map; Crack density, For physical encoders, For the background crack opening displacement, It is taken from the five-dimensional background tensor; Taken from the constructed temporal fracture state feature sequence middle; This refers to the rock mass shear modulus. Specifically, the physical encoder uses a 3-layer perceptron (MLP) to process the 7-dimensional physical parameters. Mapped to 3D physical embedding vector (usually) The embedding vector is mapped to the channel dimension C through a fully connected layer and then broadcast as a static feature. Figure 1 Achieving spatial resolution and generating three-dimensional physical feature maps. The 1D temporal feature vector is mapped using a time-space mapping module. After mapping to channel dimension C, the broadcast is the same as the static feature. Figure 1 To achieve the desired spatial resolution, generate temporal feature maps. .

[0042] The fusion process employs a physical consistency gating mechanism, expressed as follows: in, To integrate the features of the gating mechanism; The fused feature is obtained by concatenating static feature maps and temporal feature maps in the channel dimension and then compressing them through a 1×1×1 three-dimensional convolution. Physical feature map; For the output of the gated network, where, For the Sigmoid function, This is a vector concatenation operation. and These are learnable parameters; After fusion, a cross-attention mechanism is used, where Query(Q) comes from the dynamic temporal feature branch. Key(K) and Value(V) come from the static background feature branch. ; Finally, a 3D transposed convolutional network is used to map the fused features back to the original spatial resolution, and a four-dimensional probability distribution is output through a physical correction head. probability distribution ,in, For safety probability, This represents the probability of a rockburst. This represents the probability of spontaneous combustion of the coal seam. The impact-spontaneous combustion coupling probability is shown in the model architecture diagram below. Figure 2 As shown.

[0043] Step S140: Perform multi-channel separation processing on the four-dimensional risk probability distribution tensor and physical reliability index, extract the potential disaster bodies of each channel, and finally classify the disaster risk level, output the disaster risk level, the dominant disaster type and the expected breakthrough time of the disaster fissure.

[0044] Based on four-dimensional probability distribution and physical reliability indicators, a connectivity domain analysis algorithm is used to identify potential hazard bodies and invert their evolution time. The vertical overlap of upper and lower layered fractures is calculated, and an interlocking strength index is defined and calculated to quantify spatial topological correlation. Finally, combining hazard probability and interlocking strength index, a three-level risk classification judgment rule is established, and a comprehensive assessment report including risk level, main hazard type, expected breakthrough time, and physical reliability is output.

[0045] Channel separation is performed on the four-dimensional probability tensor, namely the probability of safety, the probability of rockburst, the probability of spontaneous combustion of coal seam, and the probability of rockburst-spontaneous combustion coupling. Using a connected component analysis algorithm, continuous regions in each channel whose probability exceeds a preset threshold (0.6) are extracted and identified as potential hazard bodies. For each hazard body, its geometric centroid coordinates are calculated. The specific formula is as follows: , , in, The coordinates of the geometric centroid of the disaster body are... For the first disaster in the body The coordinates of the individual elements For disaster body The probability value of a voxel. The total number of voxels contained in the disaster body; ,in, The equivalent volume of the disaster body. The total number of voxels contained in the disaster body. The physical volume of a single voxel.

[0046] By tracking the volume change sequence of the same disaster body in multi-time prediction results, and nonlinearly fitting the crack propagation trend based on the volume evolution sequence, the predicted crack penetration time of the disaster body to reach the critical penetration volume is solved.

[0047] By quantifying the vertical overlap and interaction potential of the upper and lower layer fractures, the interlocking force index is obtained, and the specific formula is as follows: in, This refers to the bite force index. The overlapping volume of the upper and lower layered disaster bodies. The union volume of the upper and lower layered disaster bodies; The vertical distance between the centers of mass of the upper and lower stratified disaster bodies; is the standard interlayer spacing, and k is the attenuation coefficient.

[0048] It should be noted that a three-level risk classification judgment rule is established based on probability thresholds: Level 1 is extremely high risk - coupled and interconnected, Level 2 is high risk - single dominant, and Level 3 is of concern - potential risks exist.

[0049] Specifically, the three-level risk classification determination rules are as follows: When the impact-spontaneous combustion coupling probability is greater than the first probability threshold and the biting force index is greater than the first force threshold, it is judged as a level one disaster risk.

[0050] It should be noted that the first probability threshold is 0.7, and the first intensity threshold is 0.6; Level 1 risk is extremely high / coupling and connected, indicating that the upper and lower layered fracture networks have formed a vertically connected channel, and the risks of rockburst and spontaneous combustion exist simultaneously and are highly coupled. Immediate production shutdown and evacuation of personnel are required, along with the activation of emergency response control measures.

[0051] When the impact-spontaneous combustion coupling probability is less than the second probability threshold, and the rockburst probability is greater than the first probability threshold or the coal seam spontaneous combustion probability is greater than the first probability threshold, and the interlocking force index is less than the first force threshold, it is judged as a level two disaster risk.

[0052] It should be noted that the second probability threshold is 0.3; the secondary risk level is high-risk / single-dominant, indicating that the area is either a pure rockburst or a pure spontaneous combustion zone, with relatively independent disaster evolution. Specific remediation measures are implemented for each specific disaster type. At this time, it is defined as a level 2 high risk, indicating that although the fractures in the upper and lower layers are not completely connected, it is necessary to strengthen the monitoring of fracture connection on the basis of single disaster management.

[0053] When all channels—namely, the probability of safety, the probability of rockburst, the probability of spontaneous combustion of coal seam, and the probability of rockburst-spontaneous combustion coupling—are greater than the second probability threshold and less than the first probability threshold, and are between 0.3 and 0.7, or when the biting force index is less than the second force threshold, the disaster risk is determined to be a level three disaster risk.

[0054] It should be noted that when the second strength threshold is 0.5, the third level of risk is "attention" or "potential danger," indicating that the crack is in the early or stable stage of development and has not yet formed a penetrating failure; the monitoring cycle should be shortened to half of the standard third level.

[0055] The disaster risk is classified, and the final output includes the risk level, the dominant disaster type, and the expected time for the disaster fissure to be connected.

[0056] The following describes the method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst using a specific embodiment. However, it is worth noting that this specific embodiment is only for better illustrating the present invention and does not constitute an improper limitation of the present invention.

[0057] Based on historical data from the upper strata of monitored coal mines, a five-dimensional background fracture-stress tensor was constructed. Real-time acquisition of microseismic and stress data at 10 time steps in the lower layer was performed, and a time-series fracture state characteristic sequence was constructed. .

[0058] For example, using historical data from the entire mining cycle (approximately one month) of the same coal mine, a temporal fracture state feature sequence was constructed with 10 time steps (each time step window being one minute) as one sample. A joint annotation strategy of sparse hard labels and semi-supervised pseudo labels was adopted. Approximately 320 hard-labeled samples were selected from key areas and obtained through non-destructive testing with ground-penetrating radar, sparse borehole photography (one verification borehole every 50m), and high-precision stress monitoring as seed samples. Using the hard labels as seeds, and combining the temporal continuity constraints and physical consistency criteria of fracture evolution, an iterative pseudo-labeling method was used to propagate label confidence to unlabeled samples, generating pseudo labels. After manual review and physical rule filtering, approximately 4300 samples were formed (Class 0 55.8%, Class 1 20.0%, Class 2 16.7%, Class 3 7.4%).

[0059] Specifically, a dataset is constructed based on the temporal fracture state feature sequence. An iterative pseudo-label method is used to generate a core training sample set with four-class labels. The core training sample set is then divided into a training set, a validation set, and a test set through stratified sampling. Simultaneously, a physical constraint sample set is constructed based on the numerical simulation results of fracture network seepage for physical supervision training of the model. During model training, a loss function is introduced, which includes weighted cross-entropy classification loss, spatial smoothing loss, physical conservation loss, and physical consistency regularization term. The physical conservation loss includes mass conservation error term and momentum conservation error term, and the physical consistency regularization term constrains the model to maintain a cubic law relationship between the predicted equivalent permeability and fracture aperture.

[0060] Introducing weighted cross-entropy into the loss function ( , The model employs a minority class bias, assigning different weights to different crack categories using weighted cross-entropy classification loss; spatial smoothing loss constrains the rate of change of risk probability between adjacent voxels; and random oversampling and SMOTE hybrid enhancement are applied to the three classes of samples during training to ensure the model's sensitivity and recall to disaster categories. Samples are divided into training, validation, and test sets in a 7:2:1 ratio, using a stratified sampling strategy.

[0061] The weighting coefficients of the physical conservation loss adopt an adaptive scheduling strategy. In the early stage of training, data fitting is the main focus, and physical constraints are gradually strengthened in the later stage. The weighting coefficients of the adaptive scheduling strategy increase with the training rounds. The trained prediction model is the mining disaster risk prediction model.

[0062] A physical constraint sample set was constructed: the theoretical seepage flow rate was calculated based on the numerical simulation of seepage flow in the fracture network and indirectly verified by the gas concentration monitoring data of the goaf, thus constructing a mass conservation monitoring signal; the momentum conservation residual was calculated based on the stress field, thus constructing a mechanical conservation monitoring signal.

[0063] Specifically, the loss function consists of four parts: ; c is the fracture category index (c=0 for intact rock mass, c=1 for closed fracture, c=2 for open fracture, c=3 for failure / instability zone). Category weights ( , , 4.0 ), To predict the probability, satisfying .

[0064] ; For spatial smoothing loss, for The overall probability of a voxel experiencing a disaster (unsafe state); , , and These represent the difference operators in the three spatial dimensions X, Y, and Z, used to calculate the rate of change of risk probability between adjacent voxels. Smoothing factor (usually taken as...) ); For time-varying physical loss weighting coefficients, an adaptive scheduling strategy is adopted: in, For training rounds, The activation threshold round is determined by physical constraints. This is the annealing coefficient. The initial weight (usually 0.01). This is the maximum weight (usually 0.5).

[0065] in, This is the physical loss function; the first term is the mass conservation error, and the second term is the momentum conservation error. This represents the total number of sampling points; The square of the L2 norm; For the density of the rock mass, For seepage velocity field, It is the first The mass flux vector formed by the product of density and velocity at each sampling point; For mass flux in the first The divergence at each point; For the first Source and sink items at each point, For the first Stress tensor at a point It is the first The divergence at each point; For the first Volume force at each point For the first Acceleration field at each point; It is the product of density and acceleration, representing the inertial force term.

[0066] It should be noted that the acceleration field is generated by extracting the instantaneous acceleration amplitude of each event from the microseismic event location results and then spatially interpolating it to produce the full-space acceleration field. In sparse regions of microseismic events or in static analysis scenarios, a quasi-static approximation is used. At this point, the momentum conservation residual degenerates into the static equilibrium equation residual. ; Calculated by automatic differentiation.

[0067] This is used to constrain the cubic law relationship between permeability and fracture aperture; Physical consistency regularization ,in, The predicted penetration rate output by the model; This corresponds to the crack aperture; , , These represent the effective porosity, shape factor, and tortuosity of the fracture network at that voxel, respectively. In a simplified scenario with relatively uniform fracture distribution and good connectivity, the physical consistency regularization term degenerates to: ; In scenarios where crack geometry parameters are missing or calculations are simplified, it can be further degraded to ; The Adam optimizer and conventional initialization strategy were used for training. A warm-up-cosine annealing learning rate strategy was adopted, and an early stopping mechanism was introduced to prevent overfitting.

[0068] For example, training employs a two-stage progressive strategy: Phase 1 (Data-driven pre-training, epochs 1-50): Setup Activate only and ; Phase 2 (Physical Constraint Reinforcement Training, epochs 51-200): Incrementing weights according to an adaptive weighting formula The physical conservation residual constraint is gradually activated. The mean physical residual is calculated on the validation set every 10 rounds. If this happens, the early stopping mechanism is triggered, rolling back to the most recent checkpoint and reducing the learning rate. The improved experimental results are shown in Table 1.

[0069] Table 1 Comparison of Disaster Risk Prediction Performance The risk prediction accuracy of this scheme reaches 96.8%, which is 24.3% higher than the traditional empirical formula method and 12.5% ​​higher than the single deep learning model. The physical consistency error is reduced to below 1.0%, and the physical reliability PCI reaches 0.96, indicating that the inverted fracture parameters and the predicted disaster risks strictly conform to the laws of rock mass fracture mechanics. Moreover, the accuracy remains above 92% in the generalization test of new mines (not involved in training), which is significantly better than the 78% of the comparison method.

[0070] This scheme designed three groups of ablation experiments, with the control groups being: A: the basic model (trained using only a static five-dimensional background tensor); B: introducing a dynamic temporal feature inversion module; C: introducing a physical information embedding and gating fusion mechanism; and D: the present scheme. Specific data are shown in Table 2.

[0071] Table 2 Ablation Experiment After introducing the physical information embedding and gating fusion mechanism (Group C), the accuracy improved by 3.6%, and the physical consistency error decreased by 33.9%, indicating that the explicit embedding of physical priors effectively improved the model's understanding of fracture mechanical behavior. After adding physical conservation residual loss and adaptive weight scheduling (Group D), the physical consistency error further decreased to below 1.0%, verifying the key role of physical constraints in ensuring that the prediction results conform to the laws of rock mechanics.

[0072] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of the invention. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. It should be understood that in the various embodiments of the invention, the sequence numbers of the above-described processes do not imply a sequential order of execution; the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the invention. The sequence numbers of the above-described embodiments of the invention are merely descriptive and do not represent the superiority or inferiority of the embodiments.

[0073] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0074] In the several embodiments provided by this invention, it should be understood that the disclosed methods can be implemented in other ways. The methods disclosed in the several method embodiments provided by this invention can be arbitrarily combined without conflict to obtain new method embodiments. The features disclosed in the several method embodiments provided by this invention can be arbitrarily combined without conflict to obtain new method embodiments.

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

Claims

1. A method for monitoring the evolution of fracture fields in layered mining of extra-thick coal seams under rockburst, characterized in that, include: Acquire microseismic monitoring point data and stress data, process the microseismic monitoring point data and stress data to obtain the five-dimensional crack-stress tensor; Based on real-time microseismic wave velocity data and real-time fiber optic strain data, a stress correction coefficient is introduced. The wave velocity-fracture aperture is dynamically adjusted using the stress correction coefficient to obtain the dynamic fracture aperture. Acquire fiber stress data, calculate the fiber stress data to obtain real-time stress; combine the dynamic fracture aperture, real-time stress and equivalent permeability to reconstruct the time-series fracture state characteristic sequence. A mining disaster risk prediction model is constructed and trained. The five-dimensional fracture-stress tensor and the temporal fracture state feature sequence are input into the mining disaster risk prediction model. After coupling and extrapolation by the model, the four-dimensional risk probability distribution tensor and physical reliability index at future time are output. The mining disaster risk prediction model sets up a static background encoding branch and a dynamic damage encoding branch. The static background encoding branch extracts the spatial features of the five-dimensional fracture-stress tensor, and the dynamic damage encoding branch extracts the temporal evolution features of the temporal fracture state feature sequence. After the two features are encoded by the physical information embedding layer to encode the fracture geometry-seepage constitutive relationship, the data-driven features and physical constraint features are adaptively fused through a physical consistency gating network, and then the prediction result is output after correction by a physical correction head. The four-dimensional risk probability distribution tensor and the physical reliability index are subjected to multi-channel separation processing to extract the potential disaster bodies in each channel. The volume change sequence of the potential disaster bodies is tracked to calculate the interlocking force index and the expected penetration time of the disaster crack. The disaster risk level is classified according to the interlocking force index, the expected penetration time of the disaster crack and the risk probability of each channel. Finally, the disaster risk level, the dominant disaster type and the expected penetration time of the disaster crack are output.

2. The method according to claim 1, characterized in that, The process of processing the microseismic monitoring point data and stress data to obtain the five-dimensional fracture-stress tensor includes: The inverse distance weighted interpolation algorithm is used to interpolate the microseismic monitoring point data and stress data to a three-dimensional grid space to construct a continuous physical field in the entire space. The five key physical parameters, namely fracture density, fracture opening, stress concentration factor, fracture attitude tensor, and rock mass splitting index, are calculated based on the continuous physical field in the entire space. The five key physical parameters are stacked in three-dimensional space to construct a five-dimensional fracture-stress tensor.

3. The method according to claim 1, characterized in that, The stress correction coefficient, introduced based on real-time microseismic wave velocity data and real-time fiber optic strain data, includes: Based on the reference stress and combined with the lithological stress sensitivity coefficient, the stress correction coefficient is constructed according to the deviation of the real-time principal stress from the initial background stress. The formula for the stress correction factor is as follows: in, This is the stress correction factor; This is the empirical benchmark coefficient; For real-time principal stress; The initial background stress; This is the lithological stress sensitivity coefficient.

4. The method according to claim 1, characterized in that, The step of dynamically adjusting the wave velocity-fracture aperture using the stress correction coefficient to obtain the dynamic fracture aperture includes: Obtain the initial fracture aperture, the theoretical wave velocity of the intact rock mass, and the real-time monitored wave velocity; The initial fracture aperture is multiplied by the square of the ratio of the theoretical wave velocity of the intact rock mass to the real-time monitored wave velocity, and then multiplied by the stress correction coefficient to calculate the dynamic fracture aperture.

5. The method according to claim 4, characterized in that, The time-series fracture state characteristic sequence is reconstructed by combining the dynamic fracture aperture, real-time stress, and equivalent permeability, including: Based on the cubic law, the dynamic fracture aperture is converted into equivalent permeability; The fracture type is determined based on the dynamic fracture aperture, real-time stress, and equivalent permeability. If the real-time stress is greater than or equal to the low stress threshold and less than or equal to the high stress threshold, or the equivalent permeability is less than or equal to the crack density threshold and the real-time stress is greater than or equal to the low stress threshold, or the equivalent permeability is greater than or equal to the crack density threshold and the real-time stress is less than or equal to the high stress threshold, then the crack is judged to be a transitional crack. If the fracture is determined to be a transitional fracture, microseismic data and stress data of the sub-coal seam are collected for T consecutive time steps; the inversion results corresponding to each time step are reconstructed in time series to obtain the time series fracture state characteristic sequence.

6. The method according to claim 1, characterized in that, The loss function used in the training process of the mining disaster risk prediction model includes weighted cross-entropy classification loss, spatial smoothing loss, physical conservation loss, and physical consistency regularization term. Specifically, the weighted cross-entropy classification loss assigns different weights to different fracture categories; the spatial smoothing loss constrains the rate of change of risk probability between adjacent voxels; the physical conservation loss includes mass conservation error terms and momentum conservation error terms; and the physical consistency regularization term constrains the model's predicted equivalent permeability to maintain a cubic law relationship with fracture aperture. An adaptive scheduling strategy is adopted for the weight coefficients of the physical conservation loss. In the early stage of training, data fitting is the main focus, and physical constraints are gradually strengthened in the later stages. The weight coefficients of the adaptive scheduling strategy increase with the training rounds.

7. The method according to claim 1, characterized in that, The static background encoding branch uses a three-dimensional convolutional neural network, taking the five-dimensional crack-stress tensor as input, to extract spatial features and output a static feature map; the dynamic damage encoding branch uses a hybrid architecture of a one-dimensional convolutional neural network and a Transformer, taking the temporal crack state feature sequence as input, to extract the temporal evolution law and output a temporal feature vector. The five-dimensional fracture-stress tensor and the temporal fracture state feature sequence are input into the mining disaster risk prediction model. After coupling and extrapolating through the model, the four-dimensional risk probability distribution tensor and physical reliability index for future moments are output, including: The temporal feature vector is mapped to the channel dimension via a spatiotemporal mapping module and broadcast to expand it to the same spatial resolution as the static feature map to generate a temporal feature map; The static feature map and the temporal feature map are concatenated along the channel dimension and then compressed using three-dimensional convolution to obtain the attention fusion feature. The crack geometry and seepage physical constraints are incorporated through a physical information embedding layer, and then a physical consistency gating network is used to adaptively fuse the physical feature map and the attention fusion feature. After the fused features are corrected by the physical correction head, they are mapped back to the original spatial resolution through a 3D transposed convolutional network, outputting a 4D risk probability distribution tensor and a physical credibility index for future time moments. The four-dimensional risk probability distribution tensor includes four dimensions: safety probability, rockburst probability, coal seam spontaneous combustion probability, and rockburst-spontaneous combustion coupling probability.

8. The method according to claim 1, characterized in that, The process of performing multi-channel separation processing on the four-dimensional risk probability distribution tensor and the physical credibility index to extract potential disaster bodies for each channel includes: By combining the four-dimensional risk probability distribution tensor with the physical credibility index, multi-channel separation processing is performed. A connected component analysis algorithm is used to extract continuous regions in each channel whose probability exceeds a preset threshold, and these continuous regions are identified as potential disaster bodies.

9. The method according to claim 1, characterized in that, The process of tracking the volume change sequence of the potential hazard body and calculating the interlocking force index and the estimated penetration time of the hazard fracture includes: The prediction results of the same potential disaster body at different times are dynamically tracked to obtain the volume change sequence of the potential disaster body. The volume change sequence is nonlinearly fitted to characterize the crack propagation trend. The predicted breakthrough time of the disaster crack is obtained based on the fitting results. Obtain the overlapping volume, union volume, and vertical distance between the centroids of the potential hazard bodies in the upper and lower layers; obtain the standard interlayer spacing and attenuation coefficient; first calculate the ratio of the overlapping volume to the union volume, then combine the difference between the vertical distance and the standard interlayer spacing and the attenuation coefficient to perform exponential normalization correction, and comprehensively calculate the interlocking force index.

10. The method according to claim 7, characterized in that, The process involves classifying disaster risk levels based on the interlocking force index, the estimated time for the disaster fissure to penetrate, and the risk probability of each channel. The final output includes the disaster risk level, the dominant disaster type, and the estimated time for the disaster fissure to penetrate, including: When the impact-spontaneous combustion coupling probability is greater than the first probability threshold and the biting force index is greater than the first force threshold, it is classified as a level one disaster risk; When the impact-spontaneous combustion coupling probability is less than the second probability threshold, the rockburst probability is greater than the first probability threshold, or the coal seam spontaneous combustion probability is greater than the first probability threshold, and the interlocking force index is less than the first force threshold, it is classified as a level two disaster risk. When the safety probability, rockburst probability, coal seam spontaneous combustion probability, and rockburst-spontaneous combustion coupling probability are all greater than the second probability threshold and less than the first probability threshold, or when the biting force index is less than the second force threshold, the risk is classified as Level III disaster risk.

Images