Mine earthquake source inversion method based on monitoring stress wave signals

Abstract

The present invention discloses a mine earthquake source inversion method based on monitoring stress wave signals, which relates to the technical field of mine earthquake source inversion, including arranging a monitoring network, recording stress wave signals at various monitoring points, processing stress wave signals, identifying signal characteristics of different waveforms, analyzing the stratigraphic structure, lithology distribution and fault structure of the target area, establishing a multidimensional velocity model including different stratigraphic velocities and a stratigraphic model of discontinuity dip angles, and then establishing an inversion equation including stratigraphic velocity and discontinuity dip angle parameters, and using an optimization algorithm to solve the inversion equation, analyzing the best combination of stratigraphic velocity and discontinuity dip angle, evaluating the accuracy of the inversion result, iterative inversion, correcting the inversion accuracy and results of the stratigraphic model, and determining the location range of the earthquake source. The mine earthquake source inversion method of the present invention integrates stress wave signal characteristics and stratigraphic models, and uses an optimization algorithm to solve, correct and optimize the inversion equation, and the earthquake source location range is determined with high accuracy.

Description

Technical Field

[0001] The invention relates to the technical field of mine earthquake source inversion, and in particular to a mine earthquake source inversion method based on monitoring stress wave signals. Background Art

[0002] As coal mining enters deeper areas, high ground stress and strong mining disturbances often cause induced mine earthquakes, which have a significant impact on the stability of surrounding tunnels and mining areas, such as roof falls, gas outbursts, and impact instability. Accurately monitoring stress waves to invert the location of the mine earthquake source can effectively analyze the impact range of the mine earthquake, and thus play a positive role in the prevention and control of dynamic disasters in coal mining.

[0003] For example, a virtual seismic source reflection wave pre-stack inversion method with Chinese patent publication number CN113267811A first establishes a virtual reflection wave travel time equation, and then, based on the virtual reflection wave travel time equation, takes the formation average velocity and the discontinuity surface dip angle as variables to establish an inversion equation; by scanning the inversion equation, the variation law of the energy of the superimposed virtual seismic source reflection wave with the average velocity v and the discontinuity surface dip angle φ is obtained.

[0004] The existing technology obtains the variation law of the energy of the superimposed virtual seismic source reflection wave with the average velocity v and the discontinuity surface dip angle φ through scanning the inversion equation, and realizes the pre-stack inversion of the virtual seismic source reflection wave. However, the complex changes of the formation velocity and the discontinuity surface dip angle will make the inversion problem more complicated. When the underground structure is complex, there will be a large deviation between the inversion result and the actual situation. Therefore, a mine earthquake source inversion method based on monitoring stress wave signals is proposed. Summary of the invention

[0005] The purpose of the present invention is to provide a mine earthquake source inversion method based on monitoring stress wave signals to solve the problems raised in the above-mentioned background technology.

[0006] In order to solve the above technical problems, the technical solution adopted by the present invention is: a mine earthquake source inversion method based on monitoring stress wave signals, comprising the following steps: step 1, arranging a monitoring network, synchronously recording the stress wave signals of each monitoring point, and preprocessing the captured stress wave signals; step 2, identifying the signal characteristics of different waveforms from the preprocessed stress wave signals, and calculating the propagation speed, attenuation law, and phase difference parameters of the stress wave according to the waveform characteristics; step 3, combining geological exploration data, analyzing the stratigraphic structure, lithology distribution and fault structure of the target area, and establishing a system including the stratigraphic structure, lithology distribution and fault structure of the target area according to the geological survey results. A multidimensional velocity model of different formation velocities and a formation model of discontinuity surface dip angles; Step 4, based on the monitored stress wave signal characteristics and in combination with the formation model, an inversion equation including formation velocity and discontinuity surface dip angle parameters is established, and the inversion equation is solved by an optimization algorithm to analyze the best combination of formation velocity and discontinuity surface dip angle; Step 5, the inverted formation velocity and discontinuity surface dip angle are compared and verified with geological exploration data to evaluate the accuracy of the inversion results; Step 6, the formation model is corrected and optimized according to the verification results, the inversion is iteratively performed, the inversion accuracy and results are corrected, and the location range of the earthquake source is determined.

[0007] The further improvement of the technical solution of the present invention is that in the step 1, the process of acquiring the stress wave signal is as follows: step 101, according to the operation area of ​​the mining operation, the number and location of the monitoring points are determined, and stress wave monitoring equipment, including velocity sensors, acceleration sensors and data acquisition instruments, are installed at each monitoring point to ensure accurate capture of the stress wave signal, wherein the operation area of ​​the mining operation needs to analyze the specific location of the mining operation, the scope of influence and the sensitivity of the surrounding environment, so as to determine the number and location of the monitoring points, and the monitoring points fully cover the potentially affected areas. First, the specific location of the mining operation, the expected operation scale and other key information are clarified, Based on the scale and nature of the mining operation, combined with geological conditions, topography, building structure and other factors, evaluate the propagation range and impact of stress waves that may be generated by the operation, ensure that the monitoring points can fully cover the potentially affected areas, including sensitive buildings, public facilities, residential areas, etc. At the same time, consider setting up several monitoring points near the operation source to obtain more direct vibration data; Step 102, after the stress wave monitoring equipment is deployed, check it, configure the equipment parameters, use the GPS time synchronization system to synchronize the time of all monitoring points, ensure that the signal timestamps of different monitoring points are accurately corresponding to the subsequent data analysis, and check to ensure that the monitoring equipment is fixed and stable. Solid, avoid displacement or damage during operation, at the same time, check whether the power supply, communication lines, etc. of the equipment are connected normally to ensure that the data can be transmitted to the data center in real time; step 103, when the mining operation is in progress, use the stress wave monitoring equipment to record and store the stress wave signal, the stress wave signal includes vibration velocity and acceleration parameters, analyze the vibration velocity to reflect the transmission of higher-frequency vibration energy in the medium and the impact on the structure, intuitively reflect the intensity of instantaneous vibration, analyze the acceleration to quickly respond to the intensity of instantaneous stratum movement under the action of low-frequency waves, and monitor the acceleration and speed to better analyze the characteristics of mine earthquake stress waves, which can reflect the impact of mining operations on the surrounding In order to avoid the influence of the environment, the collected data is transmitted to the data center in real time through the stress wave monitoring equipment; step 104, the collected stress wave signal is pre-processed, including denoising, filtering and resampling, wherein the collected stress wave signal is denoised to eliminate the interference of non-operating vibration signals of background noise and electromagnetic interference, improve the signal-to-noise ratio of the signal and the accuracy of subsequent analysis, use a digital filter to filter the stress wave signal to remove unnecessary components such as high-frequency noise and low-frequency drift, and resample the filtered signal to unify the sampling rate and time resolution of the signal to ensure subsequent data analysis and processing.

