A method for evaluating roof blasting decompression effect based on dual-source shock wave CT inversion

Abstract

The present invention discloses a method for evaluating the pressure relief effect of roof blasting based on dual-source shock wave CT inversion. First, a CT detection inversion area is determined according to the mining conditions of the working face, and the passive sources before the roof blasting are screened and the inversion conditions are judged. The passive sources that meet the conditions are CT inverted to obtain the wave velocity distribution before the roof blasting; then the roof blasting pressure relief data is recorded in the wave velocity anomaly area, the passive sources after the blasting are positioned by using a relative positioning algorithm, and then the active and passive sources are subjected to active and passive dual-source CT inversion; finally, according to the CT inversion results before and after the roof blasting, the wave velocity variation coefficient of the blasting area is calculated to evaluate the pressure relief effect of the roof blasting; this method can not only evaluate the pressure relief effect of the roof blasting area in an overall and intuitive manner, but also, due to the addition of the active source, the passive source data required for the inversion is greatly reduced, while ensuring the inversion accuracy, the inversion period is also greatly shortened, and the method has strong universality and good application value.

Description

Technical Field

[0001] The invention relates to a roof blasting pressure relief effect evaluation method based on dual-source shock wave CT inversion, and belongs to the technical field of coal mining and coal mine safety detection. Background Art

[0002] Rock burst is a dynamic disaster caused by the sudden release of energy accumulated in the surrounding rock of a roadway or in the coal and rock mass around the mining face. With the gradual increase in the mining depth and mining intensity of coal mines in my country, rock burst disasters have become more serious. Among them, the hard and thick roof has always been an important influencing factor of rock burst. The roof blasting pressure relief measure is to destroy the rock mass through the shock waves and stress waves generated by blasting to create cracks, weaken the hard roof rock structure and reduce the stress concentration of the working face. However, how to evaluate the pressure relief effect of roof blasting is still relatively difficult. Currently, the commonly used methods for evaluating the pressure relief effect mainly include borehole observation and microseismic data comparison. Among them, borehole observation requires the arrangement of a large number of detection drilling holes, which is a large construction project. In addition, borehole observation can only observe the expansion of cracks in a single borehole and cannot make a holistic evaluation of the stress changes in the pressure relief area. Microseismic data comparison can only analyze the changing trend of microseisms before and after blasting pressure relief and cannot intuitively evaluate the pressure relief effect of the blasting area.

[0003] Therefore, there is an urgent need for an evaluation method for the decompression effect of coal mine roof blasting to achieve a holistic and intuitive assessment of stress changes in the decompression area. Shockwave CT inversion technology inverts the wave velocity distribution along the propagation path based on the distance between the source and the sensor and the initial motion and travel time of the shockwave, and determines the stress distribution in the inversion area based on the correlation between wave velocity and stress, which can better meet this need. However, this technology requires a sufficient number of shockwave rays covering the detection area to achieve a relatively accurate evaluation. In addition, due to the randomness of the location of microseismic events induced by working face production, the source location calculated using the arrival time of the microseismic event has a certain error from the actual occurrence location. Therefore, using only traditional passive source CT inversion of microseismic events induced by working face production requires a sufficient number of microseismic events in the inversion area, resulting in a long inversion period and large errors in the inversion results.

[0004] In summary, how to adopt a new method to quickly realize the effect evaluation of roof blasting pressure relief and ensure the accurate inversion of the roof blasting area is the research direction required by the present invention. Summary of the Invention

[0005] In response to the problems existing in the above-mentioned prior art, the present invention provides a method for evaluating the effect of roof blasting unloading based on dual-source shock wave CT inversion, which can not only quickly realize the effect evaluation of roof blasting unloading, but also ensure the accurate inversion of the roof blasting area.

[0006] To achieve the above objectives, the present invention adopts a technical solution: a method for evaluating the roof blasting decompression effect based on dual-source shock wave CT inversion, which specifically comprises the following steps:

[0007] S1. Several microseismic sensors are arranged in the tunnels on both sides of the working face. Each microseismic sensor receives microseismic waveforms generated by active and passive sources on the working face and records them through a microseismic monitoring system. The active source is the seismic source induced by roof blasting, and the passive source is the microseismic event induced by production on the working face.

[0008] S2. Determine the CT detection inversion area based on the mining situation of the working face, screen the passive sources before roof blasting received in the detection inversion area, and identify the inversion conditions. Finally, use the known shock wave velocity inversion method to perform CT inversion on the passive sources that meet the requirements;

