Method, system, medium and equipment for advanced detection of karst cave in tunnel earthquake

Abstract

The invention discloses a tunnel earthquake karst cave advanced detection method, system, medium and equipment, and the method comprises the steps: S1, a plurality of seismic sources and a plurality of detectors are arranged on one side wall of a tunnel at equal intervals, the plurality of detectors are located between the seismic sources and a tunnel face, the plurality of detectors are arranged on the tunnel face at equal intervals, and the heights of the detectors and the seismic sources are consistent; s2, respectively connecting each seismic source and all detectors through equipment to form a seismic advanced detection system, and sequentially exciting the seismic sources to obtain actual observation data; and S3, establishing an initial velocity model, and updating the velocity model by using a Newton method to obtain a final velocity model. According to the method, a large-azimuth earthquake advanced observation system and a full-waveform inversion method are combined, the problem that the traditional zero offset transverse resolution is low is solved, meanwhile, the Newton method is used for replacing a conjugate gradient method of the traditional full-waveform inversion, the calculation efficiency of the full-waveform inversion is further improved, and the tunnel advanced detection imaging precision can be improved.

Description

Technical Field

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

[0002] The tunnel excavation stage is a stage where tunnel safety accidents frequently occur. The main reason is that the structural problems ahead are not clearly explored, among which the fault structure is the most typical. Therefore, accurate detection of faults is of great significance. Since the drilling cost is relatively high, and the reflection seismic exploration method for tunnel advance detection is more practical and convenient, the seismic reflection method is the most common in tunnel advance detection. In the seismic reflection advance detection work, how to establish an accurate fault velocity model is the key to realizing seismic advance detection. An accurate fault velocity model provides favorable conditions for realizing advance detection of faults ahead of the tunnel.

