Microseismic monitoring inversion and intelligent anomaly recognition method for dynamic and static stress fields of coal and rock strata

Abstract

The present invention discloses a method for inversion of microseismic monitoring of dynamic and static stress fields of coal-rock formations and intelligent recognition of abnormalities. Firstly, using data such as source location and coordinates of microseismic sensors, wave velocity field inversion imaging is performed on coal-rock formations; secondly, two microseismic events located in adjacent time periods and close to each other are used to perform time difference imaging on the time difference of the first arrival of the P wave reaching the sensor, and a wave velocity time difference image is obtained; on this basis, according to the microseismic signals generated by the destruction of coal-rock formations, the source location and quantitative inversion of the focal mechanism of the destruction of coal-rock formations are performed to obtain the destruction position, fracture size and crack direction of the coal-rock formations; the dynamic and static stress fields of coal-rock formations are jointly inverted by comprehensively utilizing wave velocity field imaging, wave velocity time difference imaging, source location and quantitative parameters of focal mechanism, and the dynamic and static stress fields of coal-rock formations are obtained; finally, the abnormal areas of the dynamic and static stress field images are intelligently identified based on the deep learning algorithm, providing technical support for the monitoring of the dynamic and static stress fields of coal-rock formations.

Description

Technical Field

[0001] The invention relates to a coal rock layer stress field monitoring inversion method, and in particular to a coal rock layer dynamic and static stress field microseismic monitoring inversion and abnormal intelligent identification method, belonging to the technical field of mine rock mass mechanics and mine safety engineering. Background Art

[0002] During the construction of underground projects such as mines, tunnels, and hydropower station chambers, the coal and rock strata in the mining space are in a complex stress field environment. Coal and rock strata stress includes static stresses such as the deadweight of the overlying rock strata, tectonic stress, and mining stress, as well as dynamic stresses caused by coal and rock strata breakage, fault activation, and blasting caused by mining activities. The evolution and imbalance of the dynamic and static stress field of coal and rock strata will cause problems such as difficulty in supporting the surrounding rock of the tunnel (tunnel), deformation and damage of the surrounding rock, roof caving, and face collapse, and even induce coal and rock dynamic disasters such as rock burst, rock burst, and coal and gas outburst. Monitoring of the dynamic and static stress field of coal and rock strata is the basis and prerequisite for preventing the above-mentioned engineering problems. It is of great significance to deeply understand the evolution process of the stress field of coal and rock strata and the disaster-causing mechanism, and to improve the pertinence and advancement of disaster prevention and control.

[0003] At present, the stress of coal rock formations is usually monitored by drilling holes and burying stress sensors in coal rock formations; in addition, the working resistance of hydraulic supports and the resistance of single pillars are often used on site to monitor the stress of the working face, the surrounding rock of the roadway (tunnel) and the roof of the coal rock formation. However, the above method can only monitor the static stress (pressure) of specific measuring points in the coal rock formation, but cannot monitor the dynamic stress of the coal rock formation; and due to the limitation of measuring points, it is difficult to achieve continuous and comprehensive monitoring of the regional coal rock formation stress field.

[0004] Microseismic monitoring technology has been successfully applied in the monitoring of surrounding rock stability of underground projects such as mines and tunnels, and the monitoring and early warning of coal-rock dynamic disasters. The biggest advantage of microseismic monitoring technology is that it can locate coal-rock fractures, thereby determining information such as the time, location and energy of coal-rock fractures, monitor and evaluate the stability of coal-rock formations, and monitor and early warning of the danger of coal-rock dynamic disasters. At present, domestic and foreign scholars and field technicians mainly use this technology to study the waveform signals induced by coal-rock fractures, the location of fracture sources and the source mechanism, and passive source wave velocity tomography. However, there are few studies on the monitoring and research of coal-rock fracture-induced dynamic stress using multi-frequency acoustic wave theories and technologies such as microseismic and acoustic emission.

[0005] The principles of multi-frequency acoustic wave monitoring technologies such as microseismic, geoacoustic, and acoustic emission are the same, and the main differences are in the scale of coal-rock mass fractures, waveform frequency range, etc. The current monitoring technology is mainly used for research and application in coal-rock deformation and fracture, coal-rock dynamic and static stress monitoring and early warning, etc. For example, many research results and good application effects have been achieved in the acoustic emission signal characteristics of rock sample deformation and destruction, fracture source location and fracture mechanism analysis, and microseismic monitoring and early warning of coal-rock dynamic disasters. Based on microseismic positioning monitoring, some scholars have studied the velocity imaging of vibration waves and its application in the monitoring and evaluation of coal-rock dynamic disasters. However, they rarely involve the study of the effect of velocity changes on stress changes and deformation and destruction, and there is no research on the automatic identification and interpretation of abnormal areas in the velocity (change) image. Similarly, most of the existing studies analyze the regional dynamic and static stress based on the inversion results of a single parameter, and rarely conduct joint inversion analysis of parameter results obtained by multiple methods. In addition, although microseismic, acoustic emission and other monitoring technologies have achieved a lot of research results in coal-rock fracture monitoring, there are few reports on the research on the inversion of fracture dynamic and static stress based on microseismic monitoring.

[0006] At present, there is an urgent need for a microseismic monitoring inversion and intelligent anomaly identification method for the dynamic and static stress field of coal and rock formations, which can continuously and comprehensively monitor the dynamic and static stresses of coal and rock formation fractures and their changes, and can also intelligently identify abnormal areas to quickly provide disaster warning signals. Summary of the invention