[0009] S3. Obtain the wave velocity distribution of the initial state of the working face before blasting according to step S2, perform roof blasting to relieve pressure in the wave velocity abnormality area, and record the position and time of the roof blasting and the waveforms received by each microseismic sensor;

[0010] S4. After the roof blasting is relieved, the working face continues to be mined, and the microseismic waveforms of the passive sources after the roof blasting are recorded. The active source is used as the reference source, and the passive source is used as the calculation source. The relative positioning algorithm is used to calculate the location of each passive source.

[0011] S5. Screen active and passive sources and identify inversion conditions, and use known shock wave velocity inversion methods to perform active and passive dual-source CT inversion on active and passive sources that meet the requirements;

[0012] S6. Based on the CT inversion results of the detection inversion area before and after the roof blasting, the wave velocity before and after the blasting is compared and analyzed, and the wave velocity change coefficient of the detection inversion area is calculated. Evaluate the decompression effect of roof blasting; wave velocity variation coefficient A value greater than 0 indicates that the pressure relief is effective, and the larger the value, the better the pressure relief effect; Wave velocity variation coefficient Less than 0 indicates that the stress is transferred to this area after pressure relief, increasing the degree of stress concentration; the specific formula is:

[0013]

[0014] Where, is the wave velocity variation coefficient, V p1 is the wave velocity value of the inversion area before blasting, V p2 is the wave velocity value of the detection inversion area after blasting.

[0015] Furthermore, in step S2, the passive sources before roof blasting received in the detection inversion area are screened and the inversion conditions are determined. The specific process is:

[0016] Based on microseismic positioning, the propagation paths between each microseismic event and the microseismic sensor are drawn. The source whose propagation path passes through the detection inversion area is the effective source. The propagation paths between all effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded (generally the grid spacing is 10 to 30 meters). The inversion conditions are met when the number of rays in each unit grid is greater than 60.

[0017] Furthermore, in step S4, a relative positioning algorithm is used to perform positioning calculations on each passive source, specifically:

[0018] Assuming that the parameters of the passive source to be located are W(t0, x0, y0, z0), the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W is:

[0019] t i =t0+T(D i ,V P )+ΔT(D i ,V P )+ε i (1)

[0020] Where, t i is the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W; T(D i ,V P ) is based on the P-wave velocity model V P Calculated theoretical travel time; ΔT(D i ,V P ) is the abnormal value of the theoretical travel time; ε i is the arrival reading error according to the normal distribution;

[0021] The parameters of the active source are known to be R(t R0 ,x R ,y R ,z R ), then the first arrival time of the P wave detected by the i-th microseismic sensor of the active source R is:

[0022] t Ri =t R0 +T(D Ri ,V PR )+ΔT(D Ri ,V PR )+ε Ri (2)

[0023] Using the Taylor formula to perform a first-order Taylor expansion on equation (1), the first arrival time of the P wave of the passive source W is approximately:

[0024]

[0025] Where δt0 and δx0, δy0, δz0 are the correction values ​​of the earthquake occurrence time and microseismic source coordinates relative to the active source R; the active source R and the passive source W pass through the same coal and rock mass, and the reading errors are similar, then subtracting equation (2) from equation (3) yields:

[0026]

[0027] The matrix form of formula (4) is:

[0028] δt=Aδθ (5)

[0029] In the above formula, δt is the element matrix of n-dimensional column vector δt, δt i =t i -t Ri ; δθ is the column matrix of the passive source W source parameters [t0 x0 y0 z0] T ; Matrix A is a (n×4) partial differential matrix:

[0030]

[0031] The matrix A is determined by differentiating the travel time of the active source R coordinate:

[0032]

[0033] In the formula, (x R ,y R ,z R ) is the active source coordinate; (x i ,y i ,z i ) is the coordinate of the i-th microseismic sensor; for the active source, the P-wave velocity function V P i use d Ri / (t Ri -t R0 ), then the matrix A is:

[0034]

[0035] Mode

[0036] Finally, according to Equations (5) and (8), given δt and the partial differential matrix A, the least squares method can be used to solve the seismic time t0 and the source coordinates (x0, y0, z0) of the passive source.

[0037] Furthermore, in step S5, the active sources and the passive sources are screened and the inversion conditions are determined. The specific process is as follows:

[0038] Based on the blasting positions of the active sources and the final positioning results of the passive sources, the propagation paths between each active source, passive source and microseismic sensor are drawn. The sources whose propagation paths pass through the detection inversion area are considered effective sources. The propagation paths between all effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded. The inversion conditions are met when the number of rays in each unit grid is greater than 60.

