Analysis method for fracture induced shock mechanism of giant thick interbed sandstone roof

Abstract

The invention relates to the technical field of mining engineering and rock mechanics, and discloses a huge thick interbed sandstone roof fracture induced shock mechanism analysis method, which comprises the following steps: firstly, obtaining real-time microseismic data, and quantitatively constructing a four-dimensional space-time damage field representing rock mass internal damage evolution according to the real-time microseismic data; establishing an interlayer time-lag slip model driven by the damage field to realize dynamic coupling of a damage state and mechanical parameters; then, the model is adopted for dynamic simulation, potential fracture is identified through a multi-stage early warning mechanism based on a rate index and an energy dynamic instability criterion, and graded early warning is output; and finally, introducing a closed-loop feedback correction mechanism, and iteratively correcting the model through an optimization algorithm by using the deviation between an actual seismic source inversion parameter and a prediction result. According to the method, damage evolution and mechanical behaviors are closely combined, the self-adaptive correction capacity is achieved, a shock inducing physical mechanism can be revealed, and the accuracy and reliability of analysis and prediction are improved.

Description

Technical Field

[0001] This invention relates to the fields of mining engineering and rock mechanics, specifically to an analytical method for the fracture-induced seismic mechanism of a massive interbedded sandstone roof. Background Technology

[0002] In deep coal mining, the stability of the overlying thick interbedded sandstone roof is crucial for ensuring safe production. Under mining-induced stress disturbances, the internal damage accumulation, interlayer slippage, and macroscopic fracture processes of this type of roof are interconnected, often inducing high-energy mine seismic events, posing a serious threat to personnel and equipment. Therefore, accurate analysis and early warning of the fracture-induced seismic mechanisms of this type of roof have been a long-standing technical challenge in the field of mining engineering.

[0003] Currently, the analysis methods for roof fracture-induced seismic events mainly rely on numerical simulation and microseismic monitoring. Numerical simulation methods, such as the finite element method or discrete element method, typically construct physical models based on parameters obtained from laboratory rock mechanics tests to predict stress distribution and deformation caused by mining. However, these methods have limitations in application. Rock mechanics parameters in the model, such as elastic modulus, cohesion, and friction angle, are set during model initialization and usually remain unchanged throughout the simulation. This static parameter setting method cannot reflect the actual physical process of mechanical property deterioration caused by the continuous initiation, expansion, and convergence of microcracks within the rock mass under continuous mining influence. Therefore, existing models have a biased physical basis when describing the continuous evolution process from micro-damage accumulation to macro-dynamic instability, leading to discrepancies between simulation results and engineering reality.

[0004] Meanwhile, existing early warning mechanisms largely rely on threshold judgments based on specific indicators from numerical simulation results or microseismic monitoring data. For example, an early warning is triggered when the stress concentration in a region calculated by the model exceeds a preset intensity value, or when abnormal statistical patterns appear in the monitored microseismic energy and frequency. These methods fail to deterministically link early warning criteria with the intrinsic physical mechanisms of the failure event. Early warnings are often issued not based on a direct judgment of a dynamic instability process, but rather on indirect statistical indicators or static intensity indicators. This results in insufficient quantification and clarity of physical meaning in early warning information, making it difficult to provide a phased risk assessment basis based on the evolution of physical processes.

[0005] Furthermore, existing analytical methods generally operate in an open-loop mode. That is, the model performs a one-time prediction calculation based on the initial input, and the analysis process ends once the prediction result is output. When a failure event actually occurs in the field, there is a lack of a standardized process closely integrated with the model to utilize this valuable field measurement information. The model cannot systematically correct and optimize itself based on the deviation between the prediction result and the actual event. This makes the model overly dependent on the accuracy of the initial parameters and unable to gradually improve its adaptability to specific geological environments and predictive accuracy by learning from actual events during use. Summary of the Invention

[0006] The technical problem to be solved by this invention is that the existing technology is difficult to accurately and dynamically describe the entire process of the accumulation of internal micro-damage to macro-fracture and instability of a thick interbedded sandstone roof under the influence of mining. There is a lack of an analysis and early warning method that can tightly couple real-time monitoring data with physical models and has closed-loop feedback correction capabilities.

[0007] To address the aforementioned technical problems, this invention provides an analytical method for the seismic induction mechanism of fracture in a massive interbedded sandstone roof, the method comprising the following steps: S1. Data Acquisition and Monitoring A three-dimensional geological model of the thick, interbedded sandstone roof within the target area was obtained. This model includes geological structural information such as strata distribution, thickness, occurrence, and faults. Simultaneously, a microseismic monitoring system was deployed in the mining-affected area to acquire and record microseismic event data in real time. This microseismic data includes information such as the event's time of occurrence, three-dimensional spatial coordinates, and energy.

[0008] S2, Spatiotemporal Damage Field Construction Based on the microseismic data monitored in real time during step S1, a four-dimensional (three-dimensional space + one-dimensional time) spatiotemporal damage field is constructed. Used to quantitatively characterize any location inside the massive interbedded sandstone roof. At any time The damage state. In one specific embodiment, the spatiotemporal damage field It is based on the event spatial density of the microseismic data. and accumulated energy It is constructed using quantization. A specific computational model can be represented as: ; in: In order to be in time Spatial density of microseismic events within the location neighborhood; In order to be in time Accumulated microseismic energy within the location neighborhood; The coupling coefficient is a material-related parameter. Its physical meaning is to characterize the degree of material failure caused by a unit damage factor (defined by the product of density and energy). This parameter can be determined through indoor experiments or field inversion analysis.

[0009] As a damage factor, this product term incorporates the frequency of microseismic events (by...). (characterization) and the average intensity of the event (by) (Characteristics), serving as a comprehensive indicator for measuring cumulative damage intensity; This is the power of the exponent term, and it is negative, indicating that the absolute value of the term increases as the damage factor increases; This is an exponential term with values ​​in the interval (0,1]. Physically, it can be understood as the material's... The "integrity" or "residual strength factor" at the location. As the damage factor increases, the value of this term decreases from 1 to 0; The value is obtained through The calculations show that "integrity" is transformed into a damage variable. The range of values ​​is 0 indicates that the material is undamaged, and the closer the value is to 1, the higher the degree of material damage.