[0008] A further improvement of the technical solution of the present invention is that in step 2, the process of identifying the characteristics of different waveform signals is as follows: step 201, performing characteristic analysis on the pre-processed stress wave signal, identifying the peaks and troughs, arrival time, duration and amplitude in the stress wave signal, identifying the peaks and troughs by detecting the local maximum and minimum values ​​in the signal, using the short-time Fourier transform algorithm to analyze the changes in the signal in the time-frequency domain, accurately picking up the arrival time of each waveform (longitudinal wave, transverse wave, surface wave), calculating the duration according to the start and end time of the waveform, directly reading the amplitude value of the peak or trough, and according to the shape, period, and amplitude characteristics of the waveform, the stress wave signal is identified. The stress wave signal is classified into different waveform types, namely longitudinal wave, transverse wave and surface wave. The longitudinal wave has a faster propagation speed and a shorter period, followed by the transverse wave, and the surface wave has a longer period and a complex waveform. Step 202, using the short-time Fourier transform algorithm to pick up the arrival time of the waveform, using the known distance between the monitoring points and the waveform arrival time, calculate the wave propagation speed, and establish a velocity model based on the propagation speed data of different waveforms (longitudinal wave, transverse wave, surface wave) in different media, analyze the propagation characteristics of stress waves in underground media, so as to determine the change trend of velocity with depth, wherein the expression of the short-time Fourier transform algorithm is: , is at frequency and time The short-time Fourier transform result under is a signal in the time domain, is the window function, which can be a Gaussian window or a Hanning window, is the time offset of the window function, is an imaginary unit, Between 0 and the highest frequency of the signal, is the duration of the signal; the calculation expression of the wave propagation speed is: , is the wave propagation speed, is the distance between two monitoring points, and is the arrival time of the wave at the second and first monitoring points; the calculation expression of the velocity model is: , is the distance from the source to The average speed at The medium is at a distance The speed at is the integrating variable, representing the position in the medium, is the distance from the source to the observation point, It means that as the depth increases, the speed gradually increases; Step 203, analyze the amplitude of the waveform at different distances, record the data, and apply the exponential decay model to describe the relationship between the amplitude decay and the distance to analyze the attenuation law of the amplitude with the distance. By fitting the actual data, the attenuation coefficient is solved. The expression of the exponential decay model is: , is the distance from the source The amplitude at is the initial amplitude at the source, is the attenuation coefficient, which is related to the attenuation characteristics of the medium. is the distance from the source to the observation point, the amplitude With distance Step 204, by comparing the arrival time and waveform shape of the waveform at the same monitoring point, the phase difference between different waveforms is measured, and the relationship between the phase difference and the wave propagation speed, medium density, and elastic modulus physical parameters is analyzed. The change of the phase difference can reflect the heterogeneity and anisotropy of the underground medium. The calculation expression of the phase difference is: , is the phase difference between the two waveforms, and is the phase angle between the two waveforms at corresponding frequencies, Usually in arrive Between, it indicates the phase relationship between the waveforms.

[0009] A further improvement of the technical solution of the present invention is that: in the step 3, the process of establishing a multidimensional velocity model including different formation velocities and a formation model of discontinuity surface inclination is as follows: step 301, collecting existing geological exploration data of the target area, including geological maps, drilling data, geological exploration data, and geophysical logging data, dividing the formations in the target area according to the geological data, analyzing the lithology, thickness, and distribution range of each formation to identify the main fault structure and fold geological structure in the target area, and analyzing the trend, dip, and inclination of the geological structure; step 302, classifying the rocks in the target area into sandstone, mudstone, and limestone according to the physical and mechanical properties of rock density, porosity, and elastic modulus, and drawing a lithology distribution map to determine the spatial distribution characteristics of different lithologies; step 303, determining the longitudinal and transverse wave velocities of each formation according to the geological exploration data and rock physical properties, calculating the formation velocity, and using the geological exploration results to construct a multidimensional velocity model including different formation velocities, and obtaining a multidimensional velocity evaluation index to analyze the formation velocity with increasing depth. The change trend of velocity; step 304, according to the geological structure analysis, determine the discontinuity type of the Moho surface and the Conrad surface in the target area, use the geological exploration data and the fault structure analysis results to determine the dip of the discontinuity, and integrate the discontinuity dip information into the multidimensional velocity model to construct a formation model including the discontinuity dip; step 305, use the established multidimensional velocity model including different formation velocities and the formation model including the discontinuity dip to perform simulation tests, simulate the propagation path, reflection and refraction of stress waves in the formation model, and compare the simulation results with the measured data to verify the accuracy and reliability of the formation model. According to the comparison results, make necessary adjustments and optimizations to the formation model to improve the prediction accuracy of the formation model; step 306, input the pre-processed stress wave signal into the established formation model, use the formation velocity and discontinuity information in the formation model to identify the waveform characteristics of the input stress wave signal, and distinguish between longitudinal waves, transverse waves and surface waves; step 307, calculate the propagation velocity, attenuation law and phase difference parameters of the stress wave according to the waveform characteristics.

