A method for obtaining multi-order refraction times of refracted waves at the coal-rock interface

Abstract

The present invention relates to a method for obtaining multi-order refraction times of refracted waves at the coal-rock interface. When the source wavelength is greater than the coal thickness, the refracted waves of adjacent orders are mutually aliased, and it is difficult for conventional deconvolution methods to sparsify them. By taking advantage of the characteristics of high similarity and known order among the refracted waves of each order and drawing on the orthogonal matching pursuit algorithm, the present invention can accurately extract the refraction times of each order from severely aliased refracted wave signals. The present invention takes into account the phase differences among the refracted waves of each order. By constructing a dictionary matrix containing the phase rotation angle, the refracted waves of each order can accurately match the corresponding wavelets, obtaining accurate pulse signals and improving the accuracy of the obtained time points.

Description

Technical Field

[0001] The present invention relates to a method for obtaining multi-order refraction time, belonging to the technical field of coal seam geological exploration, and specifically relates to a method for obtaining multi-order refraction time of refracted waves at the coal-rock interface. Background Art

[0002] The thickness of the coal seam is an important parameter in the coal mining process. Accurately detecting the distribution of the coal seam thickness in the working face before coal mining is the premise and foundation for the shearer to dynamically adjust its position and cutting height, so as to realize the intelligent mining of the coal mine working face, and is of great significance for improving the coal mining efficiency, reducing resource waste, and ensuring the safe production of the coal mine.

[0003] After the seismic source is excited in the coal seam, the generated seismic waves will be reflected multiple times between the two interfaces of the coal seam and the roof and floor. Only the seismic waves incident at the critical angle can generate refracted waves propagating along the interface. Therefore, when the thickness of the coal seam is stable, the refracted wave part of the received seismic signal will have obvious periodicity, that is, after the first arrival refracted wave appears at the time point corresponding to the velocity of the surrounding rock, the subsequent refracted waves will appear every fixed time period. When there are abnormal coal thickness areas such as thin coal belts in the working face, the time points at which the subsequent refracted waves appear will shift and deviate from the normal period. Therefore, by finding the time when each order of refracted waves appears, the thickness of the abnormal coal seam can be calculated therefrom.

[0004] The patent "A method and system for detecting thin coal belts based on the attenuation of refracted wave period and amplitude" (publication number: 202111434047.5) uses the periodic characteristics of the refracted waves propagating at the interface between the coal seam and the surrounding rock to detect the positions and boundaries of abnormal coal thickness areas such as thin coal belts in the coal mining working face. This patent method can only delineate the range of the thin coal belt and cannot accurately estimate the specific coal thickness and its distribution within the thin coal belt. This is because the method of deconvolving the refracted wave signal in this patent is relatively simple, and the refracted wave signal after deconvolution is not sparse, from which the specific time of multi-order refraction cannot be accurately obtained, so the coal seam thickness corresponding to each refraction point cannot be determined. Therefore, it is necessary to adopt a more effective algorithm to sparsify the refracted wave signal according to the characteristics of the refracted waves and obtain the multi-order refraction time therefrom, so as to create conditions for accurately estimating the coal thickness distribution. Summary of the Invention

[0005] The following presents a brief overview of one or more aspects to provide a basic understanding of these aspects. This overview is not an exhaustive survey of all contemplated aspects, and is neither intended to identify key or decisive elements of all aspects nor to attempt to define the scope of any or all aspects. Its sole purpose is to present some concepts of one or more aspects in a simplified form as a prelude to the more detailed description that follows.

[0006] The main object of the present invention is to solve the technical problems existing in the prior art, and provides a method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface. This method is simple and easy to use, and is suitable for the observation system and construction conditions of the existing trough wave detection project. It can extend the applicable conditions of the coal thickness detection method based on the refraction wave from the original extremely thick coal seam to the thick coal seam and medium-thick coal seam.

[0007] To solve the above problems, the solution of the present invention is as follows:

[0008] A method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, including:

