Multi-source earthquake gas advanced detection method for underground coal mine gas enrichment area

Abstract

The invention discloses a multi-source earthquake gas advanced detection method for an underground coal mine gas enrichment area, which can monitor the dynamic change of gas in real time by adopting mining-following earthquake data, and can discover the abnormal change of gas parameters in time through continuous earthquake data acquisition and analysis to provide early warning for coal mine safety. The real-time monitoring is helpful for the mine to take measures in time to prevent the occurrence of gas disasters. Compared with a traditional gas monitoring method, direct contact with gas or coal is not needed, and therefore the safety risk in the detection process is reduced. According to the method, gas advanced detection is achieved, meanwhile, safety of detection personnel is protected, possible interference to the gas occurrence state is avoided, the method has the multiple advantages of being high in precision, real-time, multi-parameter comprehensive analysis, non-contact detection, wide in applicability and the like, powerful technical guarantee is provided for safe production of a coal mine, and the method is suitable for popularization and application. The method is suitable for large-scale industrial use and popularization.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical exploration, and particularly relates to a multi-source seismic gas advanced detection method for gas enrichment areas in coal mines underground. Background Technique

[0002] With the gradual reduction of shallow coal resources, coal mining activities are increasingly developing towards the deep. Deep rock masses have the characteristics of "three highs and strong disturbances", namely high ground stress, high ground temperature, high karst water pressure, and strong disturbance phenomena during the mining process. Compared with shallow mining, deep mining activities are more easily troubled by problems such as coal and gas outbursts, rock bursts, and large deformations of roadway surrounding rocks. Effective identification and accurate early warning of coal mine gas, rock bursts, water inrusions and other disasters in all aspects are the key measures to contain these geological disasters. The continuous increase in the mining depth of the mine has led to an upward trend in the coal seam gas content and gas pressure, and thus the risk of coal and gas outburst accidents has also increased.

[0003] There are various related factors for the occurrence of coal and gas outbursts. Most outbursts occur in geological structure zones or at positions where the coal seam changes sharply, but not all geological structure zones or sharp changes in the coal seam will result in outburst accidents. Therefore, the structure zone is only an important risk factor for the occurrence of outburst accidents. Usually, under the in-situ state of high pressure, gas exists in the coal seam in an adsorbed state. After the coal seam is depressurized due to mining activities, the gas will change from the adsorbed state to the free state and be released. For coal seams with good integrity, it is not easy for gas to depressurize and the release is slow. However, once the broken coal seam is depressurized, gas is easily released rapidly, causing serious accidents of gas exceeding the limit.

[0004] In the prior art, the early warning methods and systems for gas outbursts and their compound dynamic disasters are mainly developed based on electromagnetic radiation, sound and electricity, microseismic, and gas concentration time series monitoring technologies. Due to differences in coal seam occurrence, mining conditions, ground stress, etc., the critical values of the early warning system index system are determined through experience, theory or laboratory, resulting in low accuracy of early warning results. The workload of field correction of the critical values of the index system is large, and the regional differences directly affect the on-site application effect of the early warning system. A large number of borehole explorations are time-consuming and laborious, and it is only "a view from one hole", resulting in unclear exploration of the gas enrichment areas in front of and on both sides during roadway mining, and thus unable to take more targeted prevention and control measures.

[0005] Since there is a lack of relevant methods for multi-source seismic gas advanced detection of gas enrichment areas in coal mines underground in the prior art, there is an urgent need for a new method that can timely, accurately and effectively monitor the signs of gas outbursts, prevent the occurrence of gas overrun accidents in coal mines, and keep the whole process of gas outburst gestation and formation under effective monitoring. Summary of the Invention

[0006] The object of the present invention is to provide a multi-source seismic gas advance detection method for gas-rich areas in underground coal mines, so as to solve the problem that the prior art cannot effectively predict gas overrun accidents in mines.

[0007] In order to solve the above technical problems, the present invention is implemented by adopting the following technical solutions:

[0008] A multi-source seismic gas advance detection method for gas-rich areas in underground coal mines includes the following steps:

[0009] Step 1: Collect geological data of the target area and establish a seismic wave signal acquisition system along with mining in the target area;