[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 observation 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 simulation data and the observation 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. However, the conjugate gradient method is currently used to optimize the model velocity parameters in full waveform inversion. The optimization speed of this method is relatively slow, which leads to the low computational efficiency of full waveform inversion. Summary of the invention

[0004] The first purpose of the present invention is to provide a method for tunnel seismic advance detection of caves, which combines a large-azimuth observation system with a full waveform inversion method to overcome the problem of low lateral resolution of traditional zero offset distance, and at the same time uses Newton's method to replace the conjugate gradient method of traditional full waveform inversion to further improve the computational efficiency of full waveform inversion, which helps to improve the accuracy of tunnel advance detection imaging.

[0005] The second object of the present invention is to provide a tunnel seismic advance detection system for karst caves.

[0006] A third object of the present invention is to provide a storage medium.

[0007] A fourth object of the present invention is to provide a computing device.

[0008] The purpose of the present invention is achieved through the following technical solutions: A method for tunnel seismic advance detection of karst caves, comprising the following steps:

[0009] S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, the multiple geophones are located between the seismic sources and the tunnel face, and arrange multiple geophones at equal intervals on the tunnel face, with each geophone and each seismic source located at the same height;

[0010] S2, connecting each seismic source and all geophones through equipment to form an earthquake advance detection system, exciting the seismic sources in turn, and obtaining actual observation data through the earthquake advance detection system;

[0011] S3, establishing an initial velocity model, and updating the velocity model using Newton's method to obtain a final velocity model;

[0012] The steps of updating the velocity model using Newton's method to obtain the final velocity model include:

[0013] S31, determining synthetic observation data and forward wave field data according to the velocity model;

[0014] S32, obtaining back propagation wave field data according to the velocity model and actual data;

[0015] S33, determining a target equation of the residual according to the forward wavefield data and the reverse wavefield data;

[0016] S34, calculating a gradient matrix related to wave velocity using the forward wave field data and the reverse wave field data;

[0017] S35, calculating the Hessian matrix based on the gradient matrix;

[0018] S36, updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size;

[0019] S37, repeating steps S31 to S36 until a preset convergence criterion is met to obtain a final velocity model.

[0020] Furthermore, the step of determining the synthetic observation data and the forward wavefield data according to the velocity model includes:

[0021] The observation system parameters are determined according to the earthquake advance detection system. Based on the observation system parameters, the forward wave field data are calculated using the space-time high-order finite difference algorithm:

[0022] L[v n (x,y)]p n (x,y,t)=S(x S ,y S ,t) (1);

[0023] Where L[·] is the high-order finite difference operator of acoustic wave space-time, v n (x, y) is the velocity model after the nth iteration, p n (x,y,t) is based on v n (x,y) forward wave field data, S(x S ,y S ,t) is the source wavelet, (x S ,y S ) is the earthquake source location, t is the time parameter;

[0024] The synthetic observation data is determined based on the forward wavefield data.

[0025] Furthermore, the step of obtaining back propagation wavefield data according to the velocity model and the actual observation data includes:

[0026] Taking the actual observed data as the boundary condition, the reverse wave field data is calculated using formula (2):

[0027] L[v n (x,y)]d n (x,y,t)=D(S;x1,y1;x2,y2;…x M+N ,y M+N ;t) (2);

[0028] In the formula, v n (x, y) is the velocity model after the nth iteration, d n (x,y,t) is based on v n (x, y) reverse wave field data, D(S; x1, y1; x2, y2; ...x M+N ,y M+N ; t) is the actual observed data, (x1, y1; x2, y2; ...x M+N ,y M+N ) is the location parameter of the geophone, N+M represents the total number of geophones, t is the time, and S represents the number of sources.

[0029] Further, in the step of determining the target equation of the residual according to the forward wavefield data and the backward wavefield data, the expression of the target equation is:

[0030]

[0031] In the formula, F[v n (x,y)] is the objective function value after the nth iteration, v n (x, y) is the velocity model after the nth iteration, p n (x,y,t) is based on v n (x, y) forward wave field data, d n(x,y,t) is based on v n (x, y) is the reverse wave field data, ‖·‖ is the two-norm, ζ is the regularization scale factor, W is the regularization weight operator, v apr is the expected approximate solution.

[0032] Furthermore, in the step of calculating the gradient matrix about the wave velocity using the forward wave field data and the reverse wave field data, the expression of the gradient matrix is:

[0033]

[0034] In the formula, F[v n (x,y)] about v n The first-order partial derivative of (x,y), v n (x,y) is the velocity model after the nth iteration.

[0035] Furthermore, in the step of calculating the Hessian matrix based on the gradient matrix, the expression of the Hessian matrix is:

[0036]

[0037] In the formula, k n is the gradient update amount, k n =g[v n+1 (x,y)]-g[v n (x,y)],q n For the velocity model v n The update amount of (x,y), q n =v n+1 (x,y)-v n (x,y), represents the matrix transpose, is the initial positive definite Hessian matrix of the nth iteration, For Based on this, the approximate Hessian matrix of the jth recursive iteration is performed, express About v n The first partial derivative of (x,y).

[0038] Further, the step of updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size includes:

[0039] The velocity model is updated using formula (6), which is expressed as:

[0040] v n+1 (x,y)=v n (x,y)+a n p n(6);

[0041] In the formula, a n is the iteration step length, v n (x, y) is the velocity model after the nth iteration, v n+1 (x, y) is the velocity model after the n+1th iteration, p n is the iteration direction, p n The expression is;

[0042]

[0043] Where γ is the damping coefficient and I is the unit matrix.

[0044] A tunnel seismic advance detection system for karst caves, comprising:

[0045] A detection system erection unit, used to arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are located between the seismic sources and the tunnel face, and multiple geophones are arranged at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height;

[0046] An observation data acquisition unit, used to connect each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, stimulate the seismic sources in turn, and acquire actual observation data through the earthquake advance detection system;

[0047] A velocity model updating unit, used for establishing an initial velocity model, updating the velocity model using Newton's method, and obtaining a final velocity model;

[0048] Wherein, the updating of the velocity model by using Newton's method to obtain the final velocity model includes:

[0049] A forward wavefield data acquisition module, used to determine synthetic observation data and forward wavefield data according to the velocity model;

[0050] A back propagation wavefield data acquisition module, used for obtaining back propagation wavefield data according to the velocity model and actual observation data;

[0051] An objective function building module, used for determining an objective equation of a residual according to the forward wavefield data and the reverse wavefield data;

[0052] A gradient matrix calculation module, used to calculate a gradient matrix about wave velocity using the forward wave field data and the reverse wave field data;

[0053] A Hessian matrix calculation module, used for calculating the Hessian matrix based on the gradient matrix;

[0054] A model updating module, used for updating the velocity model according to the gradient matrix, the Hessian matrix and a preset iteration step size;

[0055] The iterative judgment module is used to repeatedly execute the forward wave field data acquisition module to the model updating module until the preset convergence criterion is met to obtain the final velocity model.

[0056] A storage medium stores a program, and when the program is executed by a processor, the above-mentioned method for advance detection of karst caves by tunnel seismic detection is implemented.

[0057] A computing device includes a processor and a memory for storing a program executable by the processor. When the processor executes the program stored in the memory, the above-mentioned method for advance detection of karst caves by tunnel seismic detection is implemented.

[0058] Compared with the prior art, the present invention has the following beneficial effects:

[0059] (1) Compared with the traditional velocity analysis method, the full waveform inversion adopted by the present invention has relatively high accuracy and provides a more accurate velocity model for migration imaging;

[0060] (2) In the present invention, since the overall model of tunnel advance detection is relatively small in scale and relatively simple in structure, this provides favorable conditions for the storage and inverse matrix of the Hessian matrix. Compared with the conjugate gradient method used in the traditional full waveform inversion, the Newton method of the present invention has a faster convergence speed and a higher overall inversion speed efficiency;

[0061] (3) The present invention arranges detectors in front of the tunnel face and inside the rock mass on both sides, and forms an observation system with a relatively large offset distance with the earthquake source, which effectively solves the problem of loss of lateral resolution at zero offset distance, and fundamentally solves the difficult problem of symmetry illusion in the advance detection results. The seismic records received under the observation system have good lateral resolution, providing effective wave field information for full waveform inversion;

[0062] (4) The present invention maintains the stability and efficiency of the algorithm during the process of solving the objective function under regularization constraints. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 A flowchart of the steps of the method for advance detection of karst caves in tunnels according to the present invention;

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

[0065] Figure 3 The two-component seismic record of the numerical simulation of the acoustic wave for tunnel advance detection, Figure a is the X-separation seismic record, and Figure b is the Y-component seismic record;

[0066] Figure 4 This is the velocity result diagram of full waveform inversion.

[0067] In the figure, 1-seismic source, 2-detector, 3-tunnel, 4-first fault, 5-second fault, 6-direct longitudinal wave, 7-first reflected longitudinal wave, 8-second reflected longitudinal wave. DETAILED DESCRIPTION

[0068] 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. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0069] 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.

