Tunnel seismic cave advanced detection method

Abstract

The application provides a tunnel seismic advanced detection cave method, comprising the following steps: a plurality of seismic sources and a plurality of geophones connected with the seismic sources are arranged on a side wall of a tunnel to form a seismic advanced detection system, observation data of the seismic sources are acquired through the system, then a direct wave field of the seismic sources and a reverse wave field of a companion source are obtained through calculation, then a gradient of a uniform medium model parameter is calculated according to the direct wave field and the reverse wave field, which is beneficial to obtaining an accurate update direction of the uniform medium model parameter, finally, the Newton method is used to obtain the update direction, and then the uniform medium model parameter is iteratively updated until a target function reaches a set threshold, and the iteration and update are stopped to determine the uniform medium model parameter, at this time, full waveform inversion is carried out on the uniform medium model with the parameter. The full waveform inversion is carried out on the full wave field, the cave is imaged, the detection imaging resolution can be improved, and data support is provided for subsequent mining.

Description

Technical Field

[0001] This invention belongs to the field of tunnel technology, and in particular relates to a method for early detection of karst caves in tunnels during earthquakes. Background Technology

[0002] The tunnel excavation stage is a period prone to tunnel safety accidents, primarily due to the unclear location of karst caves ahead. Because karst caves are isolated, randomly distributed, and highly concealed, traditional drilling methods are costly and prone to inaccuracies. Limited by the tunnel's exploration environment, forward-facing reflection seismic exploration is more practical than other methods. However, the characteristic waves of karst caves are relatively chaotic, and how to pinpoint their exact location is a pressing issue.

[0003] Full waveform inversion is a waveform inversion method based on the data spatial domain. It matches forward modeling data with actual observation data, establishes an objective function by minimizing wavefield error, and then finds the optimal model parameters to achieve the best fit between the simulated and observed data. The theory of full waveform inversion was first established based on data domain fitting under the generalized least squares constraint proposed by Tarantola (1982). Tarantola (1984, 1988) applied this theory to the acoustic approximation, thus providing a complete theoretical framework for time-domain full waveform inversion. This technology is relatively mature in ground seismic exploration, and full waveform inversion has a much higher accuracy for velocity inversion than traditional velocity inversion methods. Traditional full waveform inversion usually uses direct waves and the conjugate gradient method to update the velocity model, resulting in relatively low detection imaging resolution. Summary of the Invention

[0004] The purpose of this invention is to provide a method for early detection of karst caves in tunnels by seismic advance detection. By using the full wave field for full waveform inversion, the karst caves can be classified into images, thereby improving the detection and imaging resolution.

[0005] This invention is achieved through the following technical solution:

[0006] A method for early seismic detection of karst caves in tunnels includes the following steps:

[0007] S1. Multiple seismic sources and multiple geophones are equally spaced on one side wall of the tunnel. The multiple geophones are arranged between the seismic sources and the tunnel face. Each seismic source and multiple geophones are connected to form a seismic advance detection system.

[0008] S2. Establish a two-dimensional coordinate system along the tunnel, where the X direction points to the tunnel face and the Y direction is perpendicular to the tunnel sidewall. Establish a two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and incorporate the seismic source position and the detector position into the two-dimensional coordinate system.

[0009] S3. Sequentially stimulate the seismic sources and acquire observation data through the earthquake advanced detection system;

[0010] S4. Establish an initial homogeneous medium model and update the homogeneous medium model to obtain the final homogeneous medium model.

[0011] The step of updating the homogeneous medium model to obtain the final homogeneous medium model includes:

[0012] S41. Based on the homogeneous medium model and observation data, calculate the synthetic data and the propagating wave field;

[0013] S42. Construct a convolution-based objective function F(v) p ,v s ,ρ,R r The specific definition is as follows:

[0014]

[0015] In the formula, d represents the observed data, q represents the composite data, and R0 represents the composite data. r These are the detector's position parameters, * represents the time convolution operator, and R... ref This indicates that the position parameters of the reference track need to be extracted. Represents the L2 norm;

[0016] S43. Obtain the backpropagating wave field of the accompanying source based on the observation data;

[0017] S44. Calculate the gradient of the parameters of the homogeneous medium model using the forward propagation wave field of the source and the reverse propagation wave field of the accompanying source.

