Method for advanced detection of karst cave in tunnel earthquake

Abstract

The invention provides a tunnel earthquake advanced karst cave detection method, which comprises the following steps of: arranging a plurality of seismic sources and a plurality of detectors connected with the seismic sources on one side wall of a tunnel to form an earthquake advanced detection system, acquiring observation data of the seismic sources through the system, then calculating to obtain a forward transmission wave field of the seismic sources and a backward transmission wave field of an accompanying source, and calculating the forward transmission wave field of the seismic sources and the backward transmission wave field of the accompanying source; then, according to the forward propagation wave field and the backward propagation wave field, the gradient of the homogeneous medium model parameters is calculated, the accurate updating direction of the homogeneous medium model parameters can be obtained, finally, after the updating direction is solved through the Newton method, iteration updating is conducted on the homogeneous medium model parameters till the target function reaches a set threshold value, and the target function is obtained. And stopping iterative updating to determine parameters of the homogeneous medium model, and performing full-waveform inversion by adopting the homogeneous medium model of the parameters. According to the method, full waveform inversion is carried out by using the full wave field, the karst cave is subjected to classification imaging, the detection imaging resolution can be improved, and thus data support is provided for subsequent mining.

Description

Technical Field

[0001] The invention belongs to the technical field of tunnels, and in particular relates to a method for advanced detection of karst caves through tunnel seismic advance detection. Background Art

[0002] The tunnel excavation stage is a stage where tunnel safety accidents frequently occur, and the main reason is that the location of the cave ahead is unclear. Due to the characteristics of isolated, randomly distributed, and highly concealed caves, traditional drilling methods are costly and blind. Limited by the detection environment of the tunnel, the tunnel forward detection reflection seismic detection method is more practical than other detection methods, but the characteristic waves of the cave are relatively messy. How to return the characteristic waves of the cave to the specific location of the cave is an urgent problem to be solved.

[0003] Full waveform inversion is a waveform inversion method based on the data space domain. It uses the data from the forward simulation to match the actual observed data, establishes the objective function by minimizing the wave field error, and then finds the optimal model parameters to achieve the best match between the simulated data and the observed data. The theory of full waveform inversion was first established on the basis of data domain fitting under the generalized least squares constraint proposed by Tarantola (1982). Tarantola (1984, 1988) applied this theory to acoustic wave approximation and gave a complete theoretical framework for full waveform inversion in the time domain. This technology is relatively mature in ground seismic exploration, and the velocity inversion accuracy of full waveform inversion is much higher than that of traditional velocity inversion methods. Traditional full waveform inversion usually uses direct waves and conjugate gradient method to update the velocity model, and the detection imaging resolution is low. Summary of the invention

[0004] The purpose of the present invention is to provide a method for advanced detection of karst caves by tunnel seismic. By using the full wave field to perform full waveform inversion and classifying the karst caves as images, the detection imaging resolution can be improved.

[0005] The present invention is achieved through the following technical solutions:

[0006] A method for tunnel seismic advance detection of karst caves, comprising the following steps:

[0007] S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are arranged between the seismic sources and the tunnel face, and each of the seismic sources is connected to the multiple geophones to form an earthquake advance detection system;

[0008] S2, establishing a two-dimensional coordinate system along the tunnel, wherein the X direction points to the tunnel face, the Y direction is perpendicular to the tunnel sidewall, and the two-dimensional coordinate system is established with the center of the tunnel face on the model boundary as the origin, and the source position and the detector position are included in the two-dimensional coordinate system;

[0009] S3, exciting the seismic sources in sequence, and acquiring observation data through the earthquake advance detection system;

[0010] S4, establishing an initial uniform medium model, and updating the uniform medium model to obtain a final uniform medium model;

[0011] Wherein, the step of updating the uniform medium model to obtain the final uniform medium model includes:

[0012] S41. Calculate synthetic data and forward wave field based on the uniform medium model and observation data;

[0013] S42, construct a convolution-based objective function F(v p ,v s ,ρ,R r ), which is defined as follows:

[0014]

[0015] In the formula, d is the observed data, q is the synthetic data, and R r is the location parameter of the detector, * is the time convolution operator, R ref Indicates that the position parameters of the reference track are to be extracted. represents the two-norm;

[0016] S43, acquiring the reverse propagation field of the companion source based on the observation data;

[0017] S44, calculating the gradient of the uniform medium model parameter using the forward wave field of the earthquake source and the reverse wave field of the companion source;