[0039] Compared with the prior art, the present invention first determines the CT detection inversion area according to the working face mining situation, and screens and identifies the inversion conditions of the working face production-induced microseismic events (i.e., passive sources) before the roof blasting, and performs passive CT inversion on the passive sources that meet the conditions to obtain the wave velocity distribution before the roof blasting; then, the roof blasting is performed to unload the wave velocity anomaly area, and the blasting time, microseismic events, and the vibration wave waveform received by the sensor are recorded. According to the position, time, and waveform data of the roof blasting, a relative positioning algorithm is adopted, and the active source is used as the reference source to locate the passive source after the blasting; then, active and passive dual-source CT inversion is performed on the roof blasting-induced source (i.e., active source) and the working face production-induced microseismic events (i.e., passive source) after the blasting; finally, according to the CT inversion results before and after the roof blasting, the wave velocity change is compared and analyzed, the wave velocity change coefficient of the blasting area is calculated, and the unloading effect of the roof blasting is evaluated. The present invention utilizes roof blasting-induced microseismicity and working face production-induced microseismicity to perform active and passive dual-source CT inversion to evaluate the roof blasting decompression effect. This not only allows for a holistic and intuitive evaluation of the decompression effect of the roof blasting zone, but also, due to the inclusion of the blasting source (i.e., the active source), significantly reduces the amount of microseismic event data (i.e., the passive source) induced by working face production required for inversion. This significantly shortens the inversion cycle while ensuring inversion accuracy. Furthermore, the present invention utilizes existing microseismic data and active source data generated during roof blasting decompression, eliminating the need for additional drilling peepholes and artificial excitation sources. This results in low investment costs, high accuracy in decompression effect evaluation, and strong universality, making it highly applicable. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 It is the overall flow chart of the present invention;

[0041] Figure 2 The working face advancement and production-induced microseismic conditions before roof blasting in an embodiment of the present invention;

[0042] Figure 3 This is an inversion diagram of the shock wave velocity before roof blasting according to an embodiment of the present invention;

[0043] Figure 4 The roof blasting event and microseismic sensor arrangement of the embodiment of the present invention;

[0044] Figure 5 The working face advancement and production-induced microseismic conditions after roof blasting in an embodiment of the present invention;

[0045] Figure 6 This is an inversion diagram of the active and passive dual-source shock wave velocity after roof blasting according to an embodiment of the present invention;

[0046] Figure 7 This is a diagram of the velocity variation coefficient of the shock wave before and after the roof blasting according to an embodiment of the present invention. DETAILED DESCRIPTION

[0047] The present invention will be further described below.

[0048] like Figure 1 As shown, the roof blasting pressure relief effect of the 802 working face of a coal mine was evaluated using this embodiment. The specific steps are as follows:

[0049] S1. Several microseismic sensors are arranged in the tunnels on both sides of the working face. Each microseismic sensor receives microseismic waveforms generated by active and passive sources on the working face and records them through a microseismic monitoring system. The active source is the seismic source induced by roof blasting, and the passive source is the microseismic event induced by production on the working face.

[0050] S2. Determine the CT detection inversion area according to the mining situation of the working face. Figure 2 As shown in the figure, the passive sources before roof blasting received in the detection inversion area are screened and the inversion conditions are judged. The specific process is as follows:

[0051] According to the microseismic positioning, the propagation paths between each microseismic event and the microseismic sensor are drawn. The source whose propagation path passes through the detection inversion area is the effective source; the propagation paths between all the effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded with a grid spacing of 10 to 30 meters. When the number of rays in each unit grid is greater than 60, the inversion condition is met; finally, the known shock wave velocity inversion method is used to perform CT inversion on the passive source that meets the requirements, such as Figure 3 As shown;

[0052] S3, according to step S2, obtain the wave velocity distribution of the initial state before the blasting of the working face, perform roof blasting to relieve pressure in the wave velocity abnormality area, and record the position and time of the roof blasting and the waveform received by each microseismic sensor. Figure 4 As shown;

[0053] S4. After the roof blasting is released, the working face continues to be mined. The microseismic waveform of the passive source after the roof blasting is recorded as follows: Figure 5As shown; the active source is used as the reference source, the passive source is used as the calculation source, and the relative positioning algorithm is used to calculate the positioning of each passive source, specifically:

[0054] Assuming that the parameters of the passive source to be located are W(t0, x0, y0, z0), the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W is:

[0055] t i =t0+T(D i ,V P )+ΔT(D i ,V P )+ε i (1)

[0056] Where, t i is the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W; T(D i ,V P ) is based on the P-wave velocity model V P Calculated theoretical travel time; ΔT(D i ,V P ) is the abnormal value of the theoretical travel time; ε i is the arrival reading error according to the normal distribution;