[0009] Calculating the refraction order K of the refraction wave, and synthesizing the wavelet w n (t) into the dictionary matrix W;

[0010] Calculating the sparse signal r(t) according to the following iterative steps:

[0011] Iterative step 1, initialization, let k = 1, A0 be an empty set, and u1 = s;

[0012] Iterative step 2, finding the index i k , such that , that is, taking the inner product of all column vectors of W and u k , and finding the column number corresponding to the column with the largest inner product;

[0013] Iterative step 3, obtaining the iterated A based on the following formula k :

[0014] A k = A k-1 ∪ a ik ;

[0015] Iterative step 4, finding the least squares solution of s = A k r k , that is:

[0016] Iterative step 5, updating the residual vector u k+1 = s - A k r k ;

[0017] Iterative step 6, let k = k + 1, if k ≤ K, then return to iterative step 2, otherwise output the current vector r k as the sparse signal r(t);

[0018] Iterative step 7, picking up the time corresponding to the K pulse peaks of the sparse signal r(t) to obtain the multi-order refraction time of the refraction wave;

[0019] Among them, in each iterative step, k is the current iteration number, and r kDenote the sparse vector of the k-th iteration, u k Denote the residual vector of the k-th iteration, i k Denote the index found in the k-th iteration, where the index is the corresponding column number in the W matrix, a ik Denote the column vector in matrix W with i k as the serial number, A k is a set matrix composed of all selected column vectors a ik .

[0020] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, calculate the refraction order K based on the following formula:

[0021] From the known positions of the excitation point (x s , y s ) and the receiving point (x r , y r ), calculate the refraction order K of the refraction wave:

[0022]

[0023] In the formula, the seismic wave velocity v r , the seismic wave velocity v of the coal seam c , and T is the refraction wave period.

[0024] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, calculate the refraction wave period based on the following formula:

[0025]

[0026] Wherein, the seismic wave velocity v of the surrounding rock r , the seismic wave velocity v of the coal seam c .

[0027] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, generate a wavelet based on the following formula:

[0028]

[0029] Where t ∈ [0, NΔt], the period phase rotation angle Φ, the refraction wave period T, t represents time, t ∈ [0, NΔt], N is the length of the collected refraction wave signal, and Δt is the sampling interval.

[0030] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, generate a dictionary matrix W based on the following formula:

[0031]

[0032] Where the symbol ∪ represents the set union operation.

[0033] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, the following formula is used to obtain the when it reaches the minimum as the periodic phase rotation angle Φ:

[0034]

[0035] wherein,

[0036]

[0037]

[0038] the Ricker wavelet w0(τ) with the main frequency of f c , τ represents time, τ ∈ [-T, T], and T is the refraction wave period;

[0039] H(·) is the Hilbert transform, τ ∈ [-T, T].

[0040] Preferably, for the above method for obtaining the multi-order refraction time of the refraction wave at the coal-rock interface, the Ricker wavelet w0(τ) with the main frequency of f c is constructed based on the following formula:

[0041]

[0042] wherein, τ represents time, τ ∈ [-T, T], and T is the refraction wave period.

[0043] Therefore, compared with the prior art, the advantages of the present invention are:

[0044] (1) The present invention can realize the processing of refraction waves under the condition that the source wavelength is greater than the coal thickness. When the source wavelength is greater than the coal thickness, the refraction waves of adjacent orders are overlapped with each other, and it is difficult for the conventional deconvolution method to sparsify them. The present invention utilizes the characteristics of high similarity and known order of each order of refraction waves, and draws on the orthogonal matching pursuit algorithm, and can accurately extract the refraction time of each order from the severely overlapped refraction wave signals.

[0045] (2) The present invention takes into account the phase difference between each order of refraction waves. By constructing a dictionary matrix containing the phase rotation angle, each order of refraction waves can accurately match the corresponding wavelet, and an accurate pulse signal is obtained, thereby improving the accuracy of the obtained time points. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] The drawings incorporated herein and forming a part of the specification illustrate embodiments of the present invention, and together with the description further serve to explain the principles of the present invention and enable those skilled in the art to make and use the present disclosure.