[0018] S45, using Newton's method to obtain the update direction, iteratively updating the uniform medium model parameters according to the obtained gradient, until the objective function is less than or equal to the set threshold, stopping the iteration to obtain the final uniform medium model.

[0019] Furthermore, in the step of calculating the synthetic data and the forward wave field based on the uniform medium model and the observed data, the formula for calculating the forward wave field 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 elastic wave forward modeling operator, q(R r,t,v p ,v s ,ρ) is the forward wave field, D(S s ,t) is the earthquake source, S s are the earthquake source coordinates.

[0022] Furthermore, in the step of obtaining the back propagation wave field of the companion 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 ,ρ) is the back propagation wave field, D ′ (S s ,t) is the companion source, where the companion source D ′ (S s ,t) is:

[0025]

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

[0027] Furthermore, in the step of calculating the gradient of the uniform medium model parameter using the forward wave field of the seismic source and the reverse wave field of the companion source, the specific calculation formula is as follows:

[0028]

[0029] In the formula, is the longitudinal wave parameter gradient, is the shear wave parameter gradient, represents the forward wave field q(R r ,t,v p ,v s ,ρ) stress wave field variables, represents the back propagation wave field g(R r ,t,v p ,v s ,ρ) of stress wave variables.

[0030] Further, the Newton method is used to obtain the update direction, and the steps of iteratively updating the uniform medium model parameters according to the obtained gradient include:

[0031] Find the search direction p k :

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

[0033] Where λ is the damping coefficient, I is the unit matrix, and g k is the gradient operator of the objective function, H k is the Hessian matrix,

[0034] Find the update step size α k :

[0035]

[0036] Where T represents the transposed matrix, C k is the Fréchet derivative;

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

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

[0039] In the formula, m k+1 is the updated iterative medium model parameter, m k The medium model parameters of the current iteration.

[0040] Compared with the prior art, the beneficial effects of the present invention are as follows: firstly, the observation data of the earthquake source is obtained, and then the forward wave field of the earthquake source and the reverse wave field of the accompanying source are obtained by calculation, and then the gradient of the uniform medium model parameters is calculated according to the forward wave field and the reverse wave field, which is conducive to obtaining the accurate update direction of the uniform medium model parameters, and finally the velocity field model parameters are iteratively updated until the objective function reaches the set threshold, and the iterative update is stopped to determine the uniform medium model parameters. At this time, the uniform medium model of the parameter is used for full waveform inversion. Therefore, compared with the acoustic wave equation, the present invention uses elastic wave information for inversion, uses more effective signals, and has higher inversion accuracy, thereby improving the detection imaging resolution; and the tunnel advance detection model of the present invention is relatively simple as a whole, and the model is small, which provides an effective way for the application of Newton's method. Newton's method is second-order convergent, and the convergence speed is faster than the conjugate gradient method. The algorithm is also a local convergence algorithm, which improves the problem of slow convergence solved by the traditional full waveform inversion method. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 A schematic diagram of an earthquake advance detection system in a tunnel earthquake advance detection method for karst caves according to the present invention;

[0042] Figure 2 This is a full waveform inversion result diagram of the tunnel seismic advance detection method of the present invention.

[0043] In the figure, 1-seismic source, 2-tunnel, 3-detector, 4-trail face, 5-first cave, 6-second cave. DETAILED DESCRIPTION

[0044] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations.

[0045] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0046] It should be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are only used to distinguish the description and cannot be understood as indicating or implying relative importance.

[0047] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "comprise a ..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0048] In the description of the present invention, it should be noted that the terms "upper", "lower", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings, or are directions or positional relationships in which the product of the invention is usually placed when in use. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore should not be understood as a limitation on the present invention.

[0049] See also Figure 1 and Figure 2 , Figure 1 Schematic diagram of the earthquake advance detection system in the tunnel earthquake advance detection method of the present invention, Figure 2 The full waveform inversion result diagram of the tunnel seismic advance detection method of the present invention is shown in FIG. A tunnel seismic advance detection method of a karst cave comprises the following steps:

[0050] S1. Multiple seismic sources are arranged at equal intervals in the middle of any side wall of the tunnel. Multiple geophones are arranged at equal intervals in the middle of the tunnel side wall between the seismic source and the working face. The geophones can be dual-component geophones. The track spacing between the geophones is L meters. The distance between the seismic source and the nearest geophone is A meters. The interval between the seismic sources is b meters. Figure 1 As shown. Among them, L and A, b and other specific data can be determined according to the actual situation. Then use the seismic equipment to connect each source and detector to form an earthquake advance detection system.