[0007] The purpose of the present invention is to provide a method for microseismic monitoring inversion and intelligent anomaly identification of dynamic and static stress fields in coal and rock formations, which can analyze the waveform signal data obtained from microseismic monitoring to obtain source location, wave velocity field imaging, wave velocity time difference imaging and quantitative inversion parameters of the focal mechanism, and jointly invert them to obtain dynamic and static stress field images of coal and rock formations. On this basis, the stress anomaly areas in the monitoring area are intelligently identified through a deep learning algorithm, providing comprehensive monitoring and disaster warning for the coal and rock formation rupture process.

[0008] The present invention adopts the following technical solution:

[0009] The present invention provides a method for microseismic monitoring inversion and abnormal intelligent identification of dynamic and static stress fields of coal and rock strata, which specifically comprises the following steps:

[0010] a. According to the actual situation of key monitoring areas and mining space, select a reasonable number of microseismic sensors, install microseismic signal data acquisition instruments, optimize the spatial layout of microseismic sensors, determine the optimal microseismic network layout plan, and ensure comprehensive monitoring of the target area;

[0011] b. Drill long holes into the coal and rock layers in the key monitoring areas and install microseismic sensors at the bottom of the long holes. The drilling depth ranges from several meters to several hundred meters.

[0012] c. After the microseismic sensor is installed, connect the power supply, microseismic data acquisition instrument, GPS time synchronization signal and other subsystems to form the entire microseismic monitoring system. After debugging and testing, start the microseismic monitoring system to continuously collect microseismic waveform data of coal and rock strata fractures;

[0013] d. Automatically detect and pick up the collected effective waveform signals, and use the simplex and double difference joint positioning algorithm to pre-position the earthquake source caused by the coal and rock stratum fracture, and further screen to obtain high-precision positioning events;

[0014] e. Divide the detected microseismic event data set and P-wave arrival data into different time periods T1, T2, ...Ti, ...Tn according to the time window of 24h, and then further divide the time periods T1, T2, ...Ti, ...Tn according to the number of microseismic events. If the number of microseismic events in a time period TPi is ≤10, then merge the time period Ti with Ti+1 to get Ti', and divide the new time periods T1', ..., T2', Ti', ..., Tn' in the same way;

[0015] f. Determine the initial wave velocity model V by using on-site blasting test;

[0016] g. Construct the absolute arrival time, event pair arrival time difference and sensor pair arrival time difference data of each time period T1', T2', ... Ti'..., Tn' as input, perform vibration wave velocity inversion on each sub-data set with the same initial velocity model V, obtain the velocity model and crack source parameters of each time period T1', T2', ... Ti'..., Tn' through joint inversion, and finally obtain the wave velocity field imaging map of the monitoring area;

[0017] h. Use two microseismic events located in adjacent time periods T1' and T2' and close to each other to perform time difference imaging on the first arrival time difference of the P wave reaching the sensor, and obtain the wave velocity time difference imaging map of the monitoring area;

[0018] i.According to The microseismic event density is obtained, where n is the number of microseismic events, A0 is the volume of the event coverage area, A0 = (x max -x min )(y max -y min )(z max -z min ), using microseismic event density cloud map to reveal the fracture area of ​​coal and rock formations;

[0019] j. Calculate the relative acoustic emission energy E of each rupture event r, the total energy of the accumulated events and the maximum energy of the events in each time period T1, T2, ...Ti, ...Tn are counted to reflect the time distribution law of the coal-rock stratum fracture activities and infer the distribution and evolution trend of the internal stress of the coal-rock stratum;

[0020] k. Establish a crack source model involving only displacement discontinuity behavior, and quantitatively invert to obtain the tensile shear properties, rupture scale, spatial orientation and movement direction of the crack source;

[0021] l. Detect the stress drop generated when coal and rock layers break;

[0022] m. Comprehensively utilize the energy obtained, the fracture properties of the crack source, and the stress drop to analyze the fracture mechanism of the coal-rock layer and the distribution and evolution trend of the stress in the coal-rock layer, jointly characterize the stress conditions of the source point in the monitoring area, and finally obtain the dynamic stress image of the monitoring area;

[0023] n. The velocity field imaging obtained in step g reflects the stress conditions in the monitoring area to identify the high stress area; the velocity time difference imaging obtained in step h reflects the local stress change to identify the area with drastic stress change; the stress conditions at the source point are obtained according to the dynamic stress image inverted in step m to identify the stress active area; the dynamic and static stress of the coal and rock formations are jointly characterized by comprehensively utilizing the velocity field imaging, the velocity time difference imaging and the dynamic stress image, and finally the dynamic and static stress field image results of the coal and rock formations are obtained;

[0024] o. Based on the Resnet model, the dynamic and static stress field image results obtained in step n are intelligently identified and interpreted, and stress abnormal areas are automatically identified to provide a basis for monitoring, early warning and prevention of coal and rock dynamic disasters.

[0025] Preferably, in step a, the spatial layout of microseismic sensors is designed and optimized based on the geometric diffusion effect, directional control effect and influence of the directional angle on earthquake source positioning of the microseismic network.

[0026] Preferably, in step d, the threshold value method, the long and short time window method or the Akazawa information criterion method are used to automatically detect and identify the effective waveform signal, and the source positioning parameters obtained by pre-positioning the rupture source include the source time, spatial coordinates and energy. The pre-positioning results are further screened by the number of effective waveforms, event residuals, sensitivity indicators, and trigger sequence indicators in order to obtain high-precision positioning events. The accuracy of event positioning directly affects the results of processing such as energy, wave velocity imaging and wave velocity time difference imaging.

[0027] Preferably, double-difference tomography is used when performing velocity field inversion in step g.

[0028] Preferably, in step m, a wave velocity time difference imaging method is used when performing the wave velocity time difference imaging, and the specific steps include:

[0029] (1) Based on the data sets of different time periods (T1', T2', ... Ti', ... Tn') obtained in step e, the first arrival times of the P waves automatically picked up in each adjacent time period (T1' and T2') are collected;