[0070] 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.

[0071] 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.

[0072] 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.

[0073] Example 1

[0074] See also Figure 1 and Figure 2 , Figure 1 A flowchart of the steps of the method for advance detection of karst caves in tunnels according to the present invention; Figure 2 The schematic diagram of the earthquake advance detection system in the tunnel earthquake advance detection method of the present invention is as follows. A tunnel earthquake advance detection method of karst caves includes the following steps:

[0075] S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, the multiple geophones are located between the seismic sources and the tunnel face, and arrange multiple geophones at equal intervals on the tunnel face, with each geophone and each seismic source located at the same height;

[0076] S2, connecting each seismic source and all geophones through equipment to form an earthquake advance detection system, exciting the seismic sources in turn, and obtaining actual observation data through the earthquake advance detection system;

[0077] S3, establishing an initial velocity model, and updating the velocity model using Newton's method to obtain a final velocity model;

[0078] The steps of updating the velocity model using Newton's method to obtain the final velocity model include:

[0079] S31, determining synthetic observation data and forward wave field data according to the velocity model;

[0080] S32, obtaining back propagation wave field data according to the velocity model and actual observation data;

[0081] S33, determining a target equation of the residual according to the forward wavefield data and the reverse wavefield data;

[0082] S34, calculating a gradient matrix related to wave velocity using the forward wave field data and the reverse wave field data;

[0083] S35, calculating the Hessian matrix based on the gradient matrix;

[0084] S36, updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size;

[0085] S37, repeating steps S31 to S36 until a preset convergence criterion is met to obtain a final velocity model.