[0047] Figure 1 It is the general flow chart of the method for obtaining the multi - order refraction time of the refracted wave at the coal - rock interface.

[0048] Figure 2 It is the schematic diagram of the working face model. Among them, (a) is the plan view of the model and the observation system, and (b) is the cross - sectional view of the model.

[0049] Figure 3 It is the 3D elastic wave forward simulation result of the different main - frequency vibration sources at the working face. Among them, (a) is the refracted wave signal generated by the 200Hz main - frequency vibration source, (b) is the refracted wave signal generated by the 400Hz main - frequency vibration source, and (c) is the comparison chart of the refracted wave signals of the two main frequencies received by R15.

[0050] Figure 4 It is the refracted wave signal near the vibration source and its frequency spectrum. Among them, (a) is the refracted wave signal and (b) is the frequency spectrum.

[0051] Figure 5 It is to construct a Ricker wavelet with 200Hz as the main frequency.

[0052] Figure 6 It is the error function between the wavelet and the refracted wave signal

[0053] Figure 7 It is the comparison chart of the refracted wave signal s(t) of the R15 channel and the synthesized signal obtained by performing periodic phase rotation on the wavelet w0 with Φ.

[0054] Figure 8 It is the image corresponding to the dictionary matrix W.

[0055] Figure 9 It is the schematic diagram of s(t) and the obtained sparse signal r(t).

[0056] Figure 10 It is the sparse signal obtained by processing all data channel by channel.

[0057] The embodiments of the present invention will be described with reference to the accompanying drawings. Detailed implementation manners

[0058] Embodiment

[0059] When the wavelength of the vibration source is greater than the coal thickness, the refracted waves of adjacent orders are mutually aliased, and it is difficult to sparsify them by the conventional deconvolution method. Therefore, the method of this embodiment utilizes the characteristics that the refracted waves of each order have high similarity and the order number is known, draws on the orthogonal matching pursuit algorithm, and takes into account the phase difference between the refracted waves of each order. By constructing a dictionary matrix containing the phase rotation angle, the refracted waves of each order can accurately match the corresponding wavelets, obtain accurate pulse signals, and improve the accuracy of the obtained time points.

[0060] The following is a detailed description of this embodiment in conjunction with the accompanying drawings. As Figure 1 shown, the method for obtaining the multi-order refraction time of the refracted wave at the coal-rock interface in this embodiment includes:

[0061] Step 1, the refracted wave period calculation step. From the known seismic wave velocity v of the surrounding rock r , the seismic wave velocity v of the coal seam c , and the average thickness h of the coal seam, calculate the refracted wave period T:

[0062]

[0063] Step 2, from the known positions of the excitation point (x s , y s ) and the receiving point (x r , y r ), calculate the refraction order K of the refracted wave:

[0064]

[0065] Step 3, estimate the source main frequency f c .

[0066] Select the refracted wave signal s0(t) near the excitation point, where t represents time, t ∈ [0, NΔt], N is the length of the collected refracted wave signal, and Δt is the sampling interval. Use Fourier transform to obtain its amplitude spectrum, and find the frequency corresponding to the maximum amplitude, which is the source main frequency f c .

[0067] Step 4, construct the Ricker wavelet w0(τ) with the main frequency f c .

[0068]

[0069] Where τ represents time, τ ∈ [-T, T], and T is the refracted wave period.

[0070] Step 5, estimate the periodic phase rotation angle Φ, that is, the angle by which the phase of the refracted wave signal rotates every time it passes through a period, specifically including:

[0071] Step 5.1, find the first arrival time t0 of the refracted wave:

[0072]

[0073] Step 5.2, take Calculate the error function between the wavelet and the refracted wave signal

[0074]

[0075] Where:

[0076]

[0077] H(·) is the Hilbert transform, and τ ∈ [-T, T].

[0078] Step 5.3, find the when it reaches the minimum Then is the estimated periodic phase rotation angle Φ.