[0030] (2) The arrival time difference data of event pairs are constructed using the first arrival time of the P wave in adjacent time periods (T1' and T2') for input; where the arrival time difference residual of event pair i, j at sensor k is The calculation formula of the time difference data is: also The calculation formula can also be In the formula is the observed event arrival time difference, is the theoretical event arrival time difference;

[0031] (3) The velocity model (V1) of TP1' obtained in step (1) remains unchanged, and the velocity change of T2' relative to TP1' and other source parameters are obtained by joint inversion; wherein the observed time difference of event pair i, j in adjacent time periods is obtained by using the monitored absolute arrival time data Based on the current velocity model (the initial model of T1' is V1, and the initial model of T2' is also V1, because the velocity model of TP2' will not change much relative to the velocity model of T1', which has a good constraint on the inversion velocity change) and source parameters, forward simulation is performed using ray tracing methods to calculate the theoretical arrival time difference of the event pair. The observed arrival time difference is then compared with the theoretical arrival time difference to obtain the arrival time difference residual of the event pair. The velocity model and source parameters are updated through the inversion algorithm. The goal is to find a set of parameters (including velocity model parameters and source parameters) that minimize the travel time difference residual between the theory and observation. In each iteration, the inversion algorithm adjusts the velocity model and source parameters according to the travel time difference data, and then re-performs the forward simulation and calculates the new travel time difference. Through continuous iteration and parameter adjustment, the travel time difference is gradually reduced, thereby updating and optimizing the T2' velocity model and source parameters, and finally outputting the wave velocity change of T2' relative to T1'.

[0032] (4) Obtain the wave velocity change results of each adjacent time period, and finally output the wave velocity change image of the monitoring area.

[0033] Preferably, when calculating the relative energy of the rupture event in step f, the signal energy is first calculated, and the root mean square value of the acoustic emission waveform voltage signal in the ISRM standard is used to characterize the acoustic emission energy, that is, Where T1 is the signal triggering time, T2 is the ending time, T is the duration of the signal fluctuation, and V(t) represents the signal amplitude. Then, considering the geometric attenuation effect of the acoustic emission waveform propagating in the medium, that is, the farther away from the earthquake source, the smaller the waveform amplitude recorded by the sensor, so the relative acoustic emission energy of the rupture event considering energy attenuation is Where RMS i represents the acoustic emission energy collected by the i-th sensor, N is the number of sensors involved in positioning, r i Represents the spatial distance between the i-th sensor and the event location.

[0034] Preferably, when performing focal mechanism inversion to obtain crack focal parameters in step g, the first-arrival P wave data is used for inversion, a crack focal model involving only displacement discontinuity behavior is established, and the displacement discontinuity tensor constraint condition in the crack focal model is converted into a moment tensor invariant form; the Lagrange multiplier combined with the Levenberg-Marquardt optimization algorithm is used to iteratively solve the moment tensor component under the constraint condition that the intermediate eigenvalue of the displacement discontinuity tensor is zero, so that the theoretical normal displacement and the measured normal displacement of the P wave at different sensor positions are consistent. The least squares error between displacements is the smallest; the eigenvalues ​​of the displacement discontinuity tensor and its eigenvectors are calculated from the solved moment tensor components; the crack mode angle α and the newly generated crack volume |b|ΔA are calculated based on the eigenvalues ​​of the displacement discontinuity tensor; based on the crack mode angle α and the eigenvectors of the displacement discontinuity tensor, the crack normal n and the crack displacement direction b, finally the crack source parameters, namely the tension-shear properties (determined by the crack mode angle α), the breaking scale (determined by the crack volume |b|ΔA), the spatial orientation (crack normal n), and the movement direction (crack displacement direction b), are obtained.

[0035] Preferably, in step o, the abnormal area of ​​the dynamic and static stress images is intelligently identified, and the specific steps include:

[0036] (1) Data preparation stage: Collect the dynamic and static stress field cloud map data set for joint inversion, including the training set D train and the test set D test , each sample in the data set contains a dynamic and static stress field image X and a corresponding label Y, where the label represents the abnormal area in the dynamic and static stress field; the dynamic and static stress field cloud map data is preprocessed, including noise removal, missing value filling, normalization, etc. The normalization formula is: Among them, μ is the mean of the dynamic and static stress field images, and σ is the standard deviation of the dynamic and static stress field images;

[0037] (2) Model training stage: The training set D train Input into the ResNet model, and after operations such as convolution and residual connection, the output of the model is obtained Use the cross entropy loss function to train the model and optimize the model parameters. The formula of the loss function is: Among them, N is the number of training samples, C is the number of categories, and Y ij is the true label, is the model output;

[0038] (3) Model application stage: Test set D test Input into the trained ResNet model to obtain the prediction results of the dynamic and static stress field images The prediction results are converted into binary abnormal area images using thresholds or other methods, which are used to automatically identify stress abnormal areas in dynamic and static stress fields;

[0039] (4) Result analysis phase: Use evaluation indicators (such as accuracy, recall, F1 score, etc.) to evaluate the performance of the model. The calculation formula of the evaluation indicator is: Among them, TP (True Positive) represents true positive, TN (True Negative) represents true negative, FP (False Positive) represents false positive, and FN (False Negative) represents false negative;

[0040] The dynamic and static stress field cloud map is analyzed using classical methods, and compared with the results obtained by the ResNet model. Through comparative analysis, the performance and accuracy of the model can be evaluated, possible problems in the model can be discovered and improved, and the parameters of the ResNet model can be optimized, such as adjusting the network structure and the learning rate, to further improve the performance and generalization ability of the model. The performance of the model is evaluated through cross-validation and other methods, and the optimal parameter settings are selected. A dynamic and static stress field analysis method based on the ResNet model is studied to identify stress anomaly areas, and it is combined with traditional methods to improve the accuracy and reliability of identifying anomaly areas. Through continuous optimization and improvement, the model can better adapt to the needs of dynamic and static stress field data analysis under different geological conditions.