[0086] In the above step S1, in this embodiment, S seismic sources are arranged at the waistline position of the left side wall of the tunnel, and the distance between each seismic source is h meters. When arranging the seismic sources, a borehole with a depth of e meters is dug at the seismic source arrangement position on the tunnel side wall, and the seismic source is placed in the corresponding borehole. Then, along the tunnel excavation direction, N geophones are arranged behind the seismic source, and the track spacing between the geophones is d meters. The distance between the first geophone and its adjacent seismic source is g meters. The first geophone position is (N-1)*d meters away from the tunnel face position. When arranging the geophones, a borehole with a depth of a meter is dug at the geophone arrangement position on the tunnel side wall, and then the geophone is placed in the corresponding borehole. Finally, M geophones are arranged horizontally on the tunnel face, with a track spacing of d meters between the geophones. The arrangement is consistent with the arrangement of the sidewall geophones, and the distance between the geophone closest to the left wall of the tunnel on the tunnel face and the last geophone on the left wall of the tunnel is d meters. All geophones and seismic sources are located at the same height, such as Figure 2 As shown. Among them, the specific values ​​of S, h, e, N, d, g, a and M can be determined according to the actual situation. Furthermore, the detector can be a dual-component detector, the source is a Ricker wavelet, the main frequency is 60Hz, and the absorption boundary length is 50 meters.

[0087] In the above step S2, for each earthquake source, the earthquake source and all the detectors are connected through the equipment to form an earthquake advance detection system. Then, a two-dimensional coordinate system is established with the first earthquake source position as the coordinate origin, wherein the tunnel excavation direction is the X direction and the vertical direction is the Y direction, and the earthquake source position and each detector position are included in the established two-dimensional coordinate system. Then, the detector sampling interval is set to 1ms, the sampling length is c, the explosive amount is set to Q, the earthquake source is stimulated in turn, and the seismic records of all two-component detectors are obtained. Then, conventional methods are used to remove other signals except the longitudinal wave to obtain the actual multi-shot observation data D (S; x1, y1; x2, y2; ... x M+N ,y M+N ; t), N+M is the total number of detectors, (x1, y1; x2, y2; ...x M+N ,y M+N ) is the location parameter of the detector, t is the time, S is the number of sources, one source corresponds to N+M detectors, and the combination and arrangement form the actual observation data. Specifically, the sampling interval of the detector can be set to 0.001s and the sampling length can be set to 1000.

[0088] In the above step S3, an initial velocity model v0(x, y) is established, and the initial velocity model v0(x, y) = 3000 m / s. Then the velocity model is updated using Newton's method to obtain a final velocity model, and the final velocity model is used for full waveform inversion.

[0089] In the above step S31, the simulation parameters are determined according to the arrangement of the earthquake advance detection system in steps S1 and S2, and the velocity model is numerically simulated to obtain the synthetic observation data s n (x, y, t) and the forward wave field data p n (x,y,t), where the synthetic observation data s n (x, y, t) and the forward wave field data p n (x, y, t) are based on the velocity model v n (x, y), n represents the number of iterations.

[0090] Further, in step S31, the step of determining synthetic observation data and forward wavefield data according to the velocity model includes:

[0091] S311, determining observation system parameters according to the earthquake advance detection system, and calculating forward wave field data of different earthquake sources using a spatiotemporal high-order finite difference algorithm based on the observation system parameters;

[0092] L[v n (x,y)]p n (x,y,t)=S(x S ,y S ,t)(1);

[0093] Where L[·] is the high-order finite difference operator of acoustic wave space-time, v n (x, y) is the velocity model after the nth iteration, v0(x, y) is the initial input of the uniform model, and then v0(x, y) is updated to obtain v1(x, y), v2(x, y), v3(x, y)…v n (x,y),p n (x,y,t) is based on v n (x,y) forward wave field data, S(x S ,y S ,t) is the source wavelet, (x S ,y S ) is the earthquake source location, t is the time parameter;

[0094] S312: Determine synthetic observation data based on the forward wavefield data.

[0095] In the above steps S311 and S312, the observation system is set according to the arrangement of the earthquake advance detection system in step S1, that is, the arrangement of the observation system is consistent with the arrangement of the actual earthquake advance detection system, and then the observation system parameters are substituted into the time-space high-order finite difference algorithm to perform numerical simulation on the velocity model to obtain the forward wave field data p n (x,y,t), the forward wave field data p n The synthetic two-component seismic record s can be obtained by extracting all the time of the geophone position in (x, y, t). n (x, y, t), that is, the synthetic observation data s n (x, y, t). Compared with the traditional finite difference method, the high-order finite difference method in time and space has higher calculation stability and faster calculation efficiency. The forward wave field data in the synthetic observation data is calculated by formula (1).