[0057] The parameters of the active source are known to be R(t R0 ,x R ,y R ,z R ), then the first arrival time of the P wave detected by the i-th microseismic sensor of the active source R is:

[0058] t Ri =t R0 +T(D Ri ,V PR )+ΔT(D Ri ,V PR )+ε Ri (2)

[0059] Using the Taylor formula to perform a first-order Taylor expansion on equation (1), the first arrival time of the P wave of the passive source W is approximately:

[0060]

[0061] Where δt0 and δx0, δy0, δz0 are the correction values ​​of the earthquake occurrence time and microseismic source coordinates relative to the active source R; the active source R and the passive source W pass through the same coal and rock mass, and the reading errors are similar, then subtracting equation (2) from equation (3) yields:

[0062]

[0063] The matrix form of formula (4) is:

[0064] δt=Aδθ (5)

[0065] In the above formula, δt is the element matrix of n-dimensional column vector δt, δt i =t i -t Ri ; δθ is the column matrix of the passive source W source parameters [t0 x0 y0 z0] T ; Matrix A is a (n×4) partial differential matrix:

[0066]

[0067] The matrix A is determined by differentiating the travel time of the active source R coordinate:

[0068]

[0069] In the formula, (x R ,y R ,z R ) is the active source coordinate; (x i ,y i ,z i ) is the coordinate of the i-th microseismic sensor; for the active source, the P-wave velocity function V P i use d Ri / (t Ri -t R0 ), then the matrix A is:

[0070]

[0071] Mode

[0072] Finally, according to Equations (5) and (8), given δt and the partial differential matrix A, the least squares method can be used to solve the seismic time t0 and the source coordinates (x0, y0, z0) of the passive source.

[0073] S5. Screening the active and passive sources and determining the inversion conditions. The specific process is as follows: based on the blasting position of the active source and the final positioning result of the passive source, the propagation paths between each active source, passive source and microseismic sensor are drawn. The source whose propagation path passes through the detection inversion area is the effective source; the propagation paths between all the effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded. When the number of rays in each unit grid is greater than 60, the inversion condition is met;

[0074] Then, the known shock wave velocity inversion method is used to perform active and passive dual-source CT inversion on the active and passive sources that meet the requirements; Figure 6As shown in the figure, after the roof blasting and pressure relief, the high wave velocity area in the two lanes of the 802 working face shifted to the inside of the working face, reducing the possibility of rock burst in the two lanes of the working face.

[0075] S6. Based on the CT inversion results of the detection inversion area before and after the roof blasting, the wave velocity before and after the blasting is compared and analyzed, and the wave velocity change coefficient of the detection inversion area is calculated. Evaluate the decompression effect of roof blasting; wave velocity variation coefficient A value greater than 0 indicates that the pressure relief is effective, and the larger the value, the better the pressure relief effect; Wave velocity variation coefficient Less than 0 indicates that the stress is transferred to this area after pressure relief, increasing the degree of stress concentration; the specific formula is:

[0076]

[0077] Where, is the wave velocity variation coefficient, V p1 is the wave velocity value of the inversion area before blasting, V p2 is the wave velocity value of the detection inversion area after blasting.

[0078] Figure 7 The velocity variation coefficient of the vibration wave before and after the roof blasting is shown in the figure. It can be seen from the figure that the velocity variation coefficient of the two lanes in the working face is is greater than 0, and the velocity variation coefficient in the tunnel area of ​​the working face is the largest, indicating that the roof blasting in the two tunnel areas has a good pressure relief effect; the velocity variation coefficient in the middle of the working face is It is less than 0, indicating that the stress in the two tunnel areas is transferred to the middle of the working face after the roof blasting is used to relieve pressure, which increases the stress concentration in the middle of the working face, that is, the wave velocity in the middle of the working face increases.