[0079] Step 6, generate an N×N dictionary matrix W, where N is the signal length, including the following steps:

[0080] Step 6.1 Let n = 0 and w be an empty set;

[0081] Step 6.2, generate the wavelet w n (t):

[0082]

[0083] where t ∈ [0, NΔt].

[0084] Step 6.3, denote the wavelet w n (t) as the N-dimensional vector w n ;

[0085] Step 6.4, let n = n + 1, if n < N, then return to Step 6.2;

[0086] Step 6.5, combine all w n into a matrix W:

[0087]

[0088] where the symbol ∪ represents the set union operation.

[0089] Step 7, take the refraction order K as the input parameter, and use the orthogonal matching pursuit algorithm to perform sparse deconvolution on the refraction signal s(t) with the dictionary matrix W to obtain the sparse signal r(t). The specific implementation is as follows:

[0090] Denote the refraction wave signal s(t) of length N as the N-dimensional vector s, let k be the current iteration number, r k represents the sparse vector at the k-th iteration, u k represents the residual vector at the k-th iteration, i k represents the index found at the k-th iteration (the corresponding column number in the W matrix), a ik represents the column vector in the matrix W with i k as the serial number, and A k is a matrix composed of all the selected column vectors a ikThe set matrix formed. The symbol ∪ represents the union operation of sets, and the symbol <·, ·> represents the inner product of vectors. Iterate according to the following steps:

[0091] Step 7.1, Initialization: Let k = 1, A0 be an empty set, and u1 = s;

[0092] Step 7.2, Find the index i k , such that That is, take the inner product of all column vectors of W with u k and find the column number corresponding to the column that gives the maximum inner product;

[0093] Step 7.3, A k = A k-1 ∪ a ik ;

[0094] Step 7.4, Find the least - squares solution of s = A k r k :

[0095] Step 7.5, Update the residual vector u k+1 = s - A k r k ;

[0096] Step 7.6, Let k = k + 1. If k ≤ K, return to step (2); otherwise, output the current vector r k as the sparse signal r(t);

[0097] Step 8, The sparse signal r(t) has K pulses. Pick up the times corresponding to the peaks of these K pulses to obtain the multi - order refraction times of the refracted wave.

[0098] The following takes the forward - modeling data as an example to illustrate the effect of the present invention:

[0099] The model consists of three layers: the roof, the coal seam, and the floor. The lithology of the roof and the floor is the same. The longitudinal wave velocity of the surrounding rock is 4000 m / s, the transverse wave velocity of the surrounding rock is 2300 m / s, and the density is 2.56 g / cm 3 ; The longitudinal wave velocity of the coal seam is 2000 m / s, the transverse wave velocity of the coal seam is 1050 m / s, and the density is 1.4 g / cm 3 . The thickness of the coal seam is 10 m, and there are two parallel roadways in the coal seam. The distance between the two roadways is 200 m. There is a thin - coal area in the coal seam, with a size of 100 m × 50 m and a coal thickness of 5 m. Figure 2 (a) is the plan view of the model and the observation system, Figure 2 (b) is the cross - sectional view of the model.

[0100] The excitation points and the receiving points are respectively arranged in the two roadways, and the positions are as Figure 2(a) As shown. Using Ricker wavelets with different main frequencies as the source for three-dimensional elastic wave numerical simulation, the refracted wave signals generated by the 200 Hz and 400 Hz main frequency sources are respectively as Figure 3 (a) and 3(b) shown. Figure 3 (c) is a comparison graph of the refracted wave signals of two frequencies received by R15. It can be seen that the refracted waves of each order in the 400 Hz signal are separated from each other, and the refracted waves of each order can be intuitively distinguished, and the refraction time corresponding to the maximum amplitude of each order can be found from them; while the 200 Hz refracted waves are aliased with each other, and the waveforms of the refracted waves of each order cannot be intuitively distinguished, so the refraction time of each order cannot be found.

[0101] Now, according to the steps of the present invention, the refracted wave signal s(t) received by R15 is processed to find the refraction time of each order.