[0051] S2. Establish a two-dimensional coordinate system along the tunnel, where the X direction points to the tunnel face, the Y direction is perpendicular to the tunnel side wall and perpendicular to the X direction, and the two-dimensional coordinate system is established with the center of the tunnel face on the model boundary as the origin, and the source position and the detector position are included in the two-dimensional coordinate system; based on the established two-dimensional coordinate system, the source coordinates and the coordinates of each detector can be obtained;

[0052] S3, sequentially stimulate the seismic source, and obtain observation data through the earthquake advance detection system; specifically, set the sampling interval of the detector to 1 ms, the sampling time length to ts, the sampling length to c, set the explosive quantity to Q, stimulate the seismic source, obtain all the two-component seismic records, thereby obtaining the actual observation data and the source wavelet parameters, and then synthesize the observation data D (S s ,t), and set the main frequency to f Hz;

[0053] S4. Establish an initial uniform medium model. The number of grid points in the X and Y directions of the model is nx and ny respectively, and the grid spacing is dx and dy respectively. The model parameters include the longitudinal wave data v p , shear 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 uniform medium model parameters as the initial model to calculate and obtain updated parameters, thereby obtaining an updated uniform medium model;

[0054] The step of updating the uniform medium model to obtain the final uniform medium model includes:

[0055] S41. Based on the homogeneous medium model, the synthetic data q(R ref ,t,v p ,v s ,ρ) and the forward wave field q(R r ,t,v p ,v s ,ρ), where R r is the location parameter of the detector, v p 、v s and ρ are the parameters of the homogeneous medium model, v p is the longitudinal wave velocity of the medium, v s is the shear wave velocity of the medium, ρ is the density of the medium; the formula for calculating the forward 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 elastic wave forward modeling operator, q(R r ,t,v p ,v s ,ρ) is the forward wave field, D(S s ,t) is the earthquake source, i.e., the observed data, S s is the earthquake source coordinate;

[0058] The forward wave field data q(R r ,t,v p ,v s ,ρ) can be extracted to obtain the synthetic two-component seismic record q(R ref ,t,v p ,v s ,ρ), that is, the synthetic observation data q(R ref ,t,v p ,v s ,ρ); S42, construct the convolution-based objective function F(v p ,v s ,ρ,R r ), which is defined as follows:

[0059]

[0060] In the formula, d is the observed data, q is the synthetic data, and R r is the location parameter of the detector, * is the time convolution operator, R ref Indicates that the position parameters of the reference track are to be extracted. represents the two-norm;

[0061] S43. Obtain the back propagation wave field of the companion source based on the 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 ,ρ) is the back propagation wave field, D ′ (S s ,t) is the companion source, where the companion source D′ (S s ,t) is:

[0064]

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

[0066] S44. The gradient of the uniform medium model parameters is calculated using the forward wave field of the earthquake source and the reverse wave field of the companion source. The specific calculation formula is as follows:

[0067]

[0068] In the formula, is the longitudinal wave parameter gradient, is the shear wave parameter gradient, represents the forward wave field q(R r ,t,v p ,v s ,ρ) stress wave field variables, represents the back propagation wave field g(R r ,t,v p ,v s ,ρ) of stress wave variables.

[0069] S45, using Newton's method to obtain the update direction, iteratively updating the uniform medium model parameters according to the obtained gradient, until the objective function is less than or equal to the set threshold, stopping the iterative update to determine the uniform medium model parameters, and obtaining the final uniform medium model.

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

[0071] Find the search direction p k :

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

[0073] Where λ is the damping coefficient, I is the unit matrix, and g k is the gradient operator of the objective function, H k is the Hessian matrix,

[0074] Find the update step size α k :

[0075]

[0076] Where T represents the transposed matrix, C k is the Fréchet derivative;

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

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

[0079] In the formula, m k+1 is the updated iterative medium model parameter, m k The medium model parameters of the current iteration.