[0096] Further, in step S311, the observation system parameters are determined according to the earthquake advance detection system, and based on the observation system parameters, the steps of using a spatiotemporal high-order finite difference algorithm to calculate the forward wavefield data of different seismic sources include:

[0097] S3111, obtaining a velocity model;

[0098] S3112, set the forward modeling parameters, the grid spacing is dx, dy, the source position is (xs, ys), use the Ricker wavelet for forward modeling, the main frequency is f, arrange the observation system, the arrangement of the observation system is consistent with the earthquake advance detection system, and use the time-space high-order finite difference algorithm for forward simulation;

[0099] S3113, using the absorbing boundary condition to eliminate boundary reflection, the thickness is set to H meters;

[0100] S3114. Solve the wave equation using a high-order finite difference algorithm in time and space, as follows:

[0101] (1) Establish a two-dimensional pseudo-velocity-stress acoustic wave equation:

[0102]

[0103] Where λ and μ are Lamé parameters, ρ is density, λ+2μ=ρα 2 , μ=ρβ 2 , α and β are the velocity components of the longitudinal and transverse waves, vx and vy are the velocity components of the particle vibration, τ 12 , τ 111 , τ 111x and τ 111y is the stress component of the particle vibration, is the partial derivative with respect to time t;

[0104] (2) Let X = α or X = β, the time derivative The following discretization is used:

[0105]

[0106] Where p represents the wave field variable, Δt is the time discrete step length;

[0107] (3) Difference coefficient The calculation of is as follows:

[0108]

[0109] In the formula, h φ represents the spatial discrete step length in the φ direction, express The spatial discrete step length in the axial direction;

[0110] (4) The stability conditions of the difference format are as follows:

[0111]

[0112] In the formula,

[0113] (5) The following high-order time-space difference operator is used to calculate the numerical simulation of the entire wave field p n (x,y,t).

[0114]

[0115] Where φ represents the orthogonal coordinate axis other than the derivative axis r, r, φ∈(x,z), r≠φ, h r The axis represents the discrete step length in the r direction, N r Represents half the difference length.

[0116] The above process is the calculation process of the L[·] algorithm in formula (1).

[0117] In step S32, the step of obtaining reverse wavefield data based on the velocity model and actual observation data includes:

[0118] S331, using the actual observation data as the boundary condition, and using formula (2) to calculate the reverse wave field data d n (x,y,t).

[0119] L[v n (x,y)]d n (x,y,t)=D(S;x1,y1;x2,y2;…x M+N ,y M+N ;t)(2);

[0120] Where, d n (x,y,t) is based on v n (x,y) reverse wavefield data,

[0121] D(S;x1,y1;x2,y2;…x M+N ,y M+N ; t) is the actual observation data,

[0122] (x1,y1; x2,y2; ...x M+N ,y M+N ) is the location parameter of the geophone, N+M represents the total number of geophones, t is the time, and S represents the number of sources.

[0123] In the above step S33, the expression of the target equation is:

[0124]

[0125] In the formula, F[v n (x,y)] is the objective function value after the nth iteration, v n (x, y) is the velocity model after the nth iteration, p n (x,y,t) is based on v n (x, y) forward wave field data, d n (x,y,t) is based on v n (x, y) is the reverse wave field data, ‖·‖ is the two-norm, ζ is the regularization scale factor, W is the regularization weight operator, v apr is the expected approximate solution.

[0126] In the above step S34, the forward wave field data p n (x, y., t) and reverse wave field data d n (x, y., t) calculates the gradient matrix about the wave velocity. The expression of the gradient matrix is:

[0127]

[0128] In the formula, F[v n (x,y)] about v n The first-order partial derivative of (x,y), v n (x,y) is the velocity model after the nth iteration.

[0129] In the above step S35, based on the gradient matrix g[v n (x,y)] calculates the Hessian matrix, and the expression of the Hessian matrix is:

[0130]

[0131] In the formula, k n is the gradient update amount, k n =g[v n+1 (x,y)]-g[v n (x,y)],q n For the velocity model v n The update amount of (x,y), q n =v n+1 (x,y)-v n (x,y), represents the matrix transpose, is the initial positive definite Hessian matrix of the nth iteration, For Based on this, the approximate Hessian matrix of the jth recursive iteration is performed, express About v n The first-order partial derivative of (x,y), that is, F[v n (x,y)] about v n The second-order partial derivative of (x,y).