[0010] S3. Establishment of a damage-driven physical model Based on the spatiotemporal damage field constructed in step S2, a damage-driven interlayer time-delay slip model for the thick interbedded sandstone roof is established. A core feature of this model is that its key mechanical parameters are not constant values, but rather determined by the spatiotemporal damage field. Determined dynamically.

[0011] Preferably, the interlayer interface constitutive relation of the damage-driven interlayer time-delay slip model is characterized by the rate-state friction law. This law describes the interfacial friction coefficient. With slip rate and state variables The functional relationship between them can be expressed by a set of typical mathematical expressions: ; ; in: The coefficient of friction is the dynamic friction coefficient. This represents the current scrolling speed of the interface. State variables representing the contact state of the interface, in seconds; For reference slip rate The reference friction coefficient is as follows; This is a direct rate effect parameter, describing the effect of the friction coefficient changing instantaneously with the slip rate. This is the direct rate effect term; This is the state evolution effect term.

[0012] State variables Net rate of change; These are the slip weakening terms or contact update terms in the state evolution equation; The evolution effect parameter describes the effect of the friction coefficient changing with the evolution of the interface state; The critical slip distance is the characteristic slip length required for the evolution of state variables; The reference slip rate.

[0013] In one specific embodiment, the parameters in the rate-state friction law , and All are set as damage variables in the spatiotemporal damage field. The function is represented as: ; ; ; in: These represent the initial direct rate effect parameters, initial evolution effect parameters, and initial critical slip distance of the rock under undamaged conditions. For preset parameters describing the changes in damage variables The damage function of evolution.

[0014] In another specific embodiment, the key mechanical parameters of the model also include the elastic parameters of each rock layer within the thick interbedded sandstone top plate, such as the elastic modulus. The elastic modulus is also dynamically determined by the spatiotemporal damage field: ; in: The elastic modulus of the calculation element at the current time step; The damage variable is the current location of the potential fracture surface, read from the spatiotemporal damage field. This indicates that the elastic modulus is a damage variable. The function; This is the initial elastic modulus of the rock layer in its undamaged state, determined during the model initialization step. Let be a dimensionless damage evolution function, the value of which varies with . It decreases as it increases.

[0015] S4, Dynamic Simulation and Early Warning Numerical calculations and dynamic simulations were performed using the damage-driven interlayer time-delay slip model established in step S3. The purpose of the dynamic simulation was to identify potential future fracturing events and analyze the fracturing-induced seismic mechanism of the thick interlayered sandstone roof. Based on the simulation results and according to the preset multi-level dynamic fracturing criteria, graded early warning information for the potential fracturing events was output.

[0016] Preferably, the multi-level dynamic fracture criterion includes a first-level early warning criterion. The determination condition for the first-level early warning criterion is: during the dynamic simulation process, the interface slip rate in the monitoring model is... Damage evolution rate or interfacial energy release rate When any monitoring indicator exceeds its preset threshold, a Level 1 warning is triggered.

[0017] Further preferably, the multi-level dynamic failure criterion also includes a secondary early warning criterion. The secondary early warning criterion is determined based on the energy dynamic instability criterion. The specific steps for the determination are: calculating the system's energy release rate. and with the spatiotemporal damage field Determined dynamic fracture toughness Compare. When the condition is met. When the system is in a near-instability state, a level-two warning is triggered.

[0018] In one specific embodiment, the dynamic fracture toughness The value of is determined based on a damage-coupled fracture toughness model. This model links the initial fracture toughness of the material to the current damage state; a specific model can be expressed as: ; in: Dynamic fracture toughness, which takes into account the effects of damage, is physically defined as the material's performance under its current damage state. The ability to resist fracture propagation; For materials in an undamaged state ( The initial fracture toughness under the condition is determined by indoor experiments in the model initialization step; The damage variable is the current location of the potential fracture surface, read from the spatiotemporal damage field. The material integrity factor represents the proportion of the material's effective load-bearing portion lost due to damage. The damage effect index is a material-related index, a dimensionless parameter that determines the rate at which fracture toughness decreases with increasing damage. This is the function for reducing fracture toughness due to damage, and its value is within the range.

[0019] S5. Feedback Correction: After the potential rupture event actually occurs and becomes an observable rupture event, firstly, high-precision source inversion is performed on the actual source parameters of the rupture event to obtain information such as the actual rupture process and energy release. Then, the inversion results are quantitatively compared with the prediction results of the dynamic simulation for the event. Finally, based on the deviation generated by the comparison, feedback correction is performed on the parameters in the damage-driven interlayer time-delay slip model to improve the accuracy of subsequent predictions by the model.