[0041] The beneficial effects of the present invention are:

[0042] 1. The present invention jointly characterizes the earthquake source location results, energy, quantitative inversion parameters of the earthquake source mechanism, and the wave velocity field imaging and wave velocity time difference imaging results to obtain the dynamic and static stress field images of the coal rock layer. First, the stress situation of the entire area is characterized according to the wave velocity field inversion image, and then the stress situation of the local area is characterized according to the wave velocity field time difference inversion image. Finally, the stress situation of the point is characterized according to the dynamic stress image, and finally the dynamic and static stress images of the coal rock layer are obtained, which can provide more comprehensive and accurate information related to earthquake events and coal rock layer stress fields, and is of great significance for predicting and preventing related geological disasters and ensuring safe production in coal mines.

[0043] 2. The present invention further performs a joint inversion of the dynamic and static stress fields on the coal-rock strata fractures based on the use of microseismic monitoring technology to monitor the coal-rock strata fractures. By comprehensively utilizing the real-time monitoring of the coal-rock strata fractures and the dynamic and static stress field inversion results, abnormal stress concentration or release phenomena in the coal-rock strata can be discovered in time, thereby taking preventive measures in advance to avoid accidents.

[0044] 3. The present invention adopts the wave velocity time difference imaging method, that is, directly using the travel time difference data of the earthquake pair in two adjacent time periods to the same station to obtain the wave velocity change by studying the inversion of the target function, while the traditional method only subtracts the velocity model obtained in two different time periods to obtain the velocity change. In addition, this method requires that the two events constituting the event pair are respectively in two adjacent time periods, so this method will be able to greatly reduce the impact of different data distributions in different periods on the velocity change results, which is helpful to improve the application potential of microseismic monitoring in coal mine rock burst disaster warning. Wave velocity time difference imaging can monitor the change of wave velocity in the process of coal and rock strata rupture in real time, so as to timely discover potential safety hazards. In addition, the propagation speed of seismic waves is affected by the underground stress field. The change of the underground stress field can be inferred by wave velocity time difference imaging, which is of great significance for studying seismic activity and rock mechanical properties.

[0045] 4. The present invention studies and forms a dynamic and static stress field analysis method based on the ResNet model to identify stress anomaly areas, and combines it with traditional methods to improve the accuracy and reliability of identifying anomaly areas. Intelligent recognition technology can greatly reduce the labor and time required for manual recognition, and realize real-time monitoring and analysis of stress anomaly areas in coal strata, so as to timely alarm and take measures. In addition, since the deep learning algorithm has powerful feature learning and pattern recognition capabilities, it can automatically extract key features related to rupture events from complex seismic data, which greatly improves the accuracy and reliability of stress anomaly area identification and reduces interference and misjudgment of human factors. Through continuous optimization and improvement, the model can better adapt to the needs of wave velocity field data analysis under different geological conditions. Intelligent identification of stress anomaly area technology can provide strong technical support for safe production and geological disaster prevention and control in coal mines. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0047] Figure 1 It is a flow chart of the microseismic monitoring inversion and abnormal intelligent identification method of the dynamic and static stress field of coal and rock strata of the present invention;

[0048] Figure 2 It is a joint representation diagram of wave velocity field and dynamic stress in the coal and rock layer monitoring area of ​​the present invention;

[0049] Figure 3 It is a joint characterization diagram of wave velocity time difference and dynamic stress in the coal and rock layer monitoring area of ​​the present invention;

[0050] Figure 4 This is a result diagram of the present invention for intelligently identifying abnormal areas in the coal and rock layer monitoring area; DETAILED DESCRIPTION

[0051] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] like Figure 1 As shown, the present invention provides a method for microseismic monitoring inversion and abnormal intelligent identification of dynamic and static stress fields of coal and rock formations, which specifically includes the following steps:

[0053] a. According to the actual situation of key monitoring areas and mining space, select a reasonable number of microseismic sensors to design the microseismic sensor layout plan. Under the conditions of considering the actual situation and economic factors, increase the number of sensors in the microseismic network as much as possible, and try to make the microseismic network form an irregular spatial layout, which can increase the difference between the arrival times of sensors, thereby reducing the relative arrival time error and improving the accuracy of earthquake source positioning; at the same time, in actual applications, the sensor installation site should be as far away as possible from large electrical appliances, high-voltage cables and mechanical equipment to avoid interference from electromechanical equipment; and optimize the spatial layout plan of microseismic sensors based on the geometric diffusion effect, directional control effect and the influence of azimuth angle on earthquake source positioning of the microseismic network, and determine the optimal microseismic network layout plan,

[0054] b. Use directional drilling rigs to construct long boreholes in the coal and rock layers in the key monitoring areas in mining spaces such as tunnels and stopes, and use push rods to install microseismic sensors at the bottom of the long boreholes. The drilling depth is several meters to several hundred meters, and the microseismic sensors and the bottom of the hole must be coupled through coupling agents;

[0055] c. Install microseismic sensors according to the optimal microseismic network layout plan, and then connect the power supply, microseismic data acquisition instrument, GPS time synchronization signal and other subsystems to form the entire microseismic monitoring system. After the entire microseismic monitoring system is installed, debug and test the system to ensure that all microseismic sensors can monitor effective waveform signals, and then start the microseismic monitoring system to continuously collect microseismic waveform data of coal and rock strata fractures;