[0132] In step S36, the step of updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size includes:

[0133] S361, update the velocity model using formula (6), the expression of formula (6) is:

[0134] v n+1 (x,y)=v n (x,y)+a n p n (6);

[0135] In the formula, a n is the iteration step length, v n (x,y) is the nth iteration speed model, v n+1 (x, y) is the velocity model after the n+1th iteration, p n is the iteration direction, p n The expression is;

[0136]

[0137] Where γ is the damping coefficient and I is the unit matrix.

[0138] In the above step S361, if the initial velocity model v0(x, y) is substituted into formula (6) for updating, n is 0 at this time, so as to obtain the updated velocity model v1(x, y), and then continuously iterate and update, v n(x,y) is given by v n-1 (x,y) iterates to get, v n+1 (x,y) is given by v n (x,y) is obtained by iteration.

[0139] In the above step S37, after the initial velocity model is updated, a new velocity model is obtained, and after the initial velocity model is replaced by the new velocity model, steps S31 to S36 are repeated until the preset convergence criterion is met to obtain the final velocity model. The preset convergence criterion may be, for example, that the objective function value F[v(x,y)] is less than or equal to a preset threshold value, and the preset threshold value may be set to a minimum value ε0, ε0=0.000001. For example, the initial velocity model v0(x,y) is updated to obtain a new velocity model v1(x,y), and then the new velocity model v1(x,y) is used as the initial velocity model. The objective function value F[v(x,y)] is obtained through steps S31 and S32, and it is determined whether the objective function value F[v(x,y)] is less than or equal to ε0. If so, the velocity model v1(x,y) at this time is used as the final velocity model output of the full waveform inversion; otherwise, continue to execute steps S33 to S36, update the velocity model v1(x,y) to obtain a new velocity model v2(x,y), and repeat the iteration until the objective function value F[v(x,y)] is less than or equal to ε0.

[0140] Example 2

[0141] A tunnel seismic advance detection system for karst caves, comprising:

[0142] A detection system erection unit, used to arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are located between the seismic sources and the tunnel face, and multiple geophones are arranged at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height;

[0143] An observation data acquisition unit, used to connect each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, stimulate the seismic sources in turn, and acquire actual observation data through the earthquake advance detection system;

[0144] A velocity model updating unit, used for establishing an initial velocity model, updating the velocity model using Newton's method, and obtaining a final velocity model;

[0145] Wherein, the updating of the velocity model by using Newton's method to obtain the final velocity model includes:

[0146] A forward wavefield data acquisition module, used to determine synthetic observation data and forward wavefield data according to the velocity model;

[0147] A back propagation wavefield data acquisition module, used for obtaining back propagation wavefield data according to the velocity model and actual observation data;

[0148] An objective function building module, used for determining an objective equation of a residual according to the forward wavefield data and the reverse wavefield data;

[0149] A gradient matrix calculation module, used to calculate a gradient matrix about wave velocity using the forward wave field data and the reverse wave field data;

[0150] A Hessian matrix calculation module, used for calculating the Hessian matrix based on the gradient matrix;

[0151] A model updating module, used for updating the velocity model according to the gradient matrix, the Hessian matrix and a preset iteration step size;

[0152] The iterative judgment module is used to repeatedly execute the forward wave field data acquisition module to the model updating module until the preset convergence criterion is met to obtain the final velocity model.

[0153] Example 3

[0154] A storage medium stores a program, and when the program is executed by a processor, the method for advanced detection of karst caves in a tunnel seismic detection described in Example 1 is implemented as follows:

[0155] S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are located between the seismic sources and the tunnel face, and arrange multiple geophones at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height;

[0156] S2, connecting each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, exciting the seismic sources in turn, and obtaining actual observation data through the earthquake advance detection system;

[0157] S3, establishing an initial velocity model, and updating the velocity model using Newton's method to obtain the final velocity model;

[0158] The step of updating the velocity model using Newton's method to obtain the final velocity model includes:

[0159] S31, determining synthetic observation data and forward wavefield data according to the velocity model;

[0160] S32, obtaining back propagation wavefield data according to the velocity model and actual observation data;

[0161] S33, determining a target equation of the residual according to the forward wavefield data and the reverse wavefield data;

[0162] S34, calculating a gradient matrix related to wave velocity using the forward wave field data and the reverse wave field data;