[0010] A further improvement of the technical solution of the present invention is that the calculation expression of the formation velocity is: .

[0011] in, is the formation velocity, indicating the The wave velocity of the layer, For the The density of the strata, For the The elastic modulus of the strata, For the Poisson's ratio of the strata, is the number of strata.

[0012] The expression of the multidimensional speed evaluation index is: .

[0013] in, For the depth The formation velocity at For the The thickness of the strata, It represents the limit when the thickness of each layer approaches infinity.

[0014] A further improvement of the technical solution of the present invention is that in step 4, the process of establishing an inversion equation containing formation velocity and discontinuity surface dip parameters is as follows: step 401, extracting the monitored stress wave signal characteristics of wave crests, wave troughs, arrival times, durations and amplitudes, and dividing the formation structure into multiple formations according to geological exploration data, each formation having specific physical properties, setting preliminary physical parameters for each formation, including velocity, density and elastic modulus, to initialize the formation model; step 402, based on the initialized formation model, combining the difference between the forward simulation of the wave equation and the formation model and the observed data, establishing an inversion equation using an objective function method, linking the observed signal characteristics with the formation velocity and the discontinuity surface dip; step 403, parameterizing the formation velocity and the discontinuity surface dip as model parameters, setting initial estimated values ​​for the formation velocity and the discontinuity surface dip, using a genetic algorithm to iteratively solve the inversion equation, and presetting a convergence threshold , adjust the values ​​of formation velocity and discontinuity surface dip to output the final formation velocity and discontinuity surface dip combination as the optimal formation velocity and discontinuity surface dip combination; step 404, analyze the optimal formation velocity and discontinuity surface dip combination searched by the optimization algorithm, verify the consistency between the formation model prediction and the actual observation data, and update the formation model parameters according to the verification results, and output the final formation velocity and discontinuity surface dip model.

[0015] A further improvement of the technical solution of the present invention is that the inversion equation is expressed as: .

[0016] in, is the inversion equation, constructed using the objective function, which represents the total error between the observed data and the calculated data. is the formation velocity, is the discontinuity surface inclination angle, is the number of stress wave data points, is the number of observations corresponding to each data point, For the Data point The arrival time of the observation, For the Data point The calculated arrival time.

[0017] The further improvement of the technical solution of the present invention is that: in the step 5, the process of evaluating the accuracy of the inversion result is as follows: step 501, collating the formation velocity and discontinuity dip angle data obtained by inversion to form a first data set, and extracting corresponding geological exploration data, including borehole logging data, stress wave velocity test results, and geological maps, and screening the formation velocity and discontinuity dip angle data of the same depth point as the inversion result from the geological exploration data to form a second data set; step 502, constructing a comparison table, including depth, inversion formation velocity, exploration formation velocity, inversion discontinuity dip angle and exploration discontinuity dip angle, filling the data in the first data set and the second data set into the comparison table to ensure that the depth corresponds accurately; step 503, comparing the inversion result with the geological exploration data, analyzing the difference in formation velocity between the inversion result and the geological exploration data at different formations and depths, calculating the difference between the inversion formation velocity and the geological exploration velocity, performing statistical analysis, obtaining the error function, calculating the average error rate of the comprehensive data points, and the calculation expression of the difference between the inversion formation velocity and the geological exploration velocity is: , is the difference between the inverted formation velocity and the geological exploration velocity, is the inverted formation velocity, is the velocity in the geological exploration data; step 504, respectively calculating the average values ​​of the inverted formation velocity and the geological exploration velocity, determining the correlation coefficient between the inverted formation velocity and the geological exploration velocity, and analyzing the error trend of the inversion result compared with the geological exploration data, and comprehensively analyzing to evaluate the accuracy of the inversion result, the expression of the error trend analysis is: ,in, is the error trend, the error trend It should be close to zero, indicating that there is no systematic deviation between the inversion result and the actual data.

[0018] A further improvement of the technical solution of the present invention is that the expression of the error function is: .

[0019] in, is the error function, which represents the average error rate, is the number of data points, For the The speed difference of the data points is For the The geological exploration speed of data points, It should be as small as possible and close to zero; the expression of the correlation coefficient between the inverted formation velocity and the geological exploration velocity is: ;

[0020] ;

[0021] .

[0022] in, is the correlation coefficient between the inverted formation velocity and the geological exploration velocity, and are the average values ​​of the inverted formation velocity and geological exploration velocity, is the number of data points, For the The inverted formation velocity of data points is For the The geological exploration speed of each data point.

[0023] The further improvement of the technical solution of the present invention is that: in the step 6, the process of iterative inversion and correction of inversion accuracy and results is as follows: step 601, analyzing the difference in formation velocity between the inversion result and the geological exploration data at different formations and depths, and drawing an error distribution diagram, determining the error between the inversion velocity and the geological exploration velocity, and intuitively displaying the distribution of the error; step 602, adjusting the parameters of the formation model according to the result of the error analysis, including modifying the average velocity, thickness, position and inclination of the formation, so that the inversion result is closer to the geological exploration data; step 603, using the corrected formation model A new inversion calculation is performed on the model. During the inversion process, the changes of the inversion results are continuously monitored and compared with the geological exploration data. An error threshold is set. When the error between the inversion result and the geological exploration data is less than the error threshold, the convergence result of the inversion process is output and the iteration is stopped. Step 604: preset the earthquake source location positioning accuracy, apply the earthquake positioning algorithm of multi-angle measurement, determine the earthquake source location range based on the revised stratigraphic model, repeat the inversion calculation and result verification, until the positioning of the earthquake source location range reaches the preset accuracy, and output the final earthquake source location range inversion result.