[0056] d. Automatically detect and identify the collected effective waveform signals by using the threshold value method, long-short time window method or Akazawa information criterion method, and pre-locate the earthquake source generated by coal-rock stratum fracture using the simplex and double-difference joint positioning algorithm. The earthquake source positioning parameters obtained by pre-locating the fracture source include the earthquake source time and spatial coordinates; further screen the pre-location results by the number of effective waveforms, event residuals, sensitivity indexes, and trigger sequence indexes in order to obtain high-precision positioning events;

[0057] e. Divide the detected microseismic event data set and P-wave arrival data into different time periods (TP1, TP2, ...TPi, ...TPn) according to the time window of 24h, and then further divide the time periods (TP1, TP2, ...TPi, ...TPn) according to the number of microseismic events. If the number of microseismic events in a time period TPi is ≤10, then merge the time period TPi with TPi+1 to obtain TPi', and divide the new time periods (TP1', ..., TP2', TPi', ..., TPn') in this way;

[0058] f. Determine the initial wave velocity model V by using on-site blasting test;

[0059] g. Collect the P-wave arrival time in each time period (TP1', TP2', ...TPi'..., TPn'), construct the absolute arrival time, event pair arrival time difference and sensor pair arrival time difference data of each time period (TP1', TP2', ...TPi'..., TPn') as input, the two microseismic events constituting the event pair should be close to each other and received by a sufficient number of stations, the event pair data can appropriately eliminate the ray path outside the source area, to a certain extent, the fine velocity structure of the source area can be obtained, and the waveform cross-correlation technology can be used to more accurately obtain the observation travel time of the sensor pair, thereby effectively eliminating the problem of inaccurate arrival time picking or no clear arrival time of some low signal-to-noise ratio signals; among them, the travel time residual calculation formula of a single event is In the formula represents the travel time residual from earthquake event i to station k, Δτ i represents the disturbance at the time of occurrence of earthquake event i, represents the disturbance of earthquake position, δu represents the disturbance of slowness, S k represents the disturbance of the sensor; the calculation formula for the time difference data of the story event pair is In the formula represents the residual of the arrival time difference of event pair i, j at sensor k; The calculation formula can also be In the formula is the observed event arrival time difference, is the theoretical arrival time difference of the event pair; the arrival time difference calculation formula of the sensor pair is In the formula represents the residual of the arrival time difference of event i at sensor pair k, l; where, It can also be expressed as In the formula is the observed sensor time difference data, is the theoretical sensor pair arrival time difference data; the velocity model and other source parameters (seismic time and position) in each time period are jointly inverted with the same initial velocity model V; wherein, the initial wave velocity model V is adopted to calculate the theoretical travel time field by the ray tracing method, and the observed travel time field of the sensor pair and the observed travel time field of the time pair can be obtained by using the cross-correlation technology and the absolute arrival time of event picking, and then the residual of the theoretical travel time field and the observed travel time field is used to update the model and other source parameters (seismic time and position) based on the linear optimization or nonlinear optimization method until the residual converges; according to the obtained wave velocity model (V1, V2, ...Vi..., Vn) of each time period, the wave velocity field image of the monitoring area is finally obtained as shown in Figure 2 As shown, the stress in the region with larger wave velocity is also greater.

[0060] h. Collect the P wave first arrival time automatically picked up in each adjacent time period (TP1' and TP2'), and use the P wave first arrival time of the adjacent time period (TP1' and TP2') to construct the arrival time difference data of the event pair for input, where the calculation formula for the arrival time difference residual of event pair i, j to sensor k is: also The calculation formula can also be In the formula is the observed event arrival time difference, is the theoretical arrival time difference of the event pair; the velocity model (V1) of TP1' obtained from the velocity field imaging remains unchanged, and the velocity change of TP2' relative to TP1' and other source parameters are obtained through joint inversion; among them, the observed arrival time difference of event pair i, j in adjacent time periods is obtained by using the monitored absolute arrival time data Based on the current velocity model (the initial model of TP1' is V1, and the initial model of TP2' is also V1, because the velocity model of TP2' will not change much relative to the velocity model of TP1', which has a good constraint on the inversion velocity change) and source parameters, forward simulation is performed using ray tracing methods to calculate the theoretical arrival time difference of the event pair. The observed arrival time difference is then compared with the theoretical arrival time difference to obtain the arrival time difference residual of the event pair. The velocity model and source parameters are updated through the inversion algorithm. The goal is to find a set of parameters (including velocity model parameters and source parameters) that minimize the travel time difference residual between the theory and the observation. In each iteration, the inversion algorithm adjusts the velocity model and source parameters according to the travel time difference data, and then re-performs the forward simulation and calculates the new travel time difference. Through continuous iteration and parameter adjustment, the travel time difference is gradually reduced, thereby updating and optimizing the TP2' velocity model and source parameters, and finally outputting the wave velocity change of TP2' relative to TP1', and obtaining the wave velocity time difference imaging map of the monitoring area as shown in Figure 3 As shown, local wave velocity changes reflect stress changes.

[0061] i.According to The microseismic event density is obtained, where n is the number of microseismic events, A0 is the volume of the event coverage area, A0 = (x max -x min )(y max -y min )(z max -z min ), using the density cloud map of microseismic events to vividly reveal the fracture area of ​​coal and rock layers;