[0163] S35, calculating a Hessian matrix based on the gradient matrix;

[0164] S36, updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size;

[0165] S37, repeating steps S31 to S36 until a preset convergence criterion is met to obtain a final velocity model.

[0166] In the above process, the specific processing process is as described in Example 1 and will not be repeated here.

[0167] In this embodiment, the storage medium may be a magnetic disk, an optical disk, a computer memory, a read-only memory, a random access memory, a USB flash drive, a mobile hard disk, or the like.

[0168] Example 4

[0169] A computing device includes a processor and a memory for storing a program executable by the processor. When the processor executes the program stored in the memory, the method for advanced detection of karst caves in a tunnel by seismic advance detection described in Example 1 is implemented as follows:

[0170] S1. Arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are located between the seismic sources and the tunnel face, and arrange multiple geophones at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height;

[0171] S2, connecting each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, exciting the seismic sources in turn, and obtaining actual observation data through the earthquake advance detection system;

[0172] S3, establishing an initial velocity model, and updating the velocity model using Newton's method to obtain the final velocity model;

[0173] The step of updating the velocity model using Newton's method to obtain the final velocity model includes:

[0174] S31, determining synthetic observation data and forward wavefield data according to the velocity model;

[0175] S32, obtaining back propagation wavefield data according to the velocity model and actual observation data;

[0176] S33, determining a target equation of the residual according to the forward wavefield data and the reverse wavefield data;

[0177] S34, calculating a gradient matrix related to wave velocity using the forward wave field data and the reverse wave field data;

[0178] S35, calculating a Hessian matrix based on the gradient matrix;

[0179] S36, updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size;

[0180] S37, repeating steps S31 to S36 until a preset convergence criterion is met to obtain a final velocity model.

[0181] In the above process, the specific processing process is as described in Example 1 and will not be repeated here.

[0182] In this embodiment, the computing device may be a terminal device such as a desktop computer, a laptop computer, a PDA handheld terminal, a tablet computer, etc.

[0183] The above embodiments are preferred implementations of the present invention and are not intended to limit the present invention. Any other changes or other equivalent replacement methods that do not deviate from the technical solutions of the present invention are included in the protection scope 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 located between the seismic sources and the tunnel face, and arrange multiple geophones at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height; S2, connecting each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, exciting the seismic sources in turn, and acquiring actual observation data through the earthquake advance detection system; S3, establishing an initial velocity model, and updating the velocity model using Newton's method to obtain the final velocity model; The step of updating the velocity model using Newton's method to obtain the final velocity model includes: S31, determining synthetic observation data and forward wavefield data according to the velocity model; S32, obtaining back propagation wave field data according to the velocity model and actual observation data; S33, determining a target equation of residual according to the forward wavefield data and the reverse wavefield data; S34, calculating a gradient matrix related to wave velocity using the forward wave field data and the reverse wave field data; S35, calculating a Hessian matrix based on the gradient matrix; S36, updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size; S37, repeating steps S31 to S36 until a preset convergence criterion is met to obtain a final velocity model.

2. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: The step of determining synthetic observation data and forward wavefield data according to the velocity model comprises: Determine observation system parameters according to the earthquake advance detection system, and calculate the forward wave field data based on the observation system parameters using a spatiotemporal high-order finite difference algorithm: L[v n (x,y)]p n (x,y,t)=S(x S ,y S ,t) (1); Where L[·] is the high-order finite difference operator of acoustic wave space-time, v n (x, y) is the velocity model after the nth iteration, p n (x,y,t) is based on v n (x,y) forward wave field data, S(x S ,y S ,t) is the source wavelet, (x S ,y S ) is the earthquake source location, t is the time parameter; Synthetic observation data is determined based on the forward wavefield data.

3. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: The step of obtaining back propagation wavefield data according to the velocity model and actual observation data comprises: The actual observed data is used as the boundary condition and the reverse wave field data is calculated using formula (2): L[v n (x,y)]d n (x,y,t)=D(S;x1,y1;x2,y2;…x M+N ,y M+N ;t) (2); In the formula, v n (x, y) is the velocity model after the nth iteration, d n (x,y,t) is based on v n (x, y) reverse wave field data, D(S; x1, y1; x2, y2; ...x M+N ,y M+N ; t) is the actual observed data, (x1, y1; x2, y2; ...x M+N ,y M+N ) is the location parameter of the detector, N+M represents the total number of detectors, t is the time, and S represents the number of sources.