[0018] S45. Use Newton's method to find the update direction, and iteratively update the parameters of the homogeneous medium model according to the obtained gradient until the objective function is less than or equal to the set threshold, then stop the iteration to obtain the final homogeneous medium model.

[0019] Furthermore, in the step of calculating the synthesized data and the propagating wavefield based on the homogeneous medium model and observation data, the formula for calculating the propagating wavefield is as follows:

[0020] Q(v p ,v s ,ρ)q(R r ,t,v p ,v s ,ρ)=D(S s ,t) (2);

[0021] In the formula, Q(v) p ,v s ,ρ) is the forward modeling operator for elastic waves, q(R) r,t,v p ,v s ,ρ) represents the propagating wave field, D(S) s (t) represents the earthquake source, S s These are the coordinates of the earthquake source.

[0022] Furthermore, in the step of obtaining the backpropagating wavefield of the accompanying source based on the observation data, the specific calculation formula is as follows:

[0023] Q ′ (v p ,v s ,ρ)g(R r ,t,v p ,v s ,ρ)=D ′ (S s ,t) (3);

[0024] In the formula, ′ is the adjoint operator, Q ′ (v p ,v s ,ρ) is the elastic wave adjoint forward modeling operator, g(R) r ,t,v p ,v s ,ρ) represents the anti-propagating wave field, D ′ (S s Let ,t) be the adjoint source, where the adjoint source D ′ (S s ,t) is:

[0025]

[0026] In the formula, Represents the cross-correlation calculation, d(R) r ,t) represents the observed data, d(R) ref ,t) represents the observation data of the reference channel, q(R) r ,t,v p ,v s ,ρ) represents the propagating wave field, q(R) ref ,t,v p ,v s ,ρ) represents the propagating wave field at the reference trace position.

[0027] Furthermore, in the step of calculating the gradient of the homogeneous medium model parameters using the forward propagating wave field of the earthquake source and the reverse propagating wave field of the accompanying source, the specific calculation formula is as follows:

[0028]

[0029] In the formula, For the gradient of the P-wave parameters, For the gradient of the shear wave parameters, Represents the propagating wave field q(R) r ,t,v p ,v s The stress wave field variables of ρ). Represents the anti-propagating wave field g(R) r ,t,v p ,v s The stress wave variable of ρ).

[0030] Furthermore, the steps of using Newton's method to determine the update direction and iteratively updating the model parameters of the homogeneous medium based on the obtained gradient include:

[0031] Find the search direction p k :

[0032] p k =-(H k +λI) -1 g k (6);

[0033] In the formula, λ is the damping coefficient, I is the identity matrix, and g k It is the gradient operator of the objective function. H k It is the Hessian matrix.

[0034] Find the update step size α k :

[0035]

[0036] In the formula, T represents the transpose matrix, and C k For Fréchet derivative;

[0037] The initial medium model is updated using formula (8), which is as follows:

[0038] m k+1 =m k +α k p k (8);

[0039] In the formula, m k+1 For the updated iterative medium model parameters, m k These are the parameters of the medium model in the current iteration.

[0040] Compared with existing technologies, the advantages of this invention are as follows: First, the observation data of the seismic source is acquired. Then, the forward propagation wave field of the seismic source and the reverse propagation wave field of the accompanying source are calculated. Next, the gradient of the homogeneous medium model parameters is calculated based on the forward and reverse propagation wave fields, which is beneficial to obtaining the accurate update direction of the homogeneous medium model parameters. Finally, the velocity field model parameters are iteratively updated until the objective function reaches a set threshold, at which point the iterative update stops and the homogeneous medium model parameters are determined. At this point, the homogeneous medium model with these parameters is used for full waveform inversion. Therefore, compared with the acoustic wave equation, this invention uses elastic wave information for inversion, utilizes more effective signals, and has higher inversion accuracy, thereby improving the detection imaging resolution. Furthermore, the tunnel advanced detection model of this invention is relatively simple and small, providing an effective way to apply Newton's method. Newton's method is a second-order convergence method, which converges faster than the conjugate gradient method. This algorithm is also a locally convergent algorithm, which improves the problem of slow convergence in traditional full waveform inversion methods. Attached Figure Description

[0041] Figure 1 This is a schematic diagram of the earthquake advance detection system in the tunnel earthquake advance detection method for karst caves of the present invention;

[0042] Figure 2 The image shows the full waveform inversion results of the tunnel seismic advance detection method for karst caves according to the present invention.