[0080] After the initial uniform medium model is updated by Newton's method, a new uniform medium model is obtained. After the initial uniform medium model is replaced by the new uniform medium model, the objective function F(v p ,v s ,ρ,R r ). If the objective function F(v p ,v s ,ρ,R r ) is less than or equal to the preset threshold, the convergence criterion is met, and the iterative update of the uniform medium model parameters is stopped. The preset threshold can be set to the minimum value F k . Thus, the uniform medium model parameters are determined and the final uniform medium model is obtained. If the objective function F(v p ,v s ,ρ,R r ) is greater than the preset threshold, then the updated uniform medium model is used to repeat steps S41 to S45, and the iteration is repeated until the objective function F(v p ,v s ,ρ,R r ) is less than or equal to F k, stop iterating and updating the uniform medium model, and obtain the final uniform medium model. At this time, the final uniform medium model can be used for full waveform inversion. The present invention obtains the velocity parameters of the cave in front of the entire tunnel by using the full wave field for full waveform inversion, and provides accurate model parameters for later migration imaging, thereby providing data support for subsequent mining. Figure 2 As shown in the inversion result diagram of the present invention, the location of the cave is Figure 1 The locations of the caves are basically consistent, which verifies the effectiveness and imaging accuracy of the present invention.

[0081] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Therefore, any simple modification, equivalent change and modification made to the above embodiment according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A method for advanced detection of karst caves by tunnel seismic, characterized in that: The following steps are involved: S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are arranged between the seismic sources and the tunnel face, and each of the seismic sources is connected to the multiple geophones to form an earthquake advance detection system; S2, establishing a two-dimensional coordinate system along the tunnel, wherein the X direction points to the tunnel face, the Y direction is perpendicular to the tunnel sidewall, and the two-dimensional coordinate system is established with the center of the tunnel face on the model boundary as the origin, and the source position and the detector position are included in the two-dimensional coordinate system; S3, exciting the seismic sources in sequence, and acquiring observation data through the earthquake advance detection system; S4, establishing an initial uniform medium model, and updating the uniform medium model to obtain a final uniform medium model; Wherein, the step of updating the uniform medium model to obtain the final uniform medium model includes: S41, calculating synthetic data and forward wavefield based on the uniform medium model and the observed data; S42, construct a convolution-based objective function F(v p , v s ,ρ,R r ), which is defined as follows: In the formula, d is the observed data, q is the synthetic data, and R r is the location parameter of the detector, * is the time convolution operator, R ref Indicates that the position parameters of the reference track are to be extracted. represents the two-norm; S43, acquiring the reverse propagation field of the companion source based on the observation data; S44, calculating the gradient of the uniform medium model parameter using the forward wave field of the earthquake source and the reverse wave field of the companion source; S45, using Newton's method to obtain the update direction, iteratively updating the uniform medium model parameters according to the obtained gradient, until the objective function is less than or equal to the set threshold, stopping the iteration to obtain the final uniform medium model.

2. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: In the step of calculating the synthetic data and the forward wave field based on the uniform medium model and the observed data, the formula for calculating the forward wave field 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 elastic wave forward modeling operator, q(R r ,t,v p , v s , ρ) is the forward wave field, D(S s , t) is the earthquake source, S s are the earthquake source coordinates.

3. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: In the step of obtaining the back propagation wave field of the companion source based on the observation data, the specific calculation formula 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 , ρ) is the back propagation wave field, D′(S s , t) is the companion source, where the companion source D′(S s , t) is: In the formula, represents the cross-correlation calculation, d(R r , t) is the observed data, d(R ref , t) is the observation data of the reference channel, q(R r ,t,v p , v s , ρ) is the forward wave field, q(R ref ,t,v p , v s , ρ) is the forward wave field at the reference track position.

4. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: In the step of calculating the gradient of the uniform medium model parameter using the forward wave field of the seismic source and the reverse wave field of the companion source, the specific calculation formula is as follows: In the formula, is the longitudinal wave parameter gradient, is the shear wave parameter gradient, represents the forward wave field q(R r ,t,v p , v s ,ρ) stress wave field variables, represents the back propagation wave field g(R r ,t,v p , v s , ρ) of stress wave variables.

5. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: The step of using Newton's method to obtain the update direction and iteratively updating the initial model according to the obtained gradient includes: Find the search direction p k : p k =-(H k +λI) -1 g k (6); Where λ is the damping coefficient, I is the unit matrix, and g k is the gradient operator of the objective function, H k is the Hessian matrix, Find the update step size α k : Where T represents the transposed matrix, C k is the Fréchet derivative; The initial medium model is updated using formula (8), which is specifically: m k+1 =m k +α k p k (8); In the formula, m k+1 is the updated iterative medium model parameter, m k The medium model parameters of the current iteration.

Images