[0062] j. Calculate the relative acoustic emission energy E of each rupture event r , first calculate the signal energy, and use the ISRM standard to characterize the acoustic emission energy by the root mean square value of the acoustic emission waveform voltage signal, that is, Where T1 is the signal triggering time, T2 is the ending time, T is the duration of the signal fluctuation, and V(t) represents the signal amplitude. Then, considering the geometric attenuation effect of the acoustic emission waveform propagating in the medium, that is, the farther away from the earthquake source, the smaller the waveform amplitude recorded by the sensor, so the relative acoustic emission energy of the rupture event considering energy attenuation is Where RMS i represents the acoustic emission energy collected by the i-th sensor, N is the number of sensors involved in positioning, r iIndicates the spatial distance between the ith sensor and the event location; statistics the total event energy and the maximum event energy in each time period (TP1, TP2, ...TPi, ...TPn), reflects the temporal distribution law of coal-rock stratum fracture activities, and infers the distribution and evolution trend of internal stress in coal-rock strata;

[0063] k. Through quantitative inversion of the focal mechanism, the direct parameters of the crack source, such as the tensile shear properties, rupture scale, spatial orientation and movement direction, are obtained. The first P-wave data are used for inversion to establish a crack source model involving only the displacement discontinuity (DD) behavior (the constraint condition that the intermediate eigenvalue of the DD tensor is zero must be satisfied). The moment tensor is solved using the least squares method based on the normal displacement of the P wave at different sensor locations and the measured normal displacement data. In the formula is the normal displacement measured by the i-th sensor, which is represented by the scalar sensitive parameter S A and measure the P wave onset amplitude Multiplying them together, is the moment tensor M at the i-th sensor position x jk The calculated point source P wave theoretical normal displacement is the radial displacement u i =G ij,k (x,t)*M jk (x 0 ,t) and reflection coefficient The product of , β is the angle between the incident P wave and the sensor normal, is the ratio of the P-wave velocity to the S-wave velocity of the sample; combined with the constraint that the intermediate eigenvalue of the DD tensor is zero, the Lagrange multiplier combined with the Levenberg-Marquardt method is used to solve the moment tensor M under the constraint to obtain the 6 components m1~m6 of the moment tensor; the eigenvalue and eigenvector of the DD tensor ψ are calculated according to the solved moment tensor components, where the eigenvector of the DD tensor coincides with the moment tensor M, and the calculation formula of the eigenvalue is: Where λ and μ are Lamé constants, M i is the eigenvalue of the matrix quantity, and then according to The crack mode angle α is calculated by |bΔA=ψ1-ψ3 to determine the tensile shear property of the crack and the crack volume |bΔA to reflect the fracture scale. Finally, the crack spatial orientation, crack normal n and crack movement direction b are obtained according to the relationship between the crack mode angle α and the characteristic vector of the DD tensor ψ. Among them, the tensile shear property reveals the nature of the fracture event, so as to understand the fracture mechanism and stress state of the coal rock formation; the fracture scale provides the size information of the fracture area to reflect the range of stress accumulation and release; the parameters of spatial orientation and movement direction reveal the directionality and expansion trend of the fracture event to understand the transmission path and direction of stress in the coal rock formation;

[0064] l. Detect the stress drop generated when coal and rock layers break;

[0065] m. Calculate the energy of each earthquake source point. In addition, calculate the crack mode angle α according to the eigenvalue of the displacement discontinuity tensor to determine the rupture type of large cracks, i.e. tension or shear. Calculate the rupture scale according to the crack volume |b|ΔA. According to the crack mode angle α and the eigenvector of the displacement discontinuity tensor, the crack normal n, i.e. the spatial orientation of the crack, and the crack displacement direction b, i.e. the movement direction, and according to the size of these calculated parameters, the crack earthquake source point is indicated in the figure, and finally the dynamic stress image of the monitoring area is obtained. Figure 2 and Figure 3 ;

[0066] n. Reflect the stress conditions in the monitoring area based on the velocity field imaging to identify the high stress area; reflect the local stress changes based on the velocity time difference image to identify the area with drastic stress changes; obtain the stress conditions at the source point based on the inverted dynamic stress image to identify the stress activity area; comprehensively use the velocity field imaging, velocity time difference imaging and dynamic stress image to jointly characterize the dynamic and static stress of the coal and rock formations, and finally obtain the dynamic and static stress field image results of the coal and rock formations, such as Figure 2 and Figure 3 shown.

[0067] o. Based on the Resnet model, the obtained dynamic and static stress field images are intelligently identified with stress anomaly areas, and the dynamic and static stress field cloud image data sets of joint inversion are collected, including the training set D train and the test set D test Each sample in the dataset

[0068] Contains dynamic and static stress field images X and corresponding labels Y; preprocessing of dynamic and static stress field cloud map data, including noise removal, missing value filling, normalization, etc. The normalization formula is: The training set D train Input into the ResNet model, and after operations such as convolution and residual connection, the output of the model is obtained Use the cross entropy loss function to train the model and optimize the model parameters. The formula of the loss function is: The test set D test Input into the trained ResNet model to obtain the prediction results of the dynamic and static stress field images The prediction results are converted into binary abnormal area images using thresholds or other methods for automatic identification of stress abnormal areas in dynamic and static stress fields. Evaluation indicators (such as accuracy, recall rate, F1 score, etc.) are used to evaluate the performance of the model. The calculation formula of the evaluation indicators is:

[0069]

[0070] Among them, TP (True Positive) represents true positive, TN (True Negative) represents true negative, FP (False Positive) represents false positive, and FN (False Negative) represents false negative;

[0071] The dynamic and static stress field cloud map is analyzed using classical methods, and compared with the results obtained by the ResNet model. The performance and accuracy of the model are evaluated through comparative analysis, possible problems in the model are discovered and improved, and the parameters of the ResNet model are optimized. The performance and generalization ability of the model are further improved by adjusting the network structure and the learning rate. The performance of the model is evaluated through cross-validation and other methods, and the optimal parameter settings are selected. The dynamic and static stress field analysis method based on the ResNet model is studied to identify stress anomaly areas, and it is combined with traditional methods to improve the accuracy and reliability of identifying anomaly areas. Through continuous optimization and improvement, the model can better adapt to the dynamic and static stress field data analysis needs under different geological conditions. The final anomaly area recognition image is shown in the figure. Figure 4 As shown, based on the potential probability of the abnormal area, the area where stress concentration may exist can be identified, and the possibility of danger occurring in this area is relatively high.