4. The method for advanced detection of karst caves by tunnel seismic detection according to claim 1, characterized in that: In the step of determining the target equation of the residual according to the forward wavefield data and the backward wavefield data, the expression of the target equation is: In the formula, F[v n (x,y)] is the objective function value after the nth iteration, v n (x, y) is the velocity model after the nth iteration, p n (x,y,t) is based on v n (x, y) forward wave field data, d n (x,y,t) is based on v n (x, y) is the reverse wave field data, ‖·‖ is the two-norm, ζ is the regularization scale factor, W is the regularization weight operator, v apr is the expected approximate solution.

5. The method for advanced detection of karst caves by tunnel seismic detection according to claim 4, characterized in that: In the step of calculating the gradient matrix about wave velocity using the forward wave field data and the reverse wave field data, the expression of the gradient matrix is: In the formula, F[v n (x,y)] about v n The first-order partial derivative of (x,y), v n (x,y) is the velocity model after the nth iteration.

6. The method for advanced detection of karst caves by tunnel seismic detection according to claim 5, characterized in that: In the step of calculating the Hessian matrix based on the gradient matrix, the expression of the Hessian matrix is: In the formula, k n is the gradient update amount, k n =g[v n+1 (x,y)]-g[v n (x,y)],q n For the velocity model v n The update amount of (x,y), q n =v n+1 (x,y)-v n (x,y), T represents matrix transpose, is the initial positive definite Hessian matrix of the nth iteration, For Based on this, the approximate Hessian matrix of the jth recursive iteration is performed, express About v n The first partial derivative of (x,y).

7. The method for advanced detection of karst caves by tunnel seismic detection according to claim 6, characterized in that: The step of updating the velocity model according to the gradient matrix, the Hessian matrix and the preset iteration step size comprises: The velocity model is updated using formula (6), and the expression of formula (6) is: v n+1 (x,y)=y n (x,y)+a n p n (6); In the formula, a n is the iteration step length, v n (x, y) is the velocity model after the nth iteration, v n+1 (x, y) is the velocity model after the n+1th iteration, p n is the iteration direction, p n The expression is; Where γ is the damping coefficient and I is the unit matrix.

8. A tunnel seismic advance detection system for karst caves, characterized in that: include: A detection system erection unit, used to arrange multiple seismic sources and multiple geophones at equal intervals on one side wall of the tunnel, wherein the multiple geophones are located between the seismic sources and the tunnel face, and multiple geophones are arranged at equal intervals on the tunnel face, wherein each geophone and each seismic source is located at the same height; An observation data acquisition unit, used to connect each of the seismic sources and all the geophones through devices to form an earthquake advance detection system, stimulate the seismic sources in turn, and acquire actual observation data through the earthquake advance detection system; A velocity model updating unit, used for establishing an initial velocity model, updating the velocity model using Newton's method, and obtaining a final velocity model; Wherein, the updating of the velocity model by using Newton's method to obtain the final velocity model includes: A forward wavefield data acquisition module, used for determining synthetic observation data and forward wavefield data according to the velocity model; A back propagation wavefield data acquisition module, used for obtaining back propagation wavefield data according to the velocity model and actual observation data; An objective function building module, used for determining an objective equation of a residual according to the forward wavefield data and the reverse wavefield data; A gradient matrix calculation module, used to calculate a gradient matrix about wave velocity using the forward wave field data and the reverse wave field data; A Hessian matrix calculation module, used for calculating the Hessian matrix based on the gradient matrix; A model updating module, used for updating the velocity model according to the gradient matrix, the Hessian matrix and a preset iteration step size; The iterative judgment module is used to repeatedly execute the forward wave field data acquisition module to the model updating module until the preset convergence criterion is met to obtain the final velocity model.

9. A storage medium, characterized in that: A program is stored, and when the program is executed by a processor, a method for advance detection of karst caves in tunnels by seismic advance detection of karst caves in tunnels as described in any one of claims 1 to 7 is implemented.

10. A computing device, characterized in that It comprises a processor and a memory for storing a program executable by the processor. When the processor executes the program stored in the memory, a method for advance detection of karst caves in a tunnel by seismic advance detection as described in any one of claims 1 to 7 is implemented.

Images