[0010] Step 2: During the coal mining process in the working face of the underground coal mine in the target area, use the seismic wave signal acquisition system along with mining to record all seismic wave signals along with mining as original data;

[0011] Step 3: According to the geological data of the target area and the original data, establish an observation system and evenly divide the target area into multiple grids;

[0012] Step 4: Divide the original data into multiple segments at fixed time intervals, select one segment of data, and determine the actual relative delay time of the shear wave signals of any two seismic wave signals according to the cross-correlation function of the discrete signals of these two seismic wave signals;

[0013] Using this method, traverse the discrete signals of any two seismic wave signals in each segment of data, and obtain the shear wave velocity corresponding to each grid point in the target area according to all actual relative delay times;

[0014] Step 5: According to the shear wave velocity corresponding to each grid point in the target area obtained in Step 4, use the CT imaging algorithm to establish a shear wave velocity model of the working face of the underground coal mine in the target area;

[0015] Step 6: Using the method of Steps 4 - 5, obtain the longitudinal wave velocity corresponding to each grid point in the target area and establish a longitudinal wave velocity model of the working face of the underground coal mine in the target area;

[0016] Step 7: Invert the longitudinal wave velocity model in Step 6, and determine the structural and stress anomaly areas in the target area according to the inversion results;

[0017] Step 8: Calculate the Poisson's ratio corresponding to each grid point according to the shear wave velocity and longitudinal wave velocity corresponding to each grid point in the target area;

[0018] Step 9: Interpret the Poisson's ratio anomaly area according to the Poisson's ratio to complete the detection of the gas-rich area in the coal mining face.

[0019] The present invention also has the following characteristics:

[0020] Further, when establishing the seismic wave signal acquisition system during mining in the target area in Step 1, specifically, a plurality of measuring points are arranged on one side close to the working face of the haulage and return air headings of the coal mine underground working face in the target area to form a seismic wave signal acquisition system during mining;

[0021] The distance between adjacent measuring points is 15 m, and the distance from the first measuring point to the working face does not exceed 20 m.

[0022] Further, Step 4 includes the following sub-steps:

[0023] Step 41, divide the original data into multiple segments at fixed time intervals, select a segment of data, and the cross-correlation function of the discrete signals of any two seismic wave signals in this segment of data is as follows. According to the maximum value of the cross-correlation function, determine the initial travel times of the discrete signals of the two selected seismic wave signals in this selection:

[0024]

[0025] where φ ij represents the cross-correlation function of the discrete signals of the i-th and j-th seismic wave signals;

[0026] x i and x j respectively represent the discrete signals of the i-th and j-th seismic wave signals;

[0027] t i and t j respectively represent the initial travel times of the discrete signals of the i-th and j-th seismic wave signals;

[0028] τ represents the delay time relative to the initial travel time;

[0029] T represents the length of the correlation window;

[0030] t0 represents the time between the initial travel time and the start point of the correlation window;

[0031] δt represents the sampling rate;

[0032] Step 42, use the following formula to calculate the actual relative delay time dr of the discrete signals of the i-th and j-th seismic wave signals selected in this selection in the current segment of data ij :

[0033] dr ij = t i - t j

[0034] Step 43, traverse and calculate the discrete signals of any two seismic channels in the current segment of data, and obtain the actual relative delay times of the discrete signals of all pairwise combinations of seismic wave signals in the current segment of data;

[0035] Step 44: Take the origin time and wave velocity of the seismic source as the initial values for inversion. For all the actual delay times and the theoretical delay times corresponding to each of the actual delay times one by one, take the derivative with respect to the shear wave velocity v, expand, and represent them in matrix form:

[0036] GΔv = d

[0037] where Δv represents the correction value vector composed of all the shear wave velocities v to be differentiated;

[0038] G represents a large sparse partial derivative matrix;

[0039] d represents the matrix composed of all the relative delay times;

[0040] where it is assumed that the target area includes N grid points, and there are a total of Ng arrival times of shear waves detected by geophones;

[0041] where G, d, and Δv are respectively expressed as follows:

[0042]

[0043] where represents the theoretical arrival time difference between the shear waves recorded by the first geophone and the second geophone;

[0044] represents taking the partial derivative with respect to the shear wave velocity v corresponding to the Nth grid point;