Claims

1. A method for inversion of microseismic monitoring of dynamic and static stress fields of coal and rock formations and intelligent identification of abnormalities, characterized by: The specific steps are as follows: a. According to the actual situation of key monitoring areas and mining space, select a reasonable number of microseismic sensors, install microseismic signal data acquisition instruments, optimize the spatial layout of microseismic sensors, determine the optimal microseismic network layout plan, and ensure comprehensive monitoring of the target area; b. Drill long holes into the coal and rock layers in the key monitoring areas and install microseismic sensors at the bottom of the long holes. The drilling depth ranges from several meters to several hundred meters. c. After the microseismic sensor is installed, connect the power supply, microseismic data acquisition instrument, GPS time synchronization signal and other subsystems to form the entire microseismic monitoring system. After debugging and testing, start the microseismic monitoring system to continuously collect microseismic waveform data of coal and rock strata fractures; d. Automatically detect and pick up the collected effective waveform signals, and use the simplex and double difference joint positioning algorithm to pre-position the earthquake source caused by coal and rock stratum fractures, and further screen to obtain high-precision positioning events; e. Divide the detected microseismic event data set and P-wave arrival time data into different time periods T1, T2, ...Ti, ...Tn according to the time window of 24h, and then further divide the time periods T1, T2, ...Ti, ...Tn according to the number of microseismic events. If the number of microseismic events in a time period TPi is ≤10, then merge the time period Ti with Ti+1 to get Ti', and divide the new time periods T1', ..., T2', Ti', ..., Tn' in the same way; f. Determine the initial velocity model V by using on-site blasting test; g. Construct the absolute arrival time, event pair arrival time difference and sensor pair arrival time difference data of each time period T1', T2', ... Ti'..., Tn' as input, perform vibration wave velocity inversion on each sub-data set with the same initial velocity model V, obtain the velocity model and crack source parameters of each time period T1', T2', ... Ti'..., Tn' through joint inversion, and finally obtain the wave velocity field imaging map of the monitoring area; h. Use two microseismic events located in adjacent time periods T1' and T2' and close to each other to perform time difference imaging on the first arrival time difference of the P wave reaching the sensor, and obtain the wave velocity time difference imaging map of the monitoring area; i. According to The microseismic event density is obtained as follows: is the number of microseismic events, is the event coverage area volume, , using microseismic event density cloud map to reveal the fracture area of ​​coal and rock formations; j. Calculate the relative acoustic emission energy of each rupture event , the total energy of the accumulated events and the maximum energy of the events in each time period T1, T2, ...Ti, ...Tn are counted to reflect the time distribution law of the coal-rock stratum fracture activities and infer the distribution and evolution trend of the internal stress of the coal-rock stratum; k. Establish a crack source model involving only displacement discontinuity behavior, and quantitatively invert to obtain the tensile shear properties, rupture scale, spatial orientation and movement direction of the crack source; 1. Detect the stress drop generated when coal and rock layers break; m. Comprehensively utilize the energy obtained, the fracture properties of the crack source, and the stress drop to analyze the fracture mechanism of the coal-rock layer and the distribution and evolution trend of the stress in the coal-rock layer, jointly characterize the stress situation of the source point in the monitoring area, and finally obtain the dynamic stress image of the monitoring area; n. The wave velocity field imaging obtained in step g reflects the stress conditions in the monitoring area to identify the high stress area; the wave velocity time difference imaging obtained in step h reflects the local stress change conditions to identify the area with drastic stress changes; the dynamic stress image inverted in step m obtains the stress conditions at the earthquake source to identify the stress active area; The dynamic and static stress of coal and rock formations are jointly characterized by comprehensive use of wave velocity field imaging, wave velocity time difference imaging and dynamic stress images, and finally the dynamic and static stress field image results of coal and rock formations are obtained; o. Based on the Resnet model, the dynamic and static stress field image results obtained in step n are intelligently identified and interpreted, and the stress abnormality area is automatically identified, providing a basis for monitoring, early warning and prevention of coal and rock dynamic disasters.

2. The method for microseismic monitoring, inversion and intelligent abnormality identification of dynamic and static stress fields of coal and rock formations according to claim 1 is characterized by: In step a, the spatial layout of microseismic sensors is designed and optimized based on the geometric diffusion effect, directional control effect and influence of the azimuth angle on earthquake source location of the microseismic network.

3. The method for microseismic monitoring, inversion and intelligent abnormality identification of dynamic and static stress fields of coal and rock formations according to claim 1 is characterized by: In step d, the threshold value method, the long-short time window method or the Akasawa information criterion method are used to automatically detect and identify the effective waveform signal, and the source positioning parameters obtained by pre-positioning the rupture source include the source time, spatial coordinates and energy.

4. The method for microseismic monitoring inversion and abnormal intelligent identification of dynamic and static stress fields of coal and rock formations according to claim 1 is characterized by: Double-difference tomography is used when performing velocity field inversion in step g.