[0024] Due to the adoption of the above-mentioned technical scheme, the technical progress achieved by the present invention compared with the prior art is as follows: 1. The present invention provides a mine earthquake source inversion method based on monitoring stress wave signals. Based on the monitored stress wave signal characteristics and in combination with the formation model, an inversion equation including formation velocity and discontinuity surface inclination parameters is established, an optimization algorithm is used to solve the inversion equation, the best combination of formation velocity and discontinuity surface inclination is analyzed, the inverted formation velocity and discontinuity surface inclination are compared with geological exploration data, so as to correct and optimize the formation model, iteratively invert, correct the inversion accuracy and results, and thus determine the location range of the earthquake source.

[0025] 2. The present invention provides a mine earthquake source inversion method based on monitoring stress wave signals. By analyzing the stress wave signals generated by the operation and utilizing advanced signal processing technology and complex numerical simulation, the key geological parameters of the formation velocity, density, elastic modulus and inclination of the discontinuity surface can be accurately identified and quantified. Through an iterative optimization process, the formation model is continuously adjusted and improved until the model prediction is optimally matched with the actual observation data, thereby improving the resolution and accuracy of detection.

[0026] 3. The present invention provides a mine earthquake source inversion method based on monitoring stress wave signals. Through the source inversion method, the stress wave signals generated by the operation are monitored and analyzed, so as to evaluate the operation effect and the impact on the surrounding environment. It can help engineers optimize the operation design, such as adjusting operation parameters, controlling operation risks, and reducing damage to the surrounding environment. Through precise underground structure imaging, it can guide mining operations to avoid key geological structures and reduce the negative impact of operations on underground facilities, water resources and the ecological environment. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0028] Figure 1 The figure is a flow chart of the method of the present invention.

[0029] Figure 2 A flow chart for establishing a multidimensional velocity model including different formation velocities and a formation model of discontinuity dip angles for the present invention.

[0030] Figure 3 A flow chart including the inversion equations of formation velocity and discontinuity dip parameters is established for the present invention.

[0031] Figure 4 Flow chart for evaluating the accuracy of inversion results of the present invention. DETAILED DESCRIPTION

[0032] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are 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.

[0033] Embodiment 1, as Figure 1 , Figure 2As shown, the present invention provides a mine earthquake source inversion method based on monitoring stress wave signals, comprising the following steps: Step 1, arranging a monitoring network, synchronously recording stress wave signals of each monitoring point, and preprocessing the captured stress wave signals, determining the number and location of monitoring points according to the operating area of ​​the mining operation, and installing stress wave monitoring equipment at each monitoring point, including a velocity sensor, an acceleration sensor, and a data acquisition instrument, to ensure accurate capture of stress wave signals, wherein the operating area of ​​the mining operation needs to analyze the specific location of the mining operation, the scope of influence, and the sensitivity of the surrounding environment, so as to determine the number and location of the monitoring points, and the monitoring points fully cover the potentially affected areas First, clarify the specific location of the mining operation, the expected scale of operation and other key information. Based on the scale and nature of the mining operation, combined with geological conditions, topography, building structure and other factors, evaluate the propagation range and impact of the stress waves that may be generated by the operation, and ensure that the monitoring points can fully cover the potentially affected areas, including sensitive buildings, public facilities, residential areas, etc. At the same time, consider setting up several monitoring points near the operation source to obtain more direct vibration data. After the stress wave monitoring equipment is deployed, it is inspected and the equipment parameters are configured. The GPS time synchronization system is used to synchronize the time of all monitoring points to ensure that the subsequent data analysis accurately corresponds to different The signal timestamp of the monitoring point is checked to ensure that the monitoring equipment is fixed firmly to avoid displacement or damage during operation. At the same time, check whether the power supply and communication lines of the equipment are connected normally to ensure that the data can be transmitted to the data center in real time. When mining operations are in progress, stress wave monitoring equipment is used to record and store stress wave signals. Stress wave signals include vibration velocity and acceleration parameters. Analysis of vibration velocity reflects the transmission of higher-frequency vibration energy in the medium and its impact on the structure, which intuitively reflects the intensity of instantaneous vibration. Analysis of acceleration quickly reflects the intensity of instantaneous stratum movement caused by low-frequency waves. Monitoring acceleration and velocity can better analyze the characteristics of mine earthquake stress waves. It can reflect the impact of mining operations on the surrounding environment. The collected data is transmitted to the data center in real time through the stress wave monitoring equipment, and the collected stress wave signals are pre-processed, including denoising, filtering and resampling. Among them, the collected stress wave signals are denoised to eliminate the interference of non-operating vibration signals such as background noise and electromagnetic interference, improve the signal-to-noise ratio of the signal and the accuracy of subsequent analysis, use a digital filter to filter the stress wave signal, remove unnecessary components such as high-frequency noise and low-frequency drift, and resample the filtered signal to unify the sampling rate and time resolution of the signal to ensure subsequent data analysis and processing.