[0102] Execute step 1. From the model parameters v r = 4000 m / s, v c = 2000 m / s, h = 10 m, the refraction wave period T = 86.6 ms is obtained.

[0103] Execute step 2. From the position of the excitation point (x s , y s ) = (400, 55) and the position of the receiving point R15 (x r , y r ) = (140, 250), K = 5 is obtained.

[0104] Execute step 3. The refracted wave signal s0(t) near the source and its spectrum are respectively as Figure 4 (a) and 4(b) shown, and the main frequency f c is 200 Hz.

[0105] Execute step 4. Construct a Ricker wavelet w0 with 200 Hz as the main frequency, as Figure 5 shown.

[0106] Execute step 5(1), and calculate the refraction wave arrival time t0 = 71 ms.

[0107] Execute step 5(2). The error function between the wavelet and the refracted wave signal is as Figure 6 shown.

[0108] Execute step 5(3). From Figure 5 the period phase rotation angle Φ = 1.24 is estimated. Figure 7 is a comparison graph of the refracted wave signal s(t) of the R15 trace and the synthesized signal obtained by performing period phase rotation on the wavelet w0 with Φ.

[0109] Execute step 6 to synthesize the dictionary matrix W. Since the length of the intercepted refracted wave signal is 1600 sampling points, the size of the generated W matrix is 1600×1600. Taking each element in W as each point in the image, the W matrix can be output as a 1600×1600 image. As Figure 8 shown, the abscissa is the row of the matrix, the ordinate is the column of the matrix, and the values in the image are the colors corresponding to the elements in the matrix (such as grayscale values).

[0110] Execute step 7 to process the refracted wave signal s(t) as follows:

[0111] Execute step 7(1) to complete the initialization. Let k = 1, set A0 as an empty set, and set the data in the signal s(t) as u1;

[0112] Execute step 7(2). Take the inner product of each column in W and u1, and find that the column with the largest inner product is the 1108th column, then i k = 1108;

[0113] Execute step 7(3) to merge A0 and the 1108th column in W into A1;

[0114] Execute step 7(4) to calculate r1;

[0115] Execute step 7(5) to calculate u1;

[0116] Execute step 7(7). Since k ≤ K, let k = k + 1, and return to step 7(2) until k > K; s(t) and the obtained sparse signal r(t) are as Figure 9 shown.

[0117] Execute step 8. From Figure 9 the sparse signals, extract the times corresponding to the maximum value points respectively, which are (76, 849, 938, 102.4, 1108, 1174) ms.

[0118] It can be seen that although there is aliasing between the refracted waves of each order, the refraction times of each order obtained by using the present invention differ by basically 86 ms, that is, one refraction wave period. Therefore, the present invention can accurately obtain the refraction times of each order. Processing all data channel by channel according to the above steps, the obtained sparse signals are as Figure 10 shown. Comparing Figure 3 with a, it can be clearly seen that there are abnormalities in the refraction times of each order in the circled area on the right. This kind of abnormality is caused by the thin coal area in the model.

[0119] Generally speaking, the present invention uses the periodicity of the refracted wave propagating on the coal-rock interface to perform sparse deconvolution on it, and obtains multi-order refraction times from the sparse signal. Accurate multi-order refraction times are a prerequisite for delineating the abnormal coal thickness area and estimating the coal seam thickness in the abnormal area, and are of great significance for ensuring the intelligent mining of coal mine working faces, improving the recovery efficiency and extraction rate, and reducing resource waste.

[0120] In this embodiment, although the above method is illustrated and described as a series of actions for simplicity of explanation, it should be understood and appreciated that these methods are not limited by the order of the actions, because according to one or more embodiments, some actions may occur in a different order and / or concurrently with other actions that are illustrated and described herein or that are not illustrated and described herein but are understood by those skilled in the art.