[0045] represents the theoretical arrival time difference between the (Ng - 1)th shear wave and the Ngth shear wave;

[0046] Step 45: Use the least squares method to solve the equations in Step 44 to obtain the shear wave velocity corresponding to each grid point in the target area.

[0047] Furthermore, in Step 7, based on all the P-wave velocities, establish a P-wave velocity model from the stress dynamic imaging results obtained by CT inversion, and interpret the structure and stress anomaly areas in the exploration area according to the P-wave velocity model.

[0048] Furthermore, use the following formula to calculate the Poisson's ratio at each grid in the target area:

[0049]

[0050] where μ i represents the Poisson's ratio corresponding to the ith grid point;

[0051] V p i represents the P-wave velocity corresponding to the ith grid point;

[0052] V s i represents the shear wave velocity corresponding to the i-th grid point.

[0053] Furthermore, in step 9, when interpreting the Poisson ratio anomaly area based on the Poisson ratio, according to the combination of structure and stress anomaly area determined in step 6, explore the gas enrichment area in the coal mining face. The high-value area of the Poisson ratio indicates that the anomaly corresponds to the structure and stress anomaly area, forming a corroborative relationship, and thus completing the prediction of the gas enrichment area.

[0054] Furthermore, the time interval in step 4 is 5 min.

[0055] Compared with the prior art, the present invention has the following technical effects:

[0056] The multi-source seismic gas advance detection method for gas enrichment areas in underground coal mines of the present invention can use the seismic data during mining to monitor the dynamic changes of gas in real time. Through continuous seismic data acquisition and analysis, abnormal changes in gas parameters can be detected in a timely manner, providing an early warning for coal mine safety. This real-time monitoring helps the mine operators take measures in a timely manner to prevent the occurrence of gas disasters.

[0057] In addition, the seismic detection technology during mining is adopted to obtain multiple gas-related parameters simultaneously, such as the imaging of bad geological structures causing gas enrichment, the dynamic changes of coal seam stress, and the Poisson ratio of coal bodies. The comprehensive analysis of these parameters helps to more comprehensively understand the occurrence state and outburst risk of gas, improving the accuracy and reliability of prediction.

[0058] Compared with traditional gas monitoring methods, there is no need to directly contact the gas or coal body, thus reducing the safety risks during the detection process. While realizing the advance detection of gas, this method protects the safety of the detection personnel and avoids possible interference with the occurrence state of gas. It has multiple advantages such as high precision, real-time performance, multi-parameter comprehensive analysis, non-contact detection, and wide applicability, providing a strong technical guarantee for coal mine safety production and being suitable for large-scale industrial use and promotion. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 is a flow chart of the multi-source seismic gas advance detection method for gas enrichment areas in underground coal mines of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0060] It should be noted that all components in the present invention, unless otherwise specified, are all components known in the prior art. For example, the coal wall bolts adopt the commonly known devices in the prior art.

[0061] All methods in the present invention, unless otherwise specified, all adopt methods known in the prior art. For example, CT inversion is a method known in the prior art.

[0062] The following are specific embodiments of the present invention. It should be noted that the present invention is not limited to the following specific embodiments, and all equivalent transformations made on the basis of the technical solutions of this application fall within the protection scope of the present invention.

[0063] As Figure 1 shown, a multi-source seismic gas advanced detection method for gas enrichment areas in coal mines includes the following steps:

[0064] Step 1, collect geological data of the target area and establish a seismic wave signal acquisition system during mining in the target area;

[0065] Step 2, during the coal mining face mining process in the target area, use the seismic wave signal acquisition system during mining to record all seismic wave signals during mining as original data;

[0066] Step 3, according to the geological data of the target area and the original data, establish an observation system and evenly divide the target area into multiple grids;

[0067] Among them, when dividing the target area into multiple grids, it is determined according to the actual situation of the target area combined with conventional technical means in the field. There is no fixed division method, and those skilled in the art are also capable of dividing. For example, for an area of 100m * 100m, it can be divided at intervals of 1m into 1000 grids; when divided at intervals of 5m, it is divided into 200 grids.