[0034] Step 2, identify the signal characteristics of different waveforms from the preprocessed stress wave signal, and calculate the propagation speed, attenuation law, and phase difference parameters of the stress wave based on the waveform characteristics, perform feature analysis on the preprocessed stress wave signal, identify the peaks and troughs, arrival time, duration, and amplitude in the stress wave signal, identify the peaks and troughs by detecting the local maximum and minimum values ​​in the signal, use the short-time Fourier transform algorithm to analyze the changes in the signal in the time-frequency domain, accurately pick up the arrival time of each waveform (longitudinal wave, transverse wave, surface wave), calculate the duration based on the start and end time of the waveform, directly read the amplitude value of the peak or trough, and calculate the amplitude value based on the shape of the waveform, Based on the period and amplitude characteristics, the stress wave signal is classified into different waveform types, namely longitudinal wave, transverse wave and surface wave. The longitudinal wave has a faster propagation speed and a shorter period, followed by the transverse wave. The surface wave has a longer period and a complex waveform. The short-time Fourier transform algorithm is used to pick up the arrival time of the waveform. The wave propagation speed is calculated using the known distance between the monitoring points and the waveform arrival time. According to the propagation speed data of different waveforms (longitudinal wave, transverse wave, surface wave) in different media, a velocity model is established to analyze the propagation characteristics of stress waves in underground media to determine the velocity change trend with depth. The expression of the short-time Fourier transform algorithm is: , is at frequency and time The short-time Fourier transform result under is a signal in the time domain, is the window function, which can be a Gaussian window or a Hanning window, is the time offset of the window function, is an imaginary unit, Between 0 and the highest frequency of the signal, is the duration of the signal; the calculation expression of the wave propagation speed is: , is the wave propagation speed, is the distance between two monitoring points, and is the arrival time of the wave at the second and first monitoring points; the calculation expression of the velocity model is: , is the distance from the source to The average speed at The medium is at a distance The speed at is the integrating variable, representing the position in the medium, is the distance from the source to the observation point, It means that as the depth increases, the speed gradually increases. The amplitude of the waveform at different distances is analyzed, the data is recorded, and the exponential decay model is used to describe the relationship between the amplitude decay and the distance to analyze the attenuation law of the amplitude with the distance. By fitting the actual data, the attenuation coefficient is solved. The expression of the exponential decay model is: , is the distance from the source The amplitude at is the initial amplitude at the source, is the attenuation coefficient, which is related to the attenuation characteristics of the medium. is the distance from the source to the observation point, the amplitude With distance The phase difference between different waveforms is measured by comparing the arrival time and waveform shape of the waveform at the same monitoring point, and the relationship between the phase difference and the wave propagation speed, medium density, and elastic modulus physical parameters is analyzed. The change of phase difference can reflect the heterogeneity and anisotropy of the underground medium. The calculation expression of phase difference is: , is the phase difference between the two waveforms, and is the phase angle between the two waveforms at corresponding frequencies, Usually in arrive Between, it indicates the phase relationship between the waveforms.

[0035] Step 3: Analyze the stratigraphic structure, lithology distribution and fault structure of the target area in combination with geological exploration data, and establish a multidimensional velocity model containing different stratigraphic velocities and a stratigraphic model of discontinuity dip angles based on geological survey results. Collect existing geological exploration data of the target area, including geological maps, drilling data, geological exploration data, and geophysical logging data. According to the geological data, divide the stratigraphic areas of the target area, analyze the lithology, thickness and distribution range of each stratigraphic area to identify the main fault structures and fold geological structures in the target area, analyze the strike, dip and dip of the geological structure, classify the rocks in the target area into sandstone, mudstone and limestone according to the physical and mechanical properties of rock density, porosity and elastic modulus, and draw a lithology distribution map to determine the spatial distribution characteristics of different lithologies. According to the geological exploration data and rock physical properties, determine the P-wave and S-wave velocities of each stratigraphic area, calculate the stratigraphic velocity, and use the geological exploration results to construct a multidimensional velocity model containing different stratigraphic velocities, obtain a multidimensional velocity evaluation index, and analyze the influence of depth on the stratigraphic velocity. , the change trend of formation velocity, according to the geological structure analysis, determine the discontinuity type of the Moho surface and Conrad surface in the target area, use the geological exploration data and the results of fault structure analysis to determine the dip of the discontinuity, and integrate the discontinuity dip information into the multidimensional velocity model to build a formation model containing the discontinuity dip, use the established multidimensional velocity model containing different formation velocities and the formation model containing the discontinuity dip to carry out simulation tests, simulate the propagation path, reflection and refraction of stress waves in the formation model, and compare the simulation results with the measured data to verify the accuracy and reliability of the formation model. According to the comparison results, make necessary adjustments and optimizations to the formation model to improve the prediction accuracy of the formation model, input the preprocessed stress wave signal into the established formation model, use the formation velocity and discontinuity information in the formation model to identify the waveform characteristics of the input stress wave signal, distinguish between longitudinal waves, shear waves and surface waves, and calculate the propagation velocity, attenuation law and phase difference parameters of the stress wave according to the waveform characteristics; further, the calculation expression of the formation velocity is: ;in, is the formation velocity, indicating the The wave velocity of the strata, For the The density of the layers, For the The elastic modulus of the strata, For the Poisson's ratio of the strata, is the number of strata; the expression of the multidimensional velocity evaluation index is: ;in, For the depth The formation velocity at For the The thickness of the strata, represents the limit when the thickness of each layer tends to infinity, With depth The increase in the formation velocity is linear or nonlinear. When there is little variation between layers, With depth Linear increase or decrease; if the formation velocity There are significant differences between the different layers. With depth The change of shows nonlinear characteristics. If the velocity of a layer is significantly higher or lower than that of the adjacent layer, it will lead to At this depth, a peak or valley occurs; at a geological interface, such as a fault or unconformity, the formation velocity mutations will occur, leading to At this depth, abrupt changes occur; in formations containing multiple rock types, such as alternating layers of sandstone and shale, It fluctuates with depth, reflecting the physical properties of different rock types.

[0036] Step 4: Based on the monitored stress wave signal characteristics and in combination with the formation model, an inversion equation including formation velocity and discontinuity dip parameters is established, and an optimization algorithm is used to solve the inversion equation to analyze the best combination of formation velocity and discontinuity dip.

[0037] Step 5: Compare and verify the inverted formation velocity and discontinuity surface dip with the geological exploration data to evaluate the accuracy of the inversion results.

[0038] Step 6: Modify and optimize the formation model based on the verification results, iterate the inversion, modify the inversion accuracy and results, and determine the location range of the earthquake source.