[0043] In the diagram, 1-seismic source, 2-tunnel, 3-geometer, 4-working face, 5-first karst cave, 6-second karst cave. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0045] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0046] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0047] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0048] In the description of this invention, it should be noted that the terms "upper," "lower," "inner," "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship in which the product of this invention is usually placed when in use. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting this invention.

[0049] Please see Figure 1 and Figure 2 , Figure 1 This is a schematic diagram of the earthquake advance detection system in the tunnel earthquake advance detection method for karst caves of the present invention. Figure 2 This image shows the full waveform inversion result of the tunnel seismic advance detection method for karst caves according to the present invention. The method includes the following steps:

[0050] S1. Multiple seismic sources are arranged at equal intervals at the midpoint of any sidewall of the tunnel. Multiple geophones are arranged at equal intervals at the midpoint of the tunnel sidewall between the seismic sources and the tunnel face. The geophones can be dual-component geophones. The channel spacing between the geophones is L meters, the distance between the seismic source and its nearest geophone is A meters, and the interval between the seismic sources is b meters. Figure 1 As shown. The specific data for L, A, b, etc., can be determined based on the actual situation. Then, seismic equipment is used to connect each seismic source and detector to form a seismic advance detection system.

[0051] S2. Establish a two-dimensional coordinate system along the tunnel, where the X direction points towards the tunnel face and the Y direction is perpendicular to the tunnel sidewall and perpendicular to the X direction. Establish the two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and incorporate the location of the seismic source and the location of the geophones into the two-dimensional coordinate system. Based on the established two-dimensional coordinate system, the coordinates of the seismic source and the coordinates of each geophone can be obtained.

[0052] S3. Sequentially excite the seismic source and acquire observational data through the earthquake early detection system; specifically, set the geophone sampling interval to l ms, the sampling time length to ts, the sampling length to c samples, and the explosive charge to Q, excite the seismic source, acquire all two-component seismic records, thereby obtaining actual observational data and source wavelet parameters, and then synthesize the observational data D(S) using the actual observational data and source wavelet parameters. s ,t), and set the main frequency to f Hz;

[0053] S4. Establish an initial homogeneous medium model. The number of grid points in the X and Y directions are nx and ny, respectively, and the grid spacing is dx and dy, respectively. The model parameters include the P-wave data v. p transverse wave velocity v s And density ρ, the above parameters can be set according to actual needs. In this embodiment, the longitudinal wave velocity v p =2000m / s, shear wave velocity v s =1000m / s, density ρ=1.0kg / cm³ 3 Using the above-mentioned homogeneous medium model parameters as the initial model, calculations are performed to obtain the updated parameters, thereby obtaining the updated homogeneous medium model.

[0054] The steps for updating the homogeneous medium model to obtain the final homogeneous medium model include:

[0055] S41. Based on the homogeneous medium model, calculate the synthesized data q(R) using the finite difference algorithm. ref ,t,v p ,v s ,ρ) and propagating wave field q(R) r ,t,v p ,v s ,ρ), where R r These are the position parameters of the detector, v p v s ρ and v are parameters of the homogeneous medium model. p It is the longitudinal wave velocity of the medium, v s ρ is the transverse wave velocity of the medium, and ρ is the density of the medium; the formula for calculating the forward propagation wave field is as follows:

[0056] Q(v p ,v s ,ρ)q(R r ,t,v p ,v s ,ρ)=D(S s ,t) (2);

[0057] In the formula, Q(v) p ,v s,ρ) is the forward modeling operator for elastic waves, q(R) r ,t,v p ,v s ,ρ) represents the propagating wave field, D(S) s (t) represents the earthquake source, i.e., the observed data, S s These are the coordinates of the earthquake source;

[0058] The forward propagation wavefield data q(R) r ,t,v p ,v s By extracting all the time values ​​of the detector positions in (ρ), the synthesized two-component seismic record q(R) can be obtained. ref ,t,v p ,v s ,ρ), that is, to obtain the synthetic observation data q(R) ref ,t,v p ,v s S42, Construct a convolution-based objective function F(v,ρ); p ,v s ,ρ,R r The specific definition is as follows:

[0059]

[0060] In the formula, d represents the observed data, q represents the composite data, and R0 represents the composite data. r These are the detector's position parameters, * represents the time convolution operator, and R... ref This indicates that the position parameters of the reference track need to be extracted. Represents the L2 norm;