[0068] Step 4, divide the original data into multiple segments at fixed time intervals, select a segment of data, and determine the actual relative delay time of the shear wave signals of any two seismic wave signals according to the cross-correlation function of the discrete signals of any two seismic wave signals in this segment of data;

[0069] Using this method, traverse the discrete signals of any two seismic wave signals in each segment of data, and obtain the shear wave velocity corresponding to each grid point in the target area according to all actual relative delay times;

[0070] Step 5, according to the shear wave velocity corresponding to each grid point in the target area obtained in Step 4, use the CT imaging algorithm to establish a shear wave velocity model of the coal mining face in the target area;

[0071] Step 6, use the methods in Steps 4-5 to obtain the longitudinal wave velocity corresponding to each grid point in the target area and establish a longitudinal wave velocity model of the coal mining face in the target area;

[0072] Step 7: Invert the longitudinal wave velocity model in Step 6, and determine the structure and stress anomaly areas in the target area according to the inversion results;

[0073] Step 8: Calculate the Poisson's ratio corresponding to each grid point in the target area according to the shear wave velocity and longitudinal wave velocity corresponding to each grid point in the target area;

[0074] Step 9: Interpret the Poisson's ratio anomaly area according to the Poisson's ratio to complete the detection of the gas enrichment area in the coal mining face.

[0075] As a specific implementation method, in Step 1, when arranging multiple measuring points on one side close to the working face in the haulage and return air headings of the coal mining face in the target area, the distance between adjacent measuring points is set to 15 m, and the distance from the first measuring point to the working face does not exceed 20 m.

[0076] The distance between adjacent measuring points is 15 m, and the distance from the first measuring point to the working face does not exceed 20 m.

[0077] Specifically, Step 4 includes the following sub-steps:

[0078] Step 41: Divide the original data into multiple segments at fixed time intervals, select a segment of data, and the cross-correlation function of the discrete signals of any two seismic wave signals in this segment of data is as follows. According to the maximum value of the cross-correlation function, determine the initial travel times of the discrete signals of the two selected seismic wave signals:

[0079] x i and x j respectively represent the discrete signals of the i-th and j-th seismic wave signals. After performing correlation on these two, a new seismic record can be obtained. The new seismic record is also composed of one sample value after another. Read the maximum value among these sample values. After obtaining the maximum value of the cross-correlation function, the initial travel times of the discrete signals of the i-th and j-th seismic wave signals can be obtained by backtracking. Among them, calculating the initial travel times of the discrete signals of the seismic wave signals according to the maximum value of the cross-correlation function is a well-known method in the art;

[0080]

[0081] Among them, φ ij represents the cross-correlation function of the discrete signals of the i-th and j-th seismic wave signals;

[0082] x i and x j respectively represent the discrete signals of the i-th and j-th seismic wave signals;

[0083] t i and t j respectively represent the initial travel times of the discrete signals of the i-th and j-th seismic wave signals;

[0084] τ represents the delay time relative to the initial travel time;

[0085] T represents the length of the correlation window;

[0086] t0 represents the time between the initial travel time and the start point of the correlation window;

[0087] δt represents the sampling rate;

[0088] Step 42, calculate the actual relative delay time dr of the discrete signals of the i-th and j-th seismic wave signals selected this time in the current segment of data using the following formula ij :

[0089] dr ij = t i - t j

[0090] Step 43, traverse the discrete signals of any two seismic channels in each segment of data to obtain the actual relative delay times of the discrete signals of all pairwise combinations of seismic wave signals in the current segment of data;

[0091] Step 44, use the origin time and wave velocity of the seismic source as the initial inversion values, and take the derivative and expand with respect to the shear wave velocity v for all the actual delay times and the theoretical delay times corresponding one by one to all the actual delay times;

[0092] Specifically, use the origin time and wave velocity of the seismic source as the initial inversion values. Assume that the theoretical delay times detected by the i-th and j-th detectors corresponding to the two selected discrete signals are T i and T j , then the difference in the theoretical arrival times measured by the i-th and j-th detectors is ΔT ij ; here, the seismic channels and detectors correspond one by one, that is, one detector corresponds to one seismic channel;

[0093] This theoretical arrival time difference ΔT ij differs from the relative delay time in Step 42. Take the derivative and expand the theoretical arrival time difference ΔT ij and the relative delay time in Step 42 with respect to the shear wave velocity v as follows