[0039] Embodiment 2, as Figure 3 , Figure 4 As shown, on the basis of Example 1, the present invention provides a technical solution: preferably, in step 4, the process of establishing the inversion equation containing the formation velocity and discontinuity surface dip parameters is: extracting the monitored stress wave signal characteristics of the peak, trough, arrival time, duration and amplitude, and dividing the formation structure into multiple formations according to the geological exploration data, each formation has specific physical properties, setting preliminary physical parameters for each formation, including velocity, density and elastic modulus, to initialize the formation model, based on the initialized formation model, combining the difference between the forward simulation of the wave equation and the formation model and the observed data, using the objective function method to establish the inversion equation, linking the observed signal characteristics with the formation velocity and the discontinuity surface dip, wherein the expression of the inversion equation is: .

[0040] in, is the inversion equation, constructed using the objective function, which represents the total error between the observed data and the calculated data. is the formation velocity, is the discontinuity surface inclination angle, is the number of stress wave data points, is the number of observations corresponding to each data point, For the Data point The arrival time of the observation, For the Data point The calculated arrival time depends on the formation velocity and the discontinuity angle , The value of decreases gradually with the iteration process and tends to zero. The formation velocity and discontinuity dip are parameterized as model parameters, and initial estimates are set for the formation velocity and discontinuity dip. The inversion equation is solved iteratively using a genetic algorithm, and a convergence threshold is preset. , adjust the values ​​of formation velocity and discontinuity dip angle to output the final combination of formation velocity and discontinuity dip angle as the optimal combination of formation velocity and discontinuity dip angle. The objective function form of the genetic algorithm is: , where minimizing This is achieved through the selection, crossover and mutation operations in the genetic algorithm, specifically:

[0041] Select an action: ;

[0042] Crossover operation: ;

[0043] Mutation operation: .

[0044] in, is the probability of selection, It is a crossover operation used to generate new combinations of formation velocity and discontinuity dip angle. is a mutation operation that introduces random perturbations. is the mutation rate; in each iteration, the formation velocity and discontinuity dip are updated: ,in, represents the number of iterations; when , the algorithm stops, is the preset convergence threshold, and the final expression of the combination of formation velocity and discontinuity surface dip is: , Represents the inversion equation The parameter value to be minimized; the optimal combination of formation velocity and discontinuity dip angle searched by the optimization algorithm is analyzed, the consistency between the formation model prediction and the actual observation data is verified, and the formation model parameters are updated according to the verification results, and the final formation velocity and discontinuity dip angle model is output; in step 5, the process of evaluating the accuracy of the inversion result is as follows: the formation velocity and discontinuity dip angle data obtained by inversion are sorted to form a first data set, and the corresponding geological exploration data, including borehole logging data, stress wave velocity test results, and geological maps, are extracted, and the formation velocity and discontinuity dip angle data at the same depth as the inversion result are screened out from the geological exploration data. According to the data, a second data set is formed, and a comparison table is constructed, including depth, inversion formation velocity, exploration formation velocity, inversion discontinuity dip and exploration discontinuity dip. The data in the first data set and the second data set are filled in the comparison table to ensure that the depth corresponds accurately. The inversion results are compared with the geological exploration data, and the differences in the formation velocities between the inversion results and the geological exploration data at different strata and depths are analyzed. The difference between the inversion formation velocity and the geological exploration velocity is calculated, and statistical analysis is performed to obtain the error function, and the average error rate of the comprehensive data points is calculated. The calculation expression for the difference between the inversion formation velocity and the geological exploration velocity is: , is the difference between the inverted formation velocity and the geological exploration velocity, is the inverted formation velocity, is the velocity in geological exploration data, where the error function is expressed as: .

[0045] in, is the error function, which represents the average error rate, is the number of data points, For the The speed difference of the data points is For the The geological exploration speed of data points, It should be as small as possible, close to zero, and the average values ​​of the inverted formation velocity and the geological exploration velocity are calculated respectively, and the correlation coefficient between the inverted formation velocity and the geological exploration velocity is determined. At the same time, the error trend of the inversion result compared with the geological exploration data is analyzed, and a comprehensive analysis is performed to evaluate the accuracy of the inversion result. The expression of the correlation coefficient between the inverted formation velocity and the geological exploration velocity is:

[0046] ;

[0047] ;

[0048] .

[0049] in, is the correlation coefficient between the inverted formation velocity and the geological exploration velocity, and are the average values ​​of the inverted formation velocity and geological exploration velocity, is the number of data points, For the The inverted formation velocity of data points is For the The geological exploration speed of data points, the correlation coefficient The value range is between -1 and 1. A value close to 1 or -1 indicates a strong correlation, and a value close to 0 indicates no correlation. The expression for error trend analysis is: ,in, is the error trend, the error trend It should be close to zero, indicating that there is no systematic deviation between the inversion result and the actual data; in step 6, the process of iterative inversion and correction of inversion accuracy and results is as follows: analyze the difference in formation velocity between the inversion result and the geological exploration data at different formations and depths, and draw an error distribution diagram to determine the error between the inversion velocity and the geological exploration velocity, and intuitively display the distribution of the error. According to the results of the error analysis, adjust the parameters of the formation model, including modifying the average velocity, thickness, position and inclination of the formation, so that the inversion result is closer to the geological exploration data, and use the corrected formation model for new inversion calculations. In the inversion process, Continuously monitor the changes in the inversion results and compare them with the geological exploration data, set the error threshold, and when the error between the inversion result and the geological exploration data is less than the error threshold, output the convergence result of the inversion process, stop the iteration, and organize all the data generated in the inversion process, including the original data, inversion results, error analysis, and model correction records. Preset the earthquake source location accuracy, apply the earthquake location algorithm of multi-angle measurement, determine the earthquake source location range based on the corrected stratigraphic model, repeat the inversion calculation and result verification until the earthquake source location range is located at the preset accuracy, and output the final earthquake source location range inversion result.