[0061] S43. Obtain the backpropagating wavefield of the accompanying source based on observation data. The specific calculation formula is as follows:

[0062] Q ′ (v p ,v s ,ρ)g(R r ,t,v p ,v s ,ρ)=D ′ (S s ,t) (3);

[0063] In the formula, ′ is the adjoint operator, Q ′ (v p ,v s ,ρ) is the elastic wave adjoint forward modeling operator, g(R) r ,t,v p ,v s ,ρ) represents the anti-propagating wave field, D ′ (S s Let ,t) be the adjoint source, where the adjoint source D′ (S s ,t) is:

[0064]

[0065] In the formula, Represents the cross-correlation calculation, d(R) r ,t) represents the observed data, d(R) ref ,t) represents the observation data of the reference channel, q(R) r ,t,v p ,v s ,ρ) represents the propagating wave field, q(R) ref ,t,v p ,v s ,ρ) represents the propagating wavefield at the reference trace position;

[0066] S44. Calculate the gradient of the parameters of the homogeneous medium model using the forward propagation wave field of the earthquake source and the reverse propagation wave field of the accompanying source. The specific calculation formula is as follows:

[0067]

[0068] In the formula, For the gradient of the P-wave parameters, For the gradient of the shear wave parameters, Represents the propagating wave field q(R) r ,t,v p ,v s The stress wave field variables of ρ). Represents the anti-propagating wave field g(R) r ,t,v p ,v s The stress wave variable of ρ).

[0069] S45. Use Newton's method to find the update direction, and iteratively update the parameters of the homogeneous medium model according to the obtained gradient until the objective function is less than or equal to the set threshold. Stop the iterative update to determine the parameters of the homogeneous medium model and obtain the final homogeneous medium model.

[0070] Further, in step S45, the step of using Newton's method to determine the update direction and iteratively updating the parameters of the homogeneous medium model based on the obtained gradient includes:

[0071] Find the search direction p k :

[0072] p k =-(H k +λI) -1 g k (6);

[0073] In the formula, λ is the damping coefficient, I is the identity matrix, and g k It is the gradient operator of the objective function. H k It is the Hessian matrix.

[0074] Find the update step size α k :

[0075]

[0076] In the formula, T represents the transpose matrix, and C k For Fréchet derivative;

[0077] The initial medium model is updated using formula (8), which is as follows:

[0078] m k+1 =m k +α k p k (8);

[0079] In the formula, m k+1 For the updated iterative medium model parameters, m k These are the parameters of the medium model in the current iteration.

[0080] After updating the initial homogeneous medium model using Newton's method, a new homogeneous medium model is obtained. Replacing the initial homogeneous medium model with the new model, the objective function F(v) is calculated using the parameters of the new model. p ,v s ,ρ,R r If the objective function F(v) p ,v s ,ρ,R r If the parameter values ​​are less than or equal to a preset threshold, the convergence criterion is met, and iterative updates of the homogeneous medium model parameters cease. The preset threshold can be set to a minimum value F. k This determines the parameters of the homogeneous medium model, thus obtaining the final homogeneous medium model. If the objective function F(v) p ,v s ,ρ,R r If the value is greater than the preset threshold, then the updated uniform medium model is used to repeat steps S41 to S45, and this process is repeated iteratively until the objective function F(v) is satisfied. p ,v s ,ρ,R r ) less than or equal to F kThe iterative update of the homogeneous medium model is stopped, and the final homogeneous medium model is obtained. At this point, the final homogeneous medium model can be used for full waveform inversion. This invention obtains the velocity parameters of the karst cave in front of the entire tunnel by using the full wave field for full waveform inversion, providing accurate model parameters for subsequent migration imaging, thereby providing data support for subsequent mining. Figure 2 As shown in the inversion results diagram of this invention, the location of the karst cave is... Figure 1 The locations of the karst caves are basically consistent, which verifies the effectiveness and imaging accuracy of the invention.