[0094]

[0095] Express the above formula in matrix form:

[0096] GΔv = d

[0097] where Δv represents the correction value vector composed of all shear wave velocities v to be differentiated;

[0098] G represents a large sparse partial derivative matrix;

[0099] d represents a matrix composed of all relative delay times;

[0100] Among them, assuming that the target area includes N grid points, there are a total of Ng arrival times of shear waves detected by detectors;

[0101] Among them, G, d, and Δv are respectively represented as follows:

[0102]

[0103]

[0104] Among them, represents the theoretical arrival time difference between the shear waves recorded by the first detector and the second detector;

[0105] represents the partial derivative with respect to the shear wave velocity v corresponding to the Nth grid point;

[0106] represents the theoretical arrival time difference between the shear wave of the (Ng - 1)th trace and the shear wave of the Ngth trace;

[0107] Step 45, use the least squares method to solve the equations in Step 44 to obtain the shear wave velocity corresponding to each grid point in the target area.

[0108] Among them, all elements of the above matrix represent the parameter values of the data within the selected time period, not all the original data. Steps 41 - 45 also calculate for the selected section of data, and Step 43 traverses all pairwise combinations in this section of data. After processing this section of data with Steps 41 - 45, the shear wave and longitudinal wave velocity values of an area corresponding to this period of data are obtained.

[0109] If you also want to detect the shear wave and longitudinal wave velocity values in other areas, just take the data of the next time period and perform the same calculation.

[0110] Specifically, in Step 7, based on all longitudinal wave velocities, establish a longitudinal wave velocity model through CT inversion of the stress dynamic imaging results, and interpret the structure and stress anomaly areas in the exploration area according to the longitudinal wave velocity model.

[0111] Specifically, in Step 8, use the following formula to calculate the Poisson's ratio at each grid in the target area:

[0112]

[0113] where μ i represents the Poisson's ratio corresponding to the i-th grid point;

[0114] Vp i represents the longitudinal wave velocity corresponding to the i-th grid point;

[0115] V s i represents the shear wave velocity corresponding to the i-th grid point.

[0116] Specifically, in step 9, when interpreting the Poisson ratio anomaly area based on the Poisson ratio, according to the combination structure and stress anomaly area determined in step 6, the gas enrichment area in the coal mining face is explored. The high-value area of the Poisson ratio indicates that the anomaly corresponds to the structure and stress anomaly area, forming a corroborative relationship, and thus completing the prediction of the gas enrichment area.

[0117] Specifically, the time interval in step 4 is 5 min.

Claims

1. A multi-source seismic gas advanced detection method for gas enrichment areas in underground coal mines, characterized in that It includes the following steps: Step 1: Collect geological data of the target area and establish a seismic wave signal acquisition system during mining in the target area; Step 2: During the coal mining process in the working face of the coal mine in the target area, use the seismic wave signal acquisition system during mining to record all seismic wave signals during mining as original data; Step 3: According to the geological data of the target area and the original data, establish an observation system and evenly divide the target area into multiple grids; Step 4: Divide the original data into multiple segments at fixed time intervals, select one segment of data, and determine the actual relative delay time of the shear wave signals of any two seismic wave signals according to the cross-correlation function of the discrete signals of any two seismic wave signals in this segment of data; Using this method, traverse the discrete signals of any two seismic wave signals in each segment of data, and obtain the shear wave velocity corresponding to each grid point in the target area according to all actual relative delay times; Step 5: According to the shear wave velocity corresponding to each grid point in the target area obtained in Step 4, use the CT imaging algorithm to establish a shear wave velocity model of the coal mine working face in the target area; Step 6: Using the method of Steps 4 - 5, obtain the longitudinal wave velocity corresponding to each grid point in the target area and establish a longitudinal wave velocity model of the coal mine working face in the target area; Step 7: Invert the longitudinal wave velocity model in Step 6, and determine the structural and stress anomaly areas in the target area according to the inversion results; Step 8: Calculate the Poisson's ratio corresponding to each grid point according to the shear wave velocity and longitudinal wave velocity corresponding to each grid point in the target area; Step 9: Interpret the Poisson's ratio anomaly area according to the Poisson's ratio to complete the detection of the gas enrichment area in the coal mining face; 2. The multi-source seismic gas advanced detection method for gas enrichment areas in underground coal mines according to claim 1, wherein When establishing the seismic wave signal acquisition system during mining in the target area in Step 1, specifically, a plurality of measuring points are arranged on each side of the transportation and return air return roadways close to the working face in the coal mine working face in the target area to form a seismic wave signal acquisition system during mining; The distance between adjacent measuring points is 15m, and the distance from the first measuring point to the working face does not exceed 20m.