[0050] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art who is familiar with the present technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A mine earthquake source inversion method based on monitoring stress wave signals, characterized in that: The following steps are involved: Step 1: Arrange a monitoring network, synchronously record stress wave signals of various monitoring points, and pre-process the captured stress wave signals; Step 2, identifying signal features of different waveforms from the preprocessed stress wave signal, and calculating the propagation speed, attenuation law, and phase difference parameter of the stress wave according to the waveform features; Step 3, combining geological exploration data, analyzing the stratigraphic structure, lithology distribution and fault structure of the target area, and establishing a multidimensional velocity model including different stratigraphic velocities and a stratigraphic model of discontinuity dip angles according to the geological survey results. In step 3, the process of establishing a multidimensional velocity model including different stratigraphic velocities and a stratigraphic model of discontinuity dip angles is as follows: Step 301, collecting existing geological exploration data of the target area, dividing the strata in the target area, analyzing the lithology, thickness and distribution range of each stratum, identifying the main fault structure and fold structure in the target area, and analyzing the trend, dip and inclination of the geological structure; Step 302, classify the rocks in the target area into sandstone, mudstone, and limestone according to the physical and mechanical properties of the rock such as density, porosity, and elastic modulus, and draw a lithology distribution map to determine the spatial distribution characteristics of different lithologies; Step 303, according to the geological exploration data and rock physical properties, determine the P-wave and S-wave velocities of each formation, calculate the formation velocity, and use the geological exploration results to construct a multidimensional velocity model containing different formation velocities, obtain a multidimensional velocity evaluation index, and analyze the change trend of the formation velocity with increasing depth; Step 304, based on the geological structure analysis, determine the discontinuity type of the Moho surface and the Conrad surface in the target area, use the geological exploration data and the fault structure analysis results to determine the dip of the discontinuity, and integrate the discontinuity dip information into the multidimensional velocity model to construct a stratigraphic model including the discontinuity dip; Step 305, using the established multi-dimensional velocity model including different formation velocities and the formation model including the discontinuity surface dip angle to perform simulation testing, simulate the propagation path of the stress wave in the formation model, and compare the simulation results with the measured data, and adjust and optimize the formation model according to the comparison results; Step 306, input the pre-processed stress wave signal into the established formation model, and use the formation velocity and discontinuity information in the formation model to identify the waveform characteristics of the input stress wave signal to distinguish between longitudinal waves, shear waves and surface waves; Step 307, calculate the propagation velocity, attenuation law, and phase difference parameter of the stress wave according to the waveform characteristics. The calculation expression of the formation velocity is: ； in, is the formation velocity, indicating the The wave velocity of the strata, For the The density of the layers, For the The elastic modulus of the strata, For the Poisson's ratio of the strata, is the number of strata; The expression of the multidimensional speed evaluation index is: ； in, For the depth The formation velocity at For the The thickness of the strata, It represents the limit when the thickness of each layer approaches infinity; Step 4, based on the monitored stress wave signal characteristics, combined with the formation model, establish an inversion equation containing formation velocity and discontinuity dip parameters, and use an optimization algorithm to solve the inversion equation to analyze the best combination of formation velocity and discontinuity dip. In step 4, the process of establishing the inversion equation containing formation velocity and discontinuity dip parameters is as follows: Step 401, extracting the monitored stress wave signal characteristics of wave crest, wave trough, arrival time, duration and amplitude, and dividing the stratum structure into multiple strata according to geological exploration data, setting preliminary physical parameters for each stratum, including velocity, density and elastic modulus, to initialize the stratum model; Step 402, based on the initialized formation model, combining the difference between the forward modeling of the wave equation and the formation model and the observed data, an inversion equation is established by using an objective function method to relate the observed signal characteristics to the formation velocity and the discontinuity surface dip; Step 403: parameterize the formation velocity and discontinuity dip angle as model parameters, set initial estimated values ​​for the formation velocity and discontinuity dip angle, use a genetic algorithm to iteratively solve the inversion equation, and preset a convergence threshold. , adjusting the values ​​of formation velocity and discontinuity surface dip angle to output a final combination of formation velocity and discontinuity surface dip angle as the optimal combination of formation velocity and discontinuity surface dip angle; Step 404, analyze the best combination of formation velocity and discontinuity dip angle searched by the optimization algorithm, verify the consistency between the formation model prediction and the actual observation data, update the formation model parameters according to the verification results, and output the final formation velocity and discontinuity dip angle model. The inversion equation is expressed as: ； in, is the inversion equation, constructed using the objective function, which represents the total error between the observed data and the calculated data. is the formation velocity, is the discontinuity surface inclination angle, is the number of stress wave data points, is the number of observations corresponding to each data point, For the Data point The arrival time of the observation, For the Data point The calculated arrival time; Step 5, comparing and verifying the formation velocity and discontinuity surface dip obtained by inversion with geological exploration data to evaluate the accuracy of the inversion results; Step 6: Modify and optimize the formation model based on the verification results, iterate the inversion, modify the inversion accuracy and results, and determine the location range of the earthquake source.