[0079] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for evaluating roof blasting decompression effect based on dual-source shock wave CT inversion, characterized in that: The specific steps are: S1. Several microseismic sensors are arranged in the tunnels on both sides of the working face. Each microseismic sensor receives microseismic waveforms generated by active and passive sources on the working face and records them through a microseismic monitoring system. The active source is the seismic source induced by roof blasting, and the passive source is the microseismic event induced by production on the working face. S2. Determine the CT detection inversion area based on the mining situation of the working face, screen the passive sources before roof blasting received in the detection inversion area, and identify the inversion conditions. Finally, use the shock wave velocity inversion method to perform CT inversion on the passive sources that meet the requirements. S3. Obtain the wave velocity distribution of the initial state of the working face before blasting according to step S2, perform roof blasting to relieve pressure in the wave velocity abnormality area, and record the position and time of the roof blasting and the waveforms received by each microseismic sensor; S4. After the roof blasting is relieved, the working face continues to be mined, and the microseismic waveforms of the passive sources after the roof blasting are recorded. The active source is used as the reference source, and the passive source is used as the calculation source. The relative positioning algorithm is used to calculate the location of each passive source. S5. Screen active and passive sources and identify inversion conditions, and use the vibration wave velocity inversion method to perform active and passive dual-source CT inversion on active and passive sources that meet the requirements; S6. Based on the CT inversion results of the detection inversion area before and after the roof blasting, the wave velocity before and after the blasting is compared and analyzed, and the wave velocity change coefficient of the detection inversion area is calculated. Evaluate the decompression effect of roof blasting; wave velocity variation coefficient A value greater than 0 indicates that the pressure relief is effective, and the larger the value, the better the pressure relief effect; Wave velocity variation coefficient Less than 0 indicates that the stress is transferred to this area after pressure relief, increasing the degree of stress concentration; the specific formula is: Where, is the wave velocity variation coefficient, V p1 is the wave velocity value of the inversion area before blasting, V p2 is the wave velocity value of the detection inversion area after blasting.

2. The roof blasting decompression effect evaluation method based on dual-source shock wave CT inversion according to claim 1 is characterized in that: In step S2, the passive sources before roof blasting received in the detection inversion area are screened and the inversion conditions are determined. The specific process is: Based on microseismic positioning, the propagation paths between each microseismic event and the microseismic sensor are drawn. The source whose propagation path passes through the detection inversion area is considered to be the effective source. The propagation paths between all effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded. The inversion conditions are met when the number of rays in each unit grid is greater than 60.

3. The roof blasting decompression effect evaluation method based on dual-source shock wave CT inversion according to claim 1 is characterized in that: In step S4, a relative positioning algorithm is used to calculate the positioning of each passive source, specifically: Assuming that the parameters of the passive source to be located are W(t0, x0, y0, z0), the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W is: t i =t0+T(D i ,V P )+ΔT(D i ,V P )+ε i (1) Where, t i is the first arrival time of the P wave detected by the i-th microseismic sensor of the passive source W; T(D i ,V P ) is based on the P-wave velocity model V P Calculated theoretical travel time; ΔT(D i ,V P ) is the abnormal value of the theoretical travel time; ε i is the arrival reading error according to the normal distribution; The parameters of the active source are known to be R(t R0 ,x R ,y R ,z R ), then the first arrival time of the P wave detected by the i-th microseismic sensor of the active source R is: t Ri =t R0 +T(D Ri ,V PR )+ΔT(D Ri ,V PR )+ε Ri (2) Using the Taylor formula to perform a first-order Taylor expansion on equation (1), the first arrival time of the P wave of the passive source W is approximately: Where δt0 and δx0, δy0, δz0 are the correction values ​​of the earthquake occurrence time and microseismic source coordinates relative to the active source R; the active source R and the passive source W pass through the same coal and rock mass, and the reading errors are similar, then subtracting equation (2) from equation (3) yields: The matrix form of formula (4) is: δt=Aδθ (5) In the above formula, δt is the element matrix of n-dimensional column vector δt, δt i =t i -t Ri ; δθ is the column matrix of the passive source W source parameters [t0x0 y0 z0] T ; Matrix A is a (n×4) partial differential matrix: The matrix A is determined by differentiating the travel time of the active source R coordinate: In the formula, (x R ,y R ,z R ) is the active source coordinate; (x i ,y i ,z i ) is the coordinate of the i-th microseismic sensor; for the active source, the P-wave velocity function V P i use d Ri / (t Ri -t R0 ), then the matrix A is: Mode Finally, according to Equations (5) and (8), given δt and the partial differential matrix A, the least squares method can be used to solve the seismic time t0 and the source coordinates (x0, y0, z0) of the passive source.

4. The roof blasting decompression effect evaluation method based on dual-source shock wave CT inversion according to claim 1 is characterized in that: In step S5, the active sources and the passive sources are screened and the inversion conditions are determined. The specific process is as follows: Based on the blasting positions of the active sources and the final positioning results of the passive sources, the propagation paths between each active source, passive source and microseismic sensor are drawn. The sources whose propagation paths pass through the detection inversion area are considered effective sources. The propagation paths between all effective sources and each microseismic sensor are formed into a ray network, and the detection inversion area is gridded. The inversion conditions are met when the number of rays in each unit grid is greater than 60.

Images