3. The multi-source seismic gas advanced detection method for gas enrichment areas in underground coal mines according to claim 1, characterized in that Step 4 includes the following sub-steps: Step 41: Divide the original data into multiple segments at fixed time intervals, select one segment of data, and the cross-correlation function of the discrete signals of any two seismic wave signals in this segment of data is as follows. According to the maximum value of the cross-correlation function, determine the initial travel time of the discrete signals of the two selected seismic wave signals this time: where φ ij represents the cross-correlation function of the discrete signals of the i-th and j-th seismic wave signals; x i and x j respectively represent the discrete signals of the i-th and j-th seismic wave signals; t i and t j respectively represent the initial travel times of the discrete signals of the ith and jth seismic wave signals; τ represents the delay time relative to the initial travel time; T represents the length of the correlation window; t0 represents the time between the initial travel time and the starting point of the correlation window; δt represents the sampling rate; Step 42, calculate the actual relative delay time dr of the discrete signals of the i-th and j-th seismic wave signals selected this time in the current segment data using the following formula ij :[[]]END]] dr ij = t i -t j Step 43: Traverse and calculate the discrete signals of any two seismic channels in the current segment of data, and obtain the actual relative delay time of the discrete signals of all pairwise combinations of seismic wave signals in the current segment of data; Step 44: Take the earthquake occurrence time and wave velocity of the seismic source as the inversion initial values, and represent the derivative of all actual delay times and the theoretical delay times corresponding to all actual delay times with respect to the shear wave velocity v and expand it in matrix form: GΔv=d where, Δv represents the correction value vector composed of all shear wave velocities v to be differentiated; G represents a large sparse partial derivative matrix; d represents the matrix composed of all relative delay times; Among them, it is assumed that the target area includes N grid points, and there are a total of Ng arrival times of shear waves detected by detectors; Among them, G, d, and Δv are respectively expressed as follows: Among them, represents the theoretical wave arrival time difference recorded by the first geophone and the second geophone; Denotes taking the partial derivative with respect to the shear wave velocity v corresponding to the Nth grid point; Denote the theoretical wave arrival time difference between the (Ng - 1)-th shear wave and the Ng-th shear wave; Step 45: Use the least squares method to solve the equations in Step 44 to obtain the shear wave velocity corresponding to each grid point in the target area.

4. The multi-source seismic gas advance detection method for gas enrichment areas in coal mines as described in claim 3, characterized in that, In Step 7, based on all the P-wave velocities, establish a P-wave velocity model according to the stress dynamic imaging results obtained by CT inversion, and interpret the structure and stress anomaly areas in the exploration area according to the P-wave velocity model.

5. For the multi-source seismic gas advanced detection method in the coal mine underground gas enrichment area as described in claim 3, in Step 8, use the following formula to calculate the Poisson's ratio at each grid in the target area: Among them, μ i represents the Poisson's ratio corresponding to the i-th grid point; V p i represents the longitudinal wave velocity corresponding to the i-th grid point; V s i represents the shear wave velocity corresponding to the i-th grid point.

6. The multi-source seismic gas advance detection method for gas enrichment areas in coal mines as described in claim 3, characterized in that In Step 9, when interpreting the Poisson's ratio anomaly area according to the Poisson's ratio, according to the combined structure and stress anomaly area determined in Step 6, explore the gas enrichment area in the coal mining face. The high-value area of the Poisson's ratio indicates that the anomaly corresponds to the structure and stress anomaly area, forming a corroborative relationship, and thus completing the prediction of the gas enrichment area.

7. The multi-source seismic gas advance detection method for gas enrichment areas in underground coal mines according to claim 1, characterized in that, The time interval in Step 4 is 5 min.

Images