[0020] In one specific embodiment, the feedback correction step specifically includes correcting the damage variables in the spatiotemporal damage field. With respect to the parameters in the aforementioned rate-state friction law ( The functional relationship between () and () is that the relationship between the functions () The form or coefficients are optimized and adjusted.

[0021] Preferably, before executing step S1, a model initialization step is included: obtaining the initial mechanical parameters of the thick interbedded sandstone roof in an undamaged state by conducting indoor mechanical loading experiments (e.g., triaxial compression experiments) on rock samples retrieved from the field engineering area, including but not limited to... The initial mechanical parameters are used for initializing the damage-driven interlaminar time-delay slip model.

[0022] This invention provides an analytical method for the seismic induction mechanism of fracture in the roof of a massive interbedded sandstone formation. It offers the following advantages: 1. This invention constructs a spatiotemporal damage field using real-time monitored microseismic data, and uses the spatiotemporal damage field to dynamically determine the key mechanical parameters in the damage-driven interlayer time-delay slip model. This makes the parameters of the physical model no longer static or empirically set, but able to reflect the non-uniform evolution of the internal damage state of the rock mass caused by mining. Thus, when describing the complete physical process from micro-damage accumulation to macro-fracture, the model has time-varying characteristics corresponding to the actual physical state.

[0023] 2. This invention uses numerical dynamic simulation to identify potential fracture events and sets up multi-level dynamic fracture criteria based on physical quantities, such as a primary criterion based on interface slip rate and a secondary criterion based on energy dynamic instability criteria. The simulated physical process output (such as energy release rate) is directly associated with explicit, graded criterion thresholds (such as dynamic fracture toughness), establishing a deterministic output path from physical model to graded early warning information, providing quantitative and non-fuzzy decision-making basis for risk assessment.

[0024] 3. This invention introduces high-precision source inversion, result comparison, and model feedback correction steps after a potential failure event actually occurs. This closed-loop correction mechanism uses actual failure event data to calibrate the functional relationship between damage variables and mechanical parameters in the model, enabling the model to iteratively optimize by comparing with actual field conditions. This reduces the model's long-term dependence on initial parameter settings and allows it to gradually adapt to the response characteristics of specific geological conditions. Attached Figure Description

[0025] Figure 1 A flowchart illustrating the analytical method for the seismic induction mechanism of fracture in a massive interlayered sandstone roof, according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the functional modules of an analysis system according to an embodiment of the present invention; Figure 3 This is a schematic diagram of the model initialization and data acquisition process according to an embodiment of the present invention; Figure 4 This is a schematic diagram of the spatiotemporal damage field construction process according to an embodiment of the present invention; Figure 5 This is a flowchart of a dynamic simulation and multi-level early warning mechanism according to an embodiment of the present invention; Figure 6 This is a flowchart of closed-loop feedback and model correction according to an embodiment of the present invention.

[0026] Among them, 10 is the data acquisition module; 20 is the spatiotemporal damage field calculation module; 30 is the physical model solving module; 40 is the early warning judgment module; and 50 is the feedback correction module. Detailed Implementation

[0027] See attached document Figure 1 , Figure 1 This is a flowchart illustrating an analysis method for the seismic induction mechanism of fracture in a massive interbedded sandstone roof according to an embodiment of the present invention. The present invention provides an analysis method for the seismic induction mechanism of fracture in a massive interbedded sandstone roof, which may include: acquiring a three-dimensional geological model and real-time microseismic data; constructing a spatiotemporal damage field characterizing the internal damage evolution of the rock mass based on the microseismic data; establishing a damage-driven interlayer time-delay slip model based on the spatiotemporal damage field, wherein the key mechanical parameters of the model are dynamically determined by the spatiotemporal damage field; performing dynamic simulation using the model to identify potential fracture events and output graded early warning information; and performing feedback correction on the model after a potential fracture event occurs.

[0028] See attached document Figure 2 , Figure 2 This is a functional module diagram of an analysis system according to an embodiment of the present invention. The system is used to execute the above-described analysis method and may include: a data acquisition module 10, a spatiotemporal damage field calculation module 20, a physical model solving module 30, an early warning judgment module 40, and a feedback correction module 50.

[0029] The data acquisition module 10 is used to acquire a three-dimensional geological model of the thick interbedded sandstone roof and to receive real-time microseismic data output by the microseismic monitoring system.

[0030] The spatiotemporal damage field calculation module 20, connected to the data acquisition module 10, is used to receive microseismic data and construct the spatiotemporal damage field. In one embodiment, the spatiotemporal damage field calculation module 20 calculates using the following formula: ; in: In order to be in time Spatial density of microseismic events within the location neighborhood; In order to be in time Accumulated microseismic energy within the location neighborhood; The coupling coefficient is a material-related parameter. Its physical meaning is to characterize the degree of material failure caused by a unit damage factor (defined by the product of density and energy). This parameter can be determined through indoor experiments or field inversion analysis.

[0031] As a damage factor, this product term incorporates the frequency of microseismic events (by...). (characterization) and the average intensity of the event (by) (Characteristics), serving as a comprehensive indicator for measuring cumulative damage intensity; This is the power of the exponent term, and it is negative, indicating that the absolute value of the term increases as the damage factor increases; This is an exponential term with values ​​in the interval (0,1]. Physically, it can be understood as the material's... The "integrity" or "residual strength factor" at the location. As the damage factor increases, the value of this term decreases from 1 to 0; The value is obtained through The calculations show that "integrity" is transformed into a damage variable. The range of values ​​is 0 indicates that the material is undamaged, and the closer the value is to 1, the higher the degree of material damage.

[0032] The physical model solving module 30 is connected to the spatiotemporal damage field calculation module 20, and is used to calculate the spatiotemporal damage field. A damage-driven interlayer time-delay slip model is established and solved. In one embodiment, the constitutive relation of the interlayer interface is characterized by the rate-state friction law: ; ; in, The coefficient of friction is the dynamic friction coefficient. This represents the current scrolling speed of the interface. State variables representing the contact state of the interface, in seconds; For reference slip rate The reference friction coefficient is as follows; This is a direct rate effect parameter, describing the effect of the friction coefficient changing instantaneously with the slip rate. This is the direct rate effect term; This is the state evolution effect term.

[0033] State variables Net rate of change; These are the slip weakening terms or contact update terms in the state evolution equation; The evolution effect parameter describes the effect of the friction coefficient changing with the evolution of the interface state; The critical slip distance is the characteristic slip length required for the evolution of state variables; The reference slip rate.

[0034] The early warning determination module 40, connected to the physical model solving module 30, receives the calculation results of the dynamic simulation and makes a determination based on multi-level dynamic fracture criteria. In one embodiment, when the early warning determination module 40 executes the second-level early warning criterion, it needs to calculate the dynamic fracture toughness. Its calculation is based on the damage-coupled fracture toughness model: ; in: Dynamic fracture toughness, which takes into account the effects of damage, is physically defined as the material's performance under its current damage state. The ability to resist fracture propagation; For materials in an undamaged state ( The initial fracture toughness under the condition is determined by indoor experiments in the model initialization step; The damage variable is the current location of the potential fracture surface, read from the spatiotemporal damage field. The material integrity factor represents the proportion of the material's effective load-bearing portion lost due to damage. The damage effect index is a material-related index, a dimensionless parameter that determines the rate at which fracture toughness decreases with increasing damage. This is the function for reducing fracture toughness due to damage, and its value is within the range.

[0035] The feedback correction module 50 is connected to the data acquisition module 10 and the physical model solving module 30. After a potential failure event occurs, the feedback correction module 50 receives the high-precision source inversion results of the event, compares them with the prediction results of the physical model solving module 30, and corrects the parameter function relationships in the physical model solving module 30 based on the comparison deviation.

[0036] See attached document Figure 1 and attached Figure 2 The analysis method of this invention is a closed-loop process that includes feedforward prediction and feedback correction. The process is executed collaboratively by the data acquisition module 10, the spatiotemporal damage field calculation module 20, the physical model solving module 30, the early warning judgment module 40, and the feedback correction module 50.

[0037] The method begins with data acquisition module 10. Data acquisition module 10 is responsible for acquiring static 3D geological model data and continuously receiving dynamic microseismic data streams from the field monitoring system. The 3D geological model provides the geometric boundaries and initial geological structure for subsequent physical model establishment. Microseismic data serves as a direct input signal characterizing damage activity within the rock mass.

[0038] The data acquisition module 10 transmits the real-time acquired microseismic data to the spatiotemporal damage field calculation module 20. The spatiotemporal damage field calculation module 20 processes the input microseismic data and calculates the damage variables based on the event spatial density and cumulative energy within a preset four-dimensional grid (three-dimensional space plus one-dimensional time). This allows for the construction of a dynamically evolving spatiotemporal damage field.

[0039] The spatiotemporal damage field calculation module 20 calculates the spatiotemporal damage field. The data is transmitted to the physical model solving module 30 in real time or periodically. The physical model solving module 30 internally constructs a damage-driven interlayer time-delay slip model. Upon receiving new spatiotemporal damage field data, the physical model solving module 30 utilizes the spatiotemporal damage field... The values ​​are updated to reflect the key mechanical parameters defined within the model, such as those in the rate-state friction law. and the elastic modulus of the rock mass .

[0040] After the parameters are updated, the physical model solving module 30 performs forward dynamic simulation calculations under external load conditions such as mining. The calculation results of the dynamic simulation, including but not limited to physical quantities such as interface slip rate and system energy release rate, are continuously transmitted to the early warning determination module 40. The early warning determination module 40 compares the received physical quantities with the thresholds of internally set multi-level dynamic failure criteria. When the trigger condition of any criterion is met, the early warning determination module 40 outputs the corresponding level of early warning information, which clearly points to one or more potential failure events identified by the dynamic simulation.

[0041] When a potential fracturing event that has been flagged as a warning actually occurs in the physical world, the method's feedback correction process is initiated. The data acquisition module 10 (or an external input interface) acquires high-precision source inversion results for the actual fracturing event. The inversion results are transmitted to the feedback correction module 50. Simultaneously, the physical model solving module 30 also transmits its previous predictions for the event to the feedback correction module 50.

[0042] The feedback correction module 50 quantifies the deviation between the prediction and inversion results and calculates the correction amount for the model parameters based on a preset optimization algorithm. This correction amount is fed back to the physical model solution module 30 to adjust the functional relationship between the damage variable and the force parameters. This step completes the closed-loop operation of the method, allowing the model to be iteratively optimized using actual event data.

[0043] See attached document Figure 3 , Figure 3 This is a schematic diagram of the model initialization and data acquisition process according to an embodiment of the present invention. The model initialization and data acquisition process is the starting stage of the analysis method of the present invention, providing basic data and initial model parameters for subsequent damage field construction and dynamic simulation.

[0044] Before executing the main steps of the method, model initialization is required. The purpose of this step is to obtain the impact of the thick interbedded sandstone roof on the unmined surface, i.e., the damage variables. The initial mechanical parameters are zero or close to zero. Specifically, rock samples were obtained from representative strata constituting the thick interbedded sandstone roof through in-situ core drilling. The rock samples were prepared into standard specimens for a series of laboratory mechanical loading experiments.

[0045] The mechanical loading experiments included: performing uniaxial and triaxial compression tests, and determining the initial elastic modulus of each rock layer in the undamaged state by recording the stress-strain curves of the specimens throughout the entire process. Compared to Poisson Fracture tests were performed on notched specimens (e.g., three-point bending beams) to determine the initial fracture toughness of each rock layer. Direct shear tests under different normal stresses, especially variable-rate shear tests, were performed on rock interfaces or pre-jointed surfaces. By monitoring the relationship between shear stress and slip rate, the initial parameters in the rate-state friction law, including the initial direct rate effect parameters, were determined. Initial evolution effect parameters and initial critical slip distance All acquired initial parameters were recorded and used as baseline inputs for the subsequent physics model solver module 30.

[0046] Step S1 of the main process involves the acquisition of field data. This step comprises two parallel parts: the construction of a 3D geological model and real-time monitoring of microseismic data. The construction of the 3D geological model is based on field geological exploration data, including borehole core logs, seismic reflection profiles, and geological mapping data. By interpreting the data, the spatial distribution, geometry, thickness variations of the main rock strata, and the location and attitude of faults, folds, and other structures are determined. Subsequently, geological modeling software or numerical simulation preprocessing software is used to transform the interpreted geological information into a discretized 3D numerical model. The model represents the spatial structure of the geological body in the form of a mesh (e.g., tetrahedral or hexahedral elements), and each mesh element is assigned initial mechanical parameters determined by the model initialization step.

[0047] Real-time monitoring of microseismic data is achieved through a microseismic monitoring system deployed within and around the mining-affected area. The system consists of multiple (usually no fewer than eight) microseismic sensors (e.g., three-component geophones or accelerometers) distributed at different locations on the surface or underground, a data acquisition unit, and a central data processing server. Each sensor continuously collects vibration signals from the formation and transmits the digital signals to the central data processing server via a wired or wireless data transmission network. The server runs an event triggering and location algorithm. When multiple sensors simultaneously receive vibration signals with energy exceeding the background noise threshold, a microseismic event is determined to have occurred, and the server calculates the event's occurrence time (accurate to milliseconds), three-dimensional spatial coordinates, magnitude, or energy level, and other source parameters. These parameters constitute a microseismic data record and are continuously and in real-time output to the spatiotemporal damage field calculation module 20.

[0048] See attached document Figure 4 , Figure 4 This is a schematic diagram of the spatiotemporal damage field construction process according to an embodiment of the present invention. The spatiotemporal damage field calculation module 20 executes this process, and its core task is to transform the microseismic event data received from the data acquisition module 10, which is discretely distributed in space and time, into a continuous, four-dimensional (three-dimensional space + one-dimensional time) spatiotemporal damage field. .

[0049] The physical basis for the construction process lies in the damage evolution of rock mass under stress, where the initiation, propagation, and connection of internal microcracks release energy in the form of microseismic events. Therefore, the density and energy magnitude of microseismic events within a specific spatiotemporal region directly characterize the damage state of the rock mass in that region. This embodiment employs a quantitative model that comprehensively reflects event density and energy to construct the spatiotemporal damage field.

[0050] Before the actual calculations, the analysis region is first spatially discretized, dividing it into a series of computational grid points or volume elements. The spatial coordinates of each grid point are... The calculation of the spatiotemporal damage field will be performed at each grid point. Up, step by step over time Evolution proceeds.

[0051] For any computational grid point and at any time The spatiotemporal damage field calculation module 20 first calculates the spatial density of microseismic events at that point. and weighted cumulative energy .

[0052] Spatial density of microseismic events The calculation employs a kernel density estimation method. Specifically, for the time up to the current moment... All that has happened Each microseismic event, (Its spatial coordinates are) For calculation points The density contribution is determined by a spatial kernel function. Determined. Spatial density. The calculation formula is: ; in: From the initial time to the current time The total number of microseismic events that occurred; The summation symbol represents the summation from the first to the last. The calculation results of all microseismic events are summed up; For the first Spatial coordinates of the epicenter of a microseismic event; This is a spatial vector pointing from the event location to the calculation point; The bandwidth parameter of the kernel function, also known as the influence radius, determines the spatial range of influence of a single microseismic event; The spatial kernel function is a function that varies with distance. The weight function that increases and decreases is used to calculate the first... Micro-seismic events at point The density contribution value at that location.

[0053] In a specific embodiment, the spatial kernel function The Gaussian kernel function can be used, and its expression in d-dimensional space (d=3 in this embodiment) is as follows: ; in Indicates the calculation point With the event The Euclidean distance between them. Bandwidth parameter. The value of can be determined based on the scale of the study area and the average spacing of microseismic events. For example, it can be 2 to 5 times the average nearest neighbor distance of the microseismic events.

[0054] Weighted cumulative energy The calculation is similar to spatial density; it is performed at the calculation point. At this point, the energy of all historical microseismic events is weighted and accumulated. The weights are also determined by the spatial kernel function. Confirmed. The calculation formula is: ; in: For the first The energy released by a microseismic event; Indicates the first Microseismic events affect the calculation point Weighted energy contribution at the location.

[0055] In calculation and Subsequently, the spatiotemporal damage field calculation module 20 calculates the damage variables using the following damage evolution formula. ; in: In order to be in time Spatial density of microseismic events within the location neighborhood; In order to be in time Accumulated microseismic energy within the location neighborhood; The coupling coefficient is a material-related parameter. Its physical meaning is to characterize the degree of material failure caused by a unit damage factor (defined by the product of density and energy). This parameter can be determined through indoor experiments or field inversion analysis.

[0056] As a damage factor, this product term incorporates the frequency of microseismic events (by...). (characterization) and the average intensity of the event (by) (Characteristics), serving as a comprehensive indicator for measuring cumulative damage intensity; This is the power of the exponent term, and it is negative, indicating that the absolute value of the term increases as the damage factor increases; This is an exponential term with values ​​in the interval (0,1]. Physically, it can be understood as the material's... The "integrity" or "residual strength factor" at the location. As the damage factor increases, the value of this term decreases from 1 to 0; The value is obtained through The calculations show that "integrity" is transformed into a damage variable. The range of values ​​is 0 indicates that the material is undamaged, and the closer the value is to 1, the higher the degree of material damage.

[0057] After completing the calculation, the spatiotemporal damage field calculation module 20 outputs a four-dimensional array. The parameters are then transferred to the physical model solving module 30 for dynamic updates of the model's mechanical parameters.

[0058] The physical model solving module 30 executes this process; its function is to receive the spatiotemporal damage field output by the spatiotemporal damage field calculation module 20. This is used as the driving field to update the key mechanical parameters in a numerical physics model in real time, thereby establishing a damage-driven interlayer time-delay slip model that can reflect the dynamic evolution of damage state inside the rock mass.

[0059] In one embodiment, the physical model is implemented on a numerical computing platform (e.g., ABAQUS, FLAC3D, PFC3D, etc.) based on a numerical computation method, such as the finite element method, finite difference method, or discrete element method.

[0060] In practical implementation, the secondary development interface provided by the numerical calculation platform can be utilized. For example, when using the finite element software ABAQUS, this can be achieved by writing a User Material Subroutine (UMAT) or a User Element Subroutine (UEL). At the beginning of each incremental step, the subroutine reads the damage variable of the corresponding element from an external file. The latest value is then used to calculate the elastic modulus used within this incremental step according to the coupling relationship of the present invention. and interfacial friction parameters This allows for a tight coupling between the damage field and the mechanical behavior of the model.

[0061] The three-dimensional geological model is discretized into a series of computational units. Spatiotemporal damage field. This is applied to these computational units. Specifically, at the start of each analysis time step, each computational unit in the numerical model is activated from the spatiotemporal damage field via user-defined subroutines or built-in scripting language functions. Read the current damage variable at its corresponding spatial location. The value of the damage variable is then determined. Subsequently, the physical model solving module 30, based on the preset damage-mechanical parameter coupling relationship, utilizes the read damage variables... To update the mechanical parameters of the unit.

[0062] The damage-mechanical parameter coupling relationship includes at least the updating of the rock layer elastic parameters and the interlayer interface friction parameters.

[0063] For the elastic parameters of rock strata, the elastic modulus is used. For example, its damage-driven update relationship is represented as: ; in: The elastic modulus of the calculation element at the current time step; The damage variable is the current location of the potential fracture surface, read from the spatiotemporal damage field. This indicates that the elastic modulus is a damage variable. The function; This is the initial elastic modulus of the rock layer in its undamaged state, determined during the model initialization step. Let be a dimensionless damage evolution function, the value of which varies with . It decreases as it increases.

[0064] In one embodiment, the function can take a linear form: In another embodiment, an exponential form may be used: ,in It is the material-related damage influence coefficient.

[0065] For the friction parameters of the interlayer interface, this embodiment uses the rate-state friction law for characterization. The mathematical expression of the law is: ; ; in, The coefficient of friction is the dynamic friction coefficient. This represents the current scrolling speed of the interface. State variables representing the contact state of the interface, in seconds; For reference slip rate The reference friction coefficient is as follows; This is a direct rate effect parameter, describing the effect of the friction coefficient changing instantaneously with the slip rate. This is the direct rate effect term; This is the state evolution effect term.

[0066] State variables Net rate of change; These are the slip weakening terms or contact update terms in the state evolution equation; The evolution effect parameter describes the effect of the friction coefficient changing with the evolution of the interface state; The critical slip distance is the characteristic slip length required for the evolution of state variables; The reference slip rate.

[0067] Among them, key parameters (Direct rate effect parameter) (Evolutionary effect parameters) and (Critical slip distance) was set as a damage variable. Functions: ; ; ; in: These represent the initial direct rate effect parameters, initial evolution effect parameters, and initial critical slip distance of the rock under undamaged conditions. For preset parameters describing the changes in damage variables The damage function of evolution.

[0068] In one specific embodiment, the function may take the following form: ,in For parameters The damage impact index.

[0069] ,in For parameters The damage impact index.

[0070] ,in For parameters The damage influence coefficient, in this form, indicates that the critical slip distance increases linearly with increasing damage.

[0071] At each time step, the physical model solver 30 performs the aforementioned parameter update calculations for each computational and interface element in the model. After all parameter updates are completed, the governing equations are solved within that time step to calculate the increments in stress, strain, and displacement. This process ensures that the mechanical behavior of the model is based on the current actual damage state of the rock mass, reflected by microseismic activity, rather than on its initial, undamaged state.

[0072] See attached document Figure 5 , Figure 5 This is a flowchart of a dynamic simulation and multi-level early warning mechanism according to an embodiment of the present invention. The process is executed collaboratively by the physical model solving module 30 and the early warning determination module 40. After the physical model solving module 30 completes the update of the mechanical parameters based on the spatiotemporal damage field, it enters the dynamic simulation and early warning stage.

[0073] Dynamic simulation is a forward time-stepping calculation process. The physical model solving module 30 applies boundary conditions or load increments caused by mining activities to the established damage-driven interlayer time-delay slip model. Subsequently, the physical model solving module 30 calculates in a small time step... Iterative solutions are performed to calculate the changes in physical quantities such as stress, strain, displacement, and interface slip rate at various points within the model under external loads. This process is repeated continuously to simulate the continuous dynamic response of the thick interbedded sandstone roof under mining influence. After the calculation is completed at each time step, the physical model solution module 30 outputs the calculation results of a set of key physical quantities to the early warning judgment module 40.

[0074] The early warning judgment module 40 receives dynamic simulation results and identifies and issues early warnings for potential failure events in real time based on preset multi-level dynamic failure criteria. The criteria include at least a primary early warning criterion and a secondary early warning criterion. The primary early warning criterion is based on the rate index of the monitored physical process. The early warning judgment module 40 continuously monitors the interface slip rate in key areas of the model (e.g., in front of the goaf or stress concentration zone). Damage evolution rate or interfacial energy release rate The module internally sets its own warning thresholds for these indicators. When the value of any metric exceeds its corresponding preset threshold, for example... The early warning judgment module 40 immediately generates and outputs a Level 1 early warning message. A Level 1 early warning indicates that the deformation or damage process of the rock mass is accelerating.

[0075] The secondary early warning criterion is based on the energy dynamic instability criterion and is used to determine whether the system is nearing macroscopic instability. At each time step of the dynamic simulation, the early warning judgment module 40 performs the following calculations and comparisons for potential rupture surfaces in the model: First, calculate the energy release rate of the system. Energy release rate Defined as the decrease in the elastic energy of a system when the potential fracture surface expands by a unit area, it characterizes the energy supply rate that drives fracture expansion.

[0076] Secondly, calculate the dynamic fracture toughness determined by the current damage state. Dynamic fracture toughness is a measure of a material's ability to resist crack propagation. The early warning judgment module 40 calculates this based on the damage-coupled fracture toughness model. ; in: Dynamic fracture toughness, which takes into account the effects of damage, is physically defined as the material's fracture toughness under its current damage state. The ability to resist fracture propagation; For materials in an undamaged state ( The initial fracture toughness under the condition is determined by indoor experiments in the model initialization step; The damage variable is the current location of the potential fracture surface, read from the spatiotemporal damage field. The material integrity factor represents the proportion of the material's effective load-bearing portion lost due to damage. The damage effect index is a material-related index, a dimensionless parameter that determines the rate at which fracture toughness decreases with increasing damage. This is the function for reducing fracture toughness due to damage, and its value is within the range.

[0077] Finally, the early warning judgment module 40 compares the energy release rate. With dynamic fracture toughness (or its equivalent energy form). When the instability condition is met, for example... (in (To convert the effective elastic modulus, toughness into energy form), the early warning determination module 40 determines that the system is nearing instability at this location and generates a secondary early warning message. The secondary early warning indicates that a high-energy macroscopic failure event is about to occur.

[0078] See attached document Figure 6 , Figure 6This is a flowchart of closed-loop feedback and model correction according to an embodiment of the present invention. The process forms a closed loop in the analysis method of the present invention, the purpose of which is to iteratively optimize the damage-driven interlaminar time-delay slip model using actual fracture event data, so that its prediction results gradually approach engineering reality. The process is primarily executed by the feedback correction module 50.

[0079] When a rupture event with clearly defined source characteristics, actually monitored on-site, occurs, the closed-loop feedback correction process is triggered. First, complete waveform data of the event is acquired from the microseismic monitoring system. Using this waveform data, high-precision source inversion is performed. In one embodiment, the inversion employs the moment tensor inversion method to calculate a set of actual source parameters characterizing the physical properties of the rupture event. Including seismic moments , fracture area Mean slip and stress drop wait.

[0080] Subsequently, the feedback correction module 50 retrieves the prediction results of the previous dynamic simulation for the same event from the physical model solution module 30, namely a set of simulated source parameters. Including simulated seismic moments Simulated fracture area The feedback correction module 50 uses a preset objective function to quantify the prediction bias. To compare the two sets of parameters. In one embodiment, the objective function... Defined as the normalized weighted sum of squared residuals: ; in: The total number of physical parameters involved in the comparison; The summation symbol represents summation over all... The calculation results of each physical parameter are accumulated; For the first The actual inversion values ​​of each physical parameter; For the first Model predictions for each physical parameter; For the first The absolute deviation of a physical parameter; For the first Normalized relative deviation of each physical parameter; For the first The weighting coefficients of each physical parameter are used to adjust the importance of different parameters in the objective function.

[0081] The core step of feedback correction is to adjust the adjustable parameters in the model using an optimization algorithm to improve the objective function. Minimize the value of . Adjustable parameters are variables that describe damage. With parameters of the rate-state friction law ( The coefficients in the functional relationship between (e.g., parameters) Let these adjustable parameters form a parameter vector. .

[0082] In one embodiment, the feedback correction module 50 employs gradient descent to perform optimization. Besides gradient descent, the optimization algorithm can also be other algorithms capable of solving such parametric optimization problems. For example, for the objective function... For non-convex models or models with multiple local minima, global optimization algorithms, such as simulated annealing or genetic algorithms, can be used. These algorithms, by introducing random perturbations or population evolution mechanisms, have the ability to escape local optima, which helps in finding better model parameter vectors. In the first In this iteration, the update rule for the parameter vector is: ; in: For the first The parameter vector at the next iteration; For the first The parameter vector updated in the next iteration; For the first The learning rate in each iteration controls the step size of parameter updates; For the objective function exist For parameter vectors The gradient of the objective function. This gradient vector indicates the direction in which the objective function value increases the fastest, so updating the parameters along its negative direction will decrease the objective function value.

[0083] Feedback correction module 50 repeats the above iterative calculation until the objective function is achieved. The values ​​converge to a preset tolerance range, or reach a preset maximum number of iterations. After the optimization process is complete, a set of optimal parameter vectors is obtained. The feedback correction module 50 will then use this optimal parameter vector. The feedback is sent to the physical model solving module 30 to update the damage-mechanical parameter coupling function. Subsequently, the physical model solving module 30 will use this set of functions, corrected for actual events, in all subsequent dynamic simulations to complete the closed-loop correction of the model.

[0084] The following will provide a detailed explanation of the entire process of the analysis method for the fracture-induced seismic mechanism of the massive interbedded sandstone roof provided by the present invention through a specific, non-limiting embodiment.

[0085] The application scenario of this embodiment is the 1102 working face of a deep coal mine. The working face is covered by a very thick interbedded sandstone roof, and the risk of rockburst induced by roof failure during mining is high.

[0086] Step 1: Model Initialization and Data Acquisition Before the working face is mined, model initialization is performed first. Core samples are drilled from key rock strata (e.g., fine sandstone and medium sandstone) in the roof of the 1102 working face and prepared into standard specimens. In the laboratory, triaxial compression tests and fracture toughness tests are performed on the specimens to obtain the initial mechanical parameters of the undamaged rock mass. For example, the initial elastic modulus of the fine sandstone is measured. GPa, initial fracture toughness 2.8MPa Through variable-rate direct shear experiments, the initial parameters of the rate-state friction law at the rock strata interface were determined as follows: m.

[0087] Simultaneously, data acquisition module 10 begins on-site data acquisition. Based on the mine's geological exploration report, a three-dimensional geological numerical model covering the 1102 working face and its surrounding area is established. A microseismic monitoring system consisting of 16 three-component geophones is deployed in the underground roadways and on the surface of this area. The system begins to record microseismic events triggered by the working face mining in real time and transmits the event's time, three-dimensional coordinates, and energy data stream to the spatiotemporal damage field calculation module 20 in real time.

[0088] Step 2: Quantitative Construction of Spatiotemporal Damage Field As the 1102 working face progressed as planned, the number of microseismic events recorded by the microseismic monitoring system gradually increased. The spatiotemporal damage field calculation module 20 continuously received this data. At a certain moment... The spatiotemporal damage field calculation module 20 targets each grid point in the numerical model. Using all microseismic events up to that point, calculate the spatial density of the microseismic events. and weighted cumulative energy .

[0089] Subsequently, the damage evolution formula was substituted into the spatiotemporal damage field calculation module 20: ; Wherein, coupling coefficient Preset The module uses this formula to calculate the value of each grid node in the model. The specific numerical values ​​of the damage variables on the surface. All these damage variable values ​​at discrete spatiotemporal nodes collectively constitute a dynamically evolving spatiotemporal damage field. In a specific calculation case, the damage variable value of a certain area of ​​the roof in front of the working face accumulates from near 0 when it was not mined to 0.4. Ultimately, this complete spatiotemporal damage field, existing in the form of a four-dimensional dataset (which can be uniformly represented mathematically as a function), is... The output is sent to the physical model solver module 30.

[0090] Step 3: Establishing a damage-driven physical model After receiving the damage field data, the physical model solving module 30 updates the model parameters before the start of the next dynamic simulation time step. For damage variables... In the region where it reaches 0.4, its elastic modulus is... according to The value is updated, and a new elastic modulus is calculated. GPa GPa.

[0091] Simultaneously, the friction parameters at the interlayer interface in this region are also updated. Assume the damage evolution function is... , .

[0092] The updated parameters are: ; ; m m; The parameters of other regions in the model are also updated accordingly based on their respective damage variable values.

[0093] Step 4: Dynamic Simulation and Multi-Level Early Warning Mechanism The physical model solver module 30 uses the updated parameters to continue simulating the next mining step. Simulation results show that in the region with the highest concentration of damage, the slip rate at the interlayer interface... The speed suddenly accelerated to 0.1 mm / s. The warning determination module 40 received this speed value and determined that it exceeded the preset first-level warning threshold. mm / s. The system then output a Level 1 warning message, indicating that "deformation is accelerating in the XX area of ​​the top plate in front of the 1102 working face, posing a risk of instability."

[0094] The Level 1 warning information may specifically include: warning level (Level 1), trigger time, trigger indicator (interface scrolling rate), current value (0.1 mm / s), threshold (0.08 mm / s), and the three-dimensional spatial coordinate range of the risk area.

[0095] The physical model solver module 30 continues the simulation. Calculations show that the energy release rate in this region... It is rapidly approaching its critical value. The early warning judgment module 40 simultaneously calculates the dynamic fracture toughness of this region. Assuming a damage influence index... The current dynamic fracture toughness is MPa When the simulated energy release rate When the instability criterion is met, the early warning judgment module 40 determines that a macroscopic fracture is about to occur in the area, and then outputs a level-two early warning message, which includes: "It is expected that a large-energy fracture event will occur in the XX area of ​​the roof in front of the 1102 working face within the next 12 hours, with an estimated magnitude of M1.8."

[0096] The specific information of a Level II warning may include: warning level (Level II), release time, expected occurrence time window (within the next 12 hours), three-dimensional spatial coordinate range of the potential rupture area, estimated energy or magnitude (M1.8), and warning confidence level (e.g., 85% based on historical data backtesting).

[0097] Step 5: Closed-loop feedback and model correction Ten hours after the warning was issued, the on-site microseismic monitoring system recorded a M1.7 mine tremor, with the epicenter location largely matching the area indicated in the secondary warning. Following the event, the feedback correction module 50 was activated. First, moment tensor inversion was performed using the waveform data of the M1.7 event to obtain a set of actual epicenter parameters. .

[0098] Feedback correction module 50 Compared with the parameters previously predicted by the model Compare and calculate the objective function. The value is 0.25. The feedback correction module 50 then initiates the gradient descent optimization algorithm to optimize the exponents (e.g., in the coupling relationship between damage and friction parameters) in the damage friction parameter coupling relationship. Adjustments were made. After several iterations, when the objective function... The optimization process terminates when the value decreases to 0.05. A set of corrected parameters is then obtained, for example, The value was adjusted from the initial setting of 0.5 to 0.62.

[0099] The corrected parameters were sent back to the physical model solving module 30, replacing the original parameters. From then on, the system will use this set of model parameters, which has been calibrated with actual events and better reflects the characteristics of the top plate of the working face, in subsequent analysis and early warning for the 1102 working face, thereby completing a complete closed-loop feedback correction.

Claims

1. An analytical method for the seismic induction mechanism of fracture in a massive interbedded sandstone roof, characterized in that, Includes the following steps: S1. Obtain a three-dimensional geological model of the thick interbedded sandstone roof and monitor the microseismic data of the mining-affected area in real time; S2. Based on the microseismic data, construct a spatiotemporal damage field characterizing the internal damage evolution of the thick interbedded sandstone roof. S3. Based on the spatiotemporal damage field, establish a damage-driven interlayer time-delay slip model for the massive interlayered sandstone roof, wherein the key mechanical parameters of the damage-driven interlayer time-delay slip model for the massive interlayered sandstone roof are dynamically determined by the spatiotemporal damage field. S4. Dynamic simulation is performed using the damage-driven interlayer time-delay slip model. The dynamic simulation is used to identify potential fracture events, analyze and reveal the fracture-induced seismic mechanism of the thick interlayered sandstone roof, and output graded early warning information for the potential fracture events based on a preset multi-level dynamic fracture criterion including at least one early warning level. S5. After the potential failure event actually occurs and forms a failure event, perform high-precision source inversion on the failure event, compare the inversion results with the prediction results of the dynamic simulation, and feed back to correct the damage-driven interlayer time-delay slip model.

2. The analytical method for the fracture-induced seismic mechanism of a super-thick interbedded sandstone roof as described in claim 1, characterized in that, In step S2, the spatiotemporal damage field is constructed by quantification based on the event spatial density and cumulative energy of the microseismic data.

3. The analytical method for the fracture-induced seismic mechanism of a massive interbedded sandstone roof as described in claim 1, characterized in that, In step S3, the constitutive relation of the damage-driven interlayer time-delay slip model of the thick interbedded sandstone roof is characterized by the rate-state friction law.

4. The analytical method for the seismic induction mechanism of fracture of a super-thick interbedded sandstone roof as described in claim 3, characterized in that, The direct rate effect parameter, evolution effect parameter, and critical slip distance in the rate-state friction law are all functions of the damage variables in the spatiotemporal damage field.

5. The analytical method for the fracture-induced seismic mechanism of a super-thick interbedded sandstone roof as described in claim 1, characterized in that, In step S3, the key mechanical parameters of the damage-driven interlayer time-delay slip model of the massive interbedded sandstone roof include the elastic parameters of each rock layer within the massive interbedded sandstone roof.

6. The analytical method for the fracture-induced seismic mechanism of a super-thick interbedded sandstone roof as described in claim 1, characterized in that, In step S4, the multi-level dynamic fracture criterion includes a first-level early warning criterion, and the determination of the first-level early warning criterion includes: The interface slip rate, damage evolution rate, and interface energy release rate in the damage-driven interlayer time-delay slip model are monitored. When any of these indicators exceeds a preset threshold, a first-level warning is triggered.

7. The analytical method for the seismic induction mechanism of fracture of a super-thick interbedded sandstone roof as described in claim 1, characterized in that, In step S4, the multi-level dynamic fracture criterion further includes a secondary early warning criterion. The secondary early warning criterion is based on an energy dynamic instability criterion, and the determination of the energy dynamic instability criterion includes: Compare the system energy release rate with the dynamic fracture toughness determined by the spatiotemporal damage field.

8. The analytical method for the fracture-induced seismic mechanism of a super-thick interbedded sandstone roof as described in claim 7, characterized in that, The value of the dynamic fracture toughness is determined based on the damage-coupled fracture toughness model, according to the initial fracture toughness of the thick interbedded sandstone roof and the damage variables in the spatiotemporal damage field.

9. The analytical method for the fracture-induced seismic mechanism of a super-thick interbedded sandstone roof as described in claim 4, characterized in that, In step S5, the feedback correction of the damage-driven interlaminar time-delay slip model specifically includes: The functional relationship between the damage variables in the spatiotemporal damage field and the parameters in the rate-state friction law is corrected.

10. The analytical method for the fracture-induced seismic mechanism of a massive interbedded sandstone roof as described in claim 1, characterized in that, Before performing step S1, the method further includes: conducting a mechanical loading experiment on rock samples taken from the site to obtain the initial mechanical parameters of the massive interbedded sandstone roof in an undamaged state. The initial mechanical parameters are used to initialize the damage-driven interlayer time-delay slip model.

Images