5. The method for microseismic monitoring inversion and abnormal intelligent identification of dynamic and static stress fields of coal and rock strata according to claim 1 is characterized by: The wave velocity time difference imaging is performed in step m using a wave velocity time difference imaging method, and the specific steps include: (1) according to the data set of different time periods T1', T2', ... Ti', ... Tn' obtained in step e of claim 1, collecting the first arrival times of the P waves automatically picked up in each adjacent time period T1' and T2'; (2) The arrival time difference data of the event pair is constructed by using the first arrival time of the P wave in the adjacent time periods T1' and T2' for input; where the event pair Arrival sensor The arrival time difference residual The calculation formula is ;also The calculation formula can also be , where is the observed event arrival time difference, is the theoretical event arrival time difference; (3) The velocity model V1 obtained by velocity field imaging at T1' remains unchanged, and the velocity change of T2' relative to T1' and other source parameters are obtained by joint inversion; the observed arrival time difference of event pair i, j in adjacent time periods is obtained by using the absolute arrival time data monitored Based on the current velocity model and source parameters, the forward modeling method is used to calculate the theoretical arrival time difference of the event pair. ; Then the observed arrival time difference is compared with the theoretical arrival time difference to obtain the arrival time difference residual of the event pair, and the velocity model and source parameters are updated through the inversion algorithm. The goal is to find a set of velocity model parameters and source parameters that minimize the arrival time difference residual between the theory and the observation; in each iteration, the inversion algorithm adjusts the velocity model and source parameters according to the arrival time difference data, and then re-performs the forward simulation and calculates the new arrival time difference; through continuous iteration and parameter adjustment, the arrival time difference is gradually reduced, thereby updating and optimizing the T2' velocity model and source parameters, and finally outputting the wave velocity change of T2' relative to T1'; (4) Obtain the wave velocity variation results of each adjacent time period, and finally output the wave velocity time difference image of the monitoring area.

6. The method for microseismic monitoring, inversion and intelligent abnormality identification of dynamic and static stress fields of coal and rock formations according to claim 1 is characterized by: When calculating the relative energy of the rupture event in step f, the signal energy is calculated first, and the root mean square value of the acoustic emission waveform voltage signal in the ISRM standard is used to characterize the acoustic emission energy, that is, , where is the signal triggering moment, To end the moment, is the duration of signal fluctuation, represents the signal amplitude; then considering the geometric attenuation effect of the acoustic emission waveform propagating in the medium, that is, the farther away from the earthquake source, the smaller the waveform amplitude recorded by the sensor, so the relative acoustic emission energy of the rupture event considering energy attenuation is , where represents the acoustic emission energy collected by the i-th sensor, is the number of sensors involved in positioning, Represents the spatial distance between the i-th sensor and the event location.

7. The method for microseismic monitoring, inversion and intelligent abnormality identification of dynamic and static stress fields of coal and rock formations according to claim 1 is characterized by: When performing focal mechanism inversion to obtain crack source parameters in step g, the first P wave data is used for inversion, and a crack source model involving only displacement discontinuity behavior is established, and the displacement discontinuity tensor constraint conditions in the crack source model are converted into moment tensor invariant forms; the Lagrange multiplier is combined with the Levenberg-Marquardt optimization algorithm to iteratively solve the moment tensor components under the constraint condition that the intermediate eigenvalue of the displacement discontinuity tensor is zero, so that the least squares error between the theoretical normal displacement and the measured normal displacement of the P wave at different sensor positions is minimized; the displacement discontinuity tensor eigenvalue and its eigenvector are calculated from the solved moment tensor components; the crack mode angle α and the newly generated crack volume |b|ΔA are calculated based on the displacement discontinuity tensor eigenvalue; based on the crack mode angle α and the eigenvector crack normal n and the crack dislocation direction b of the displacement discontinuity tensor, the crack source parameters, namely the tension-shear properties, the breaking scale, the spatial orientation, and the movement direction, are finally obtained.

8. The method for microseismic monitoring, inversion and intelligent abnormality identification of dynamic and static stress fields of coal and rock strata according to claim 1 is characterized by: In step o, the abnormal area of ​​the dynamic and static stress images is intelligently identified, and the specific steps include: (1) Data preparation stage: Collect the dynamic and static stress field cloud map data sets for joint inversion, including the training set and test set , each sample in the dataset contains dynamic and static stress field images and the corresponding labels , where the label represents the abnormal area in the dynamic and static stress field; the dynamic and static stress field cloud map data is preprocessed, including removing noise, filling missing values, normalization, etc. The normalization formula is: ,in, is the mean value of the dynamic and static stress field images, is the standard deviation of the dynamic and static stress field images; (2) Model training phase: The training set Input into the ResNet model, and after operations such as convolution and residual connection, the output of the model is obtained ; Use the cross entropy loss function to train the model and optimize the model parameters. The formula of the loss function is: in, is the number of training samples, is the number of categories, is the true label, is the model output; (3) Model application stage: Test set Input into the trained ResNet model to obtain the prediction results of the dynamic and static stress field images The prediction results are converted into binary abnormal area images using thresholds or other methods, which are used to automatically identify stress abnormal areas in dynamic and static stress fields; (4) Result analysis stage: The performance of the model is evaluated using evaluation indicators, including accuracy, recall rate, and F1 score. The calculation formula of the evaluation indicators is: , , , , where TP represents true positive, TN represents true negative, FP represents false positive, and FN represents false negative; The dynamic and static stress field cloud map is analyzed using classical methods, and compared with the results obtained by the ResNet model. Through comparative analysis, the performance and accuracy of the model can be evaluated, possible problems in the model can be discovered and improved, and the parameters of the ResNet model can be optimized, such as adjusting the network structure and the learning rate, to further improve the performance and generalization ability of the model. The performance of the model is evaluated through cross-validation and other methods, and the optimal parameter settings are selected. A dynamic and static stress field analysis method based on the ResNet model is studied to identify stress anomaly areas, and it is combined with traditional methods to improve the accuracy and reliability of identifying anomaly areas. Through continuous optimization and improvement, the model can better adapt to the needs of dynamic and static stress field data analysis under different geological conditions.

Images