2. The method for inversion of mine earthquake sources based on monitoring stress wave signals according to claim 1 is characterized in that: In step 1, the process of acquiring the stress wave signal is as follows: Step 101, according to the operation area of ​​the mining operation, determine the number and location of the monitoring points, and install stress wave monitoring equipment at each monitoring point, including a velocity sensor, an acceleration sensor and a data acquisition instrument; Step 102, after the stress wave monitoring equipment is deployed, it is inspected, and equipment parameters are configured, and the GPS time synchronization system is used to synchronize the time of all monitoring points; Step 103: When mining operations are in progress, stress wave monitoring equipment is used to record and store stress wave signals. Stress wave signals include vibration velocity and acceleration parameters. Analysis of vibration velocity reflects the transmission of higher-frequency vibration energy in the medium and its impact on the structure, which directly reflects the intensity of instantaneous vibration. Analysis of acceleration quickly reflects the intensity of instantaneous stratum movement under the action of low-frequency waves. Monitoring acceleration and velocity can better analyze the characteristics of mine earthquake stress waves. The collected data is transmitted to the data center in real time through the stress wave monitoring equipment. Step 104, preprocessing the collected stress wave signal, including denoising, filtering and resampling, wherein the collected stress wave signal is denoised to eliminate background noise, electromagnetic interference and non-operating vibration signal interference, the stress wave signal is filtered using a digital filter, and the filtered signal is resampled to unify the sampling rate and time resolution of the signal.

3. The method for inversion of mine earthquake sources based on monitoring stress wave signals according to claim 2 is characterized in that: In step 2, the process of identifying the characteristics of different waveform signals is as follows: Step 201, perform feature analysis on the pre-processed stress wave signal, identify the peaks and troughs, arrival time, duration and amplitude in the stress wave signal, and analyze the stress wave signal according to the shape, period and amplitude characteristics of the waveform. Classify stress wave signals into different waveform types, namely longitudinal waves, transverse waves and surface waves; Step 202, using the short-time Fourier transform algorithm to pick up the arrival time of the waveform, using the known distance between the monitoring points and the waveform arrival time, calculate the wave propagation velocity, and establish a velocity model based on the propagation velocity data of different waveforms in different media, analyze the propagation characteristics of stress waves in underground media, and determine the velocity change trend with depth. The expression of the short-time Fourier transform algorithm is: , is at frequency and time The short-time Fourier transform result under is a signal in the time domain, is the window function, is the time offset of the window function, is an imaginary unit; Step 203, analyze the amplitude of the waveform at different distances, record the data, apply the exponential decay model to describe the relationship between the amplitude decay and the distance, so as to analyze the attenuation law of the amplitude with the distance, and solve the attenuation coefficient by fitting the actual data. The expression of the exponential decay model is: , is the distance from the source The amplitude at is the initial amplitude at the source, is the attenuation coefficient, is the distance from the source to the observation point, the amplitude With distance The increase of decreases exponentially; Step 204, by comparing the arrival time and waveform shape of the waveform at the same monitoring point, the phase difference between different waveforms is measured, and the relationship between the phase difference and the wave propagation speed, medium density, and elastic modulus physical parameters is analyzed.

4. The method for inversion of mine earthquake sources based on monitoring stress wave signals according to claim 3 is characterized in that: In step 5, the process of evaluating the accuracy of the inversion result is: Step 501, sorting the formation velocity and discontinuity dip angle data obtained by inversion to form a first data set, and extracting corresponding geological exploration data, including borehole logging data, stress wave velocity test results, and geological maps, and screening out the formation velocity and discontinuity dip angle data of the same depth point as the inversion result from the geological exploration data to form a second data set; Step 502, constructing a comparison table, including depth, inversion formation velocity, geological exploration velocity, inversion discontinuity dip angle and exploration discontinuity dip angle, and filling the data in the first data set and the second data set into the comparison table; Step 503, compare the inversion result with the geological exploration data, analyze the difference between the inversion result and the formation velocity of the geological exploration data in different formations and depths, calculate the difference between the inverted formation velocity and the geological exploration velocity, perform statistical analysis, obtain the error function, calculate the average error rate of the comprehensive data points, and the calculation expression of the difference between the inverted formation velocity and the geological exploration velocity is: , is the difference between the inverted formation velocity and the geological exploration velocity, is the inverted formation velocity, is the geological exploration speed; Step 504, respectively calculate the average values ​​of the inverted formation velocity and the geological exploration velocity, determine the correlation coefficient between the inverted formation velocity and the geological exploration velocity, and analyze the error trend of the inversion result compared with the geological exploration data, and conduct a comprehensive analysis to evaluate the accuracy of the inversion result.

5. The method for inversion of mine earthquake sources based on monitoring stress wave signals according to claim 4 is characterized in that: The expression of the error function is: ； in, G is the error function, which represents the average error rate, is the number of stress wave data points, For the The speed difference of the data points is For the The geological exploration speed of each data point; The expression of the correlation coefficient between the inverted formation velocity and the geological exploration velocity is: ； ； ； in, is the correlation coefficient between the inverted formation velocity and the geological exploration velocity, and are the average values ​​of the inverted formation velocity and geological exploration velocity, is the number of stress wave data points, For the The inverted formation velocity of data points is For the The geological exploration speed of each data point.

6. The method for inversion of mine earthquake sources based on monitoring stress wave signals according to claim 5 is characterized in that: In step 6, the process of iterative inversion and correction of inversion accuracy and results is as follows: Step 601, analyzing the difference in formation velocity between the inversion result and the geological exploration data at different formations and depths, and drawing an error distribution diagram; Step 602, adjusting the parameters of the formation model according to the results of the error analysis, including modifying the average velocity, thickness, position and dip of the formation, so that the inversion result is close to the geological exploration data; Step 603, using the modified formation model to perform a new inversion calculation, during the inversion process, continuously monitoring the changes in the inversion results, and comparing them with the geological exploration data, setting an error threshold, and when the error between the inversion result and the geological exploration data is less than the error threshold, outputting the inversion process convergence result and stopping the iteration; Step 604, preset the earthquake source location positioning accuracy, apply the earthquake positioning algorithm of multi-angle measurement, determine the earthquake source location range based on the corrected stratigraphic model, repeat the inversion calculation and result verification until the earthquake source location range positioning reaches the preset accuracy, and output the final earthquake source location range inversion result.

Images