[0081] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Therefore, any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for early seismic detection of karst caves in tunnels, characterized in that, Includes the following steps: S1. Multiple seismic sources and multiple geophones are equally spaced on one side wall of the tunnel. The multiple geophones are arranged between the seismic sources and the tunnel face. Each seismic source and multiple geophones are connected to form a seismic advance detection system. S2. Establish a two-dimensional coordinate system along the tunnel, where the X direction points to the tunnel face and the Y direction is perpendicular to the tunnel sidewall. Establish a two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and incorporate the seismic source position and the detector position into the two-dimensional coordinate system. S3. Sequentially stimulate the seismic sources and acquire observation data through the earthquake advanced detection system; S4. Establish an initial homogeneous medium model and update the homogeneous medium model to obtain the final homogeneous medium model. The step of updating the homogeneous medium model to obtain the final homogeneous medium model includes: S41. Based on the homogeneous medium model and observation data, calculate the synthetic data and the forward propagation wave field; S42. Construct a convolution-based objective function F(v) p ,v s ,ρ,R r The specific definition is as follows: In the formula, d(R) r ,t) represents the observed data, d(R) ref ,t) represents the observation data of the reference channel, q(R) ref ,t,v p ,v s ,ρ) is the synthetic data, which is also the propagating wavefield at the reference trace position, q(R r ,t,v p ,v s ,ρ) represents the propagating wave field, R r These are the position parameters of the detector, v p It is the longitudinal wave velocity of the medium, v s R is the transverse wave velocity of the medium, ρ is the density of the medium, * is the time convolution operator, and R is the transverse wave velocity of the medium. ref This indicates that the position parameters of the reference track need to be extracted. Represents the L2 norm; S43. Obtain the backpropagating wave field of the accompanying source based on the observation data; S44. Calculate the gradient of the parameters of the homogeneous medium model using the forward propagation wave field of the source and the reverse propagation wave field of the accompanying source. S45. Use Newton's method to find the update direction, and iteratively update the parameters of the homogeneous medium model according to the obtained gradient until the objective function is less than or equal to the set threshold, then stop the iteration to obtain the final homogeneous medium model. The steps of using Newton's method to determine the update direction and iteratively updating the initial model based on the obtained gradient include: Find the search direction p k : p k =-(H k +λI) -1 g k (6); In the formula, λ is the damping coefficient, I is the identity matrix, and g k It is the gradient operator of the objective function. H k It is the Hessian matrix. Find the update step size α k : In the formula, T represents the transpose matrix, and C k For Fréchet derivative; The initial medium model is updated using formula (8), which is as follows: m k+1 =m k +α k p k (8); In the formula, m k+1 For the updated iterative medium model parameters, m k These are the parameters of the medium model in the current iteration.

2. The method for advanced seismic detection of karst caves in tunnels according to claim 1, characterized in that, In the step of calculating the composite data and the forward propagation wavefield based on the homogeneous medium model and observation data, the formula for calculating the forward propagation wavefield is as follows: Q(v p ,v s ,ρ)q(R r ,t,v p ,v s ,ρ)=D(S s ,t) (2); In the formula, Q(v) p ,v s ,ρ) is the forward modeling operator for elastic waves, q(R) r ,t,v p ,v s ,ρ) represents the propagating wave field, D(S) s (t) represents the earthquake source, S s These are the coordinates of the earthquake source.

3. The method for advanced seismic detection of karst caves in tunnels according to claim 1, characterized in that, The specific calculation formula for the step of obtaining the backpropagating wavefield of the accompanying source based on the observation data is as follows: Q ′ (v p ,v s ,ρ)g(R r ,t,v p ,v s ,ρ)=D ′ (S s ,t) (3); In the formula, ′ is the adjoint operator, Q ′ (v p ,v s ,ρ) is the elastic wave adjoint forward modeling operator, g(R) r ,t,v p ,v s ,ρ) represents the anti-propagating wave field, D ′ (S s Let ,t) be the adjoint source, where the adjoint source D ′ (S s ,t) is: In the formula, S represents the cross-correlation calculation. s These are the coordinates of the earthquake source.

4. The method for advanced seismic detection of karst caves in tunnels according to claim 1, characterized in that, The specific calculation formula for calculating the gradient of the homogeneous medium model parameters using the forward propagation wavefield of the earthquake source and the reverse propagation wavefield of the accompanying source is as follows: In the formula, For the gradient of the P-wave parameters, For the gradient of the shear wave parameters, Represents the propagating wave field q(R) r ,t,v p ,v s The stress wave field variables of ρ). Represents the anti-propagating wave field g(R) r ,t,v p ,v s The stress wave variable of ρ).

Images