[0121] Note that references in the specification to "one embodiment", "an embodiment", "example embodiment", "some embodiments", etc. indicate that the described embodiments may include specific features, structures, or characteristics, but each embodiment may not necessarily include the specific features, structures, or characteristics. Moreover, such phrases do not necessarily refer to the same embodiment. In addition, when a specific feature, structure, or characteristic is described in connection with an embodiment, whether or not explicitly described, implementing such a feature, structure, or characteristic in connection with other embodiments will be within the knowledge of those skilled in the art.

[0122] The previous description of the present disclosure is provided to enable any person skilled in the art to make or use the present disclosure. Various modifications to the present disclosure will be apparent to those skilled in the art, and the general principles defined herein can be applied to other variations without departing from the spirit or scope of the present disclosure. Thus, the present disclosure is not intended to be limited to the examples and designs described herein, but should be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for obtaining multi-order refraction times of refracted waves at the coal-rock interface, characterized in that, Including: Calculate the refraction order K of the refracted wave and synthesize the wavelet w n (t) into the dictionary matrix W; Calculating the sparse signal r(t) according to the following iterative steps: Iterative step 1, initialization: Let k = 1, A0 be an empty set, and u1 = s; Iterative step 2, find index i k , such that compute the inner product of all column vectors of W with u k and find the column number corresponding to the column that gives the maximum inner product; Iterative step 3, calculate the iterated A based on the following formula k :[[-END]] A k = A k-1 ∪ a ik ; Iterative step 4, find s = A k r k The least squares solution of, i.e.: Iterative step 5, update the residual vector u k+1 = s - A k r k ; Iterative step 6, let k = k + 1. If k ≤ K, return to iterative step 2; otherwise, output the current vector r k as the sparse signal r(t); Picking up the times corresponding to the K pulse peaks of the sparse signal r(t) to obtain the multi-order refraction times of the refracted wave; Among them, in each iteration step, k is the current iteration number, and r k represents the sparse vector of the k-th iteration, u k represents the residual vector of the k-th iteration, i k represents the index found in the k-th iteration, and the index is the corresponding column number in the W matrix, a ik represents the column vector in the matrix W with i k as the serial number, A k is a set matrix composed of all the selected column vectors a ik ; Calculating the refraction order K based on the following formula: From the known excitation point positions (x s , y s ) and receiver point positions (x r , y r ), calculate the refraction order K of the refracted wave: In the formula, \(v\) is the seismic wave velocity of the surrounding rock r and \(v\) is the seismic wave velocity of the coal seam c , and \(T\) is the refraction wave period.

2. The multi - order refraction time calculation method for the refracted wave at the coal - rock interface according to claim 1, characterized in that, Calculating the refracted wave period based on the following formula: Among them, the seismic wave velocity v of the surrounding rock r and the seismic wave velocity v of the coal seam c ; the average thickness h of the coal seam.

3. A method for obtaining multi - order refraction times of a refracted wave at a coal - rock interface according to claim 1, characterized in that, Generating a wavelet based on the following formula: where t ∈ [0, NΔt], the periodic phase rotation angle Φ, the refracted wave period T, t represents time, t ∈ [0, NΔt], N is the length of the collected refracted wave signal, and Δt is the sampling interval.

4. A method for obtaining multi-order refraction times of refracted waves at a coal-rock interface according to claim 3, characterized in that, Generating a dictionary matrix W based on the following formula: where the symbol ∪ represents the set union operation.

5. A method for obtaining the multi-order refraction time of a refracted wave at a coal-rock interface according to claim 3, characterized in that Calculate, based on the following formula, the when it reaches the minimum, and use the as the periodic phase rotation angle Φ: where, The main frequency is f c Ricker wavelet , Indicates time, ∈[-T,T], T is the period of the refracted wave; H(·) is the Hilbert transform, ∈[-T, T]; the position of the excitation point (x s , y s ), the position of the receiving point (x r , y r ); is the estimated periodic phase rotation angle Φ.

6. A method for obtaining multi-order refraction time of refracted waves at the coal-rock interface according to claim 5, characterized in that, Construct a Ricker wavelet with a main frequency of f based on the following formula c as follows :[[]] Among them, represents time, ∈[-T, T], where T is the refraction wave period.

Images