A tunnel earthquake advanced detection cave method, system, medium and equipment

Abstract

The application discloses a tunnel earthquake advanced detection cave method, system, medium and equipment, and the method comprises the following steps: S1, a plurality of seismic sources and a plurality of geophones are arranged at equal intervals on one side wall of the tunnel, the plurality of geophones are located between the seismic sources and the tunnel face, and a plurality of geophones are arranged at equal intervals on the tunnel face, and the height of the geophones is consistent with that of the seismic sources; S2, each seismic source and all geophones are connected by equipment to form an earthquake advanced detection system, the seismic sources are sequentially excited, and actual observation data are acquired; S3, an initial velocity model is established, the Newton method is used to update the velocity model, and a final velocity model is obtained. The application combines the large-azimuth earthquake advanced observation system with the full waveform inversion method, overcomes the problem of low horizontal resolution of the traditional zero-offset, simultaneously uses the Newton method to replace the conjugate gradient method of the traditional full waveform inversion, further improves the calculation efficiency of the full waveform inversion, and is helpful to improve the tunnel advanced imaging precision.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of tunnel detection, and particularly relates to a tunnel seismic advanced detection cave method, system, medium and equipment. BACKGROUND

[0002] The tunnel excavation stage is a frequent stage of tunnel safety accidents, and the main reason is that the structure problem in front is not clear, among which the fault structure is the most typical, so the accurate detection of the fault has important significance. Since the drilling cost is relatively high, the tunnel advanced detection reflection seismic exploration method is relatively practical and convenient, so the seismic reflection method is the most common in the tunnel advanced detection. In the seismic reflection advanced detection work, how to establish an accurate fault velocity model is the key to realize the seismic advanced detection, and the accurate fault velocity model provides favorable conditions for realizing the advanced detection of the fault in front of the tunnel.

[0003] Full waveform inversion is a waveform inversion method based on data space domain, which matches the forward simulation data with the actual observation data, establishes an objective function through the wave field error minimization, and then finds the best model parameter to make the simulation data and the observation data best match. For the theory of full waveform inversion, it was first established on the basis of data domain fitting under the constraint of generalized least squares proposed by Tarantola (1982). Tarantola (1984, 1988) applied this theory to acoustic approximation, and then gave a complete time domain full waveform inversion theoretical framework. This technology is relatively mature in ground seismic exploration, and the full waveform inversion has much higher accuracy than the traditional velocity inversion method. However, the conjugate gradient method is currently used to optimize the model velocity parameter in the full waveform inversion, and the optimization speed of this method is relatively slow, thereby leading to the low calculation efficiency of the current full waveform inversion. SUMMARY

[0004] The first object of the present application is to provide a tunnel seismic advanced detection cave method, which combines a large azimuth observation system with a full waveform inversion method, overcomes the problem of low traditional zero offset lateral resolution, and further improves the calculation efficiency of the full waveform inversion by using the Newton method to replace the conjugate gradient method of the traditional full waveform inversion, thereby helping to improve the tunnel advanced detection imaging accuracy.

[0005] The second object of the present application is to provide a tunnel seismic advanced detection cave system.

[0006] The third object of the present application is to provide a storage medium.

[0007] The fourth object of the present application is to provide a computing device.

[0008] The purpose of the present application is achieved by the following technical solutions: A tunnel earthquake advanced detection cave method, comprising the following steps:

[0009] S1, a plurality of seismic sources are arranged at equal intervals on one side wall of the tunnel, and a plurality of geophones are arranged at equal intervals, the plurality of geophones are located between the seismic sources and the working face, a plurality of geophones are arranged at equal intervals on the working face, and each geophone and each seismic source are located at the same height;

[0010] S2, each seismic source and all geophones are connected by a device to form an earthquake advanced detection system, the seismic sources are sequentially excited, and actual observation data are obtained through the earthquake advanced detection system;

[0011] S3, an initial velocity model is established, the Newton method is used to update the velocity model, and a final velocity model is obtained;

[0012] The step of updating the velocity model by the Newton method to obtain the final velocity model comprises:

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

[0014] S32, the backward wave field data are obtained according to the velocity model and the actual data;

[0015] S33, the residual target equation is determined according to the forward wave field data and the backward wave field data;

[0016] S34, the gradient matrix about wave velocity is calculated by using the forward wave field data and the backward wave field data;

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

[0018] S36, the velocity model is updated according to the gradient matrix, the Hessian matrix and a preset iteration step;

[0019] S37, steps S31 to S36 are repeatedly executed until a preset convergence criterion is met, and the final velocity model is obtained.

[0020] Further, the step of determining the synthetic observation data and the forward wave field data according to the velocity model comprises:

[0021] The observation system parameters are determined according to the earthquake advanced detection system, the forward wave field data are calculated based on the observation system parameters by using a 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] wherein L[·] is a high-order finite difference operator of acoustic wave in space-time, v n (x,y) is a velocity model after the n th iteration, p n (x,y,t) is forward wavefield data based on v n (x,y), S(x S ,y S ,t) is a source wavelet, (x S ,y S ) is a source position, and t is a time parameter.

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

[0025] Further, the step of obtaining the backward wavefield data based on the velocity model and the actual observation data comprises:

[0026] The actual observation data is taken as a boundary condition, and the backward wavefield data is calculated by 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] wherein v n (x,y) is a velocity model after the n th iteration, d n (x,y,t) is backward wavefield data based on v n (x,y), D(S;x1,y1;x2,y2;…x M+N ,y M+N ;t) is actual observation data, (x1,y1;x2,y2;…x M+N ,y M+N ) is a position parameter of a receiver, N+M represents a total number of receivers, t is time, and S represents a number of sources.

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

[0030]

[0031] wherein F[v n (x,y)] is a target function value after the n th iteration, v n (x,y) is a velocity model after the n th iteration, p n (x,y,t) is forward wavefield data based on v n (x,y), d n(x,y,t) is based on v n (x,y), ‖·‖ is the two-norm, ζ is a regularization scaling factor, W is a regularization weight operator, v apr is the desired approximate solution.

[0032] Further, in the step of calculating the gradient matrix with respect to wave velocity using the forward wavefield data and the backward wavefield data, the expression of the gradient matrix is:

[0033]

[0034] wherein, is the first-order partial derivative of F[v n (x,y)] with respect to v n (x,y), v n (x,y) is the velocity model after the nth iteration.

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

[0036]

[0037] wherein, k n is the update amount of the gradient, k n = g[v n+1 (x,y)] - g[v n (x,y)], q n is the update amount of the velocity model v n (x,y), q n = v n+1 (x,y) - v n (x,y), T represents the matrix transpose, is the initial positive definite Hessian matrix of the nth iteration, is the approximate Hessian matrix after the jth recursive iteration based on , and represents the first-order partial derivative with respect to v n (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 comprises:

[0039] The velocity model is updated using formula (6), and the expression of formula (6) is:

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

[0041] wherein a n is the iteration step, v n (x,y) is the velocity model after the nth iteration, v n+1 (x,y) is the velocity model after the (n+1)th iteration, p n is the iteration direction, p n The expression of p is:

[0042]

[0043] wherein γ is a damping coefficient, and I is a unit matrix.

[0044] A tunnel seismic advanced detection cave system, comprising:

[0045] A detection system erection unit for arranging a plurality of seismic sources at equal intervals on a side wall of a tunnel and a plurality of geophones at equal intervals, the plurality of geophones being located between the seismic sources and a working face, a plurality of geophones being arranged at equal intervals on the working face, each geophone and each seismic source being located at the same height;

[0046] An observation data acquisition unit for connecting each seismic source and all geophones through equipment to form a seismic advanced detection system, sequentially exciting the seismic sources, and acquiring actual observation data through the seismic advanced detection system;

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

[0048] The updating of the velocity model by using the Newton method and the obtaining of the final velocity model include:

[0049] A forward wave field data acquisition module for determining synthetic observation data and forward wave field data according to the velocity model;

[0050] A reverse wave field data acquisition module for obtaining reverse wave field data according to the velocity model and actual observation data;

[0051] A target function construction module for determining a target equation of a residual according to the forward wave field data and the reverse wave field data;

[0052] A gradient matrix calculation module for calculating a gradient matrix of wave velocity by using the forward wave field data and the reverse wave field data;

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

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

[0055] an iteration judging module, configured to repeatedly execute the forward wave field data obtaining module to the model updating module until a preset convergence criterion is met, and obtain a final velocity model.

[0056] A storage medium storing a program, when executed by a processor, implements the tunnel seismic advanced detection cave method.

[0057] A computing device comprising a processor and a memory for storing a program executable by the processor, when the processor executes the program stored in the memory, implements the tunnel seismic advanced detection cave method.

[0058] Compared with the prior art, the tunnel seismic advanced detection cave method has the following advantages:

[0059] (1) Compared with the traditional velocity analysis method, the full waveform inversion method adopted in the tunnel seismic advanced detection cave method has relatively high precision, and can provide a more accurate velocity model for migration imaging;

[0060] (2) In the tunnel seismic advanced detection cave method, the overall model scale of the tunnel advanced detection is relatively small, and the structure is relatively simple, which provides favorable conditions for Hessian matrix storage and inverse matrix, and compared with the conjugate gradient method adopted in the traditional full waveform inversion method, the convergence speed of the Newton method of the tunnel seismic advanced detection cave method is faster, and the efficiency of the overall inversion speed is faster;

[0061] (3) In the tunnel seismic advanced detection cave method, the geophones are arranged in front of the tunnel face and inside the rock mass on both sides, and form a relatively large offset observation system with the seismic source, effectively solving the problem of loss of lateral resolution at zero offset, and fundamentally solving the difficult problem of symmetric false image in the advanced detection result. The seismic record received under the observation system has good lateral resolution, and provides effective wave field information for full waveform inversion;

[0062] (4) The tunnel seismic advanced detection cave method maintains the stability and efficiency of the algorithm in the process of solving the objective function under the regularization constraint. BRIEF DESCRIPTION OF DRAWINGS

[0063] Figure 1 The figure is a step flow chart of the tunnel seismic advanced detection cave method;

[0064] Figure 2 The figure is a schematic diagram of the seismic advanced detection system in the tunnel seismic advanced detection cave method;

[0065] Figure 3 The figure is a two-component seismic record of tunnel advanced detection acoustic wave numerical simulation, figure a is an X-component seismic record, and figure b is a Y-component seismic record;

[0066] Figure 4 to obtain a full waveform inversion velocity result map.

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

[0068] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some but not all of the embodiments of the present application. The components of the embodiments of the present application described and shown in the drawings herein can be arranged and designed in various different configurations.

[0069] Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present application without creative work fall within the scope of protection of the present application.

[0070] It should be noted that: similar reference numerals and letters represent similar items in the following drawings, therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings. Meanwhile, in the description of the present application, the terms “first”, “second” and the like are only used to distinguish description, and cannot be understood as indicating or implying relative importance.

[0071] It should be noted that: similar reference numerals and letters represent similar items in the following drawings, therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings. Meanwhile, in the description of the present application, the terms “first”, “second” and the like are only used to distinguish description, and cannot be understood as indicating or implying relative importance.

[0072] In the description of the present application, it should be noted that the terms "upper", "lower", "inner", "outer" and the like indicate the orientation or positional relationship shown in the drawings, or the orientation or positional relationship commonly used when the product of the present application is used, and are only for the convenience of describing the present application and simplifying the description, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as a limitation on the present application.

[0073] Embodiment 1

[0074] Please refer to Figure 1 and Figure 2 , Figure 1 is a flow chart of the steps of the tunnel seismic advanced detection cave method of the present application; Figure 2 is a schematic diagram of the seismic advanced detection system in the tunnel seismic advanced detection cave method of the present application. A tunnel seismic advanced detection cave method comprises the following steps:

[0075] S1, arranging a plurality of seismic sources at equal intervals on one side wall of the tunnel and a plurality of geophones at equal intervals, the plurality of geophones being located between the seismic sources and the working face, arranging a plurality of geophones at equal intervals on the working face, each geophone and each seismic source being located at the same height;

[0076] S2, connecting each seismic source and all geophones through equipment to form a seismic advanced detection system, sequentially exciting the seismic sources, and acquiring actual observation data through the seismic advanced detection system;

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

[0078] The step of updating the velocity model using Newton's method to obtain a final velocity model comprises:

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

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

[0081] S33, determining a residual target equation according to the forward wave field data and the backward wave field data;

[0082] S34, calculating a gradient matrix with respect to wave velocity using the forward wave field data and the backward wave field data;

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

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

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

[0086] In step S1, in the present embodiment, S sources are arranged at the waistline position of the left sidewall of the tunnel, and the distance between each source is h meters. When arranging the sources, drill holes with a depth of e meters are dug at the source arrangement positions on the sidewall of the tunnel, and the sources are placed in the corresponding drill holes. Then, N receivers are arranged behind the sources along the tunnel excavation direction, and the trace interval between the receivers is d meters. The distance between the first receiver and its adjacent source is g meters. The first receiver is located at a distance of (N-1)*d meters from the position of the tunnel face. When arranging the receivers, drill holes with a depth of a meters are dug at the receiver arrangement positions on the sidewall of the tunnel, and the receivers are placed in the corresponding drill holes. Finally, M receivers are arranged along the horizontal direction at the tunnel face, and the trace interval between the receivers is d meters. The arrangement mode is consistent with that of the sidewall receivers. The distance between the receiver closest to the left sidewall of the tunnel on the tunnel face and the last receiver on the left sidewall of the tunnel is d meters. Each receiver and each source are located at the same height, as shown in FIG. 1. In the figure, S, h, e, N, d, g, a, and M are specific values that can be determined according to actual conditions. Further, the receivers can be two-component receivers, and the sources can be Ricker wavelets with a main frequency of 60 Hz and an absorption boundary length of 50 meters. Figure 2

[0087] In step S2, for each source, the source and all receivers are connected by equipment to form a seismic advanced detection system. Then, a two-dimensional coordinate system is established with the first source position as the origin. The tunnel excavation direction is the X direction, and the vertical direction is the Y direction. The source position and each receiver position are assigned to the established two-dimensional coordinate system. Then, the receiver sampling interval is set to lms, the sampling length is set to c, the explosive quantity is set to Q, the sources are sequentially excited, the seismic records of all two-component receivers are obtained, and then other signals except for P-waves are removed by using a conventional method to obtain 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 receivers, (x1, y1; x2, y2; … x M+N ,y M+N ) is the position parameter of the receivers, t is the time, and S is the number of sources. One source corresponds to N+M receivers, which are arranged in combination to form actual observation data. Specifically, the receiver sampling interval can be set to 0.001s, and the sampling length can be set to 1000.

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

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

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

[0091] S311. Determine the observation system parameters based on the earthquake early detection system, and calculate the forward propagation wavefield data of different earthquake sources based on the observation system parameters using the spatiotemporal high-order finite difference algorithm.

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

[0093] In the formula, L[·] is the higher-order finite difference operator of sound wave spatiotemporal space, v n (x,y) represents the velocity model after the nth iteration, where v0(x,y) is the initial input, which is a uniform model. Subsequently, 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 Forward propagating wavefield data of (x,y), S(x) S ,y S ,t) is the source wavelet, (x S ,y S () represents the location of the earthquake source, and t represents the time parameter;

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

[0095] In the steps S311 and S312, the observation system is set according to the arrangement of the earthquake early warning system in step S1, that is, the arrangement of the observation system is consistent with the arrangement of the actual earthquake early warning system, then the observation system parameters are substituted into the space-time high-order finite difference algorithm, the velocity model is numerically simulated, and the forward wave field data p n (x,y,t) is obtained. n (x,y,t) is obtained. n (x,y,t) is obtained. n (x,y,t) is obtained. The space-time high-order finite difference method has higher calculation stability and faster calculation efficiency compared with the traditional finite difference method. 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 early warning system, and the steps of calculating the forward wave field data of different sources by using the space-time high-order finite difference algorithm based on the observation system parameters include:

[0097] S3111, obtaining a velocity model;

[0098] S3112, setting forward parameters, the grid spacing is dx, dy, the source position is (x s ,y s ), using a Ricker wavelet for forward, the main frequency is f, arranging the observation system, the arrangement of the observation system is consistent with the earthquake early warning system, and using the space-time high-order finite difference algorithm for forward simulation;

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

[0100] S3114, solving the wave equation by using the space-time high-order finite difference algorithm, specifically as follows:

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

[0102]

[0103] In the formula, λ and μ are Lame parameters, ρ is density, λ+2μ=ρα 2 , μ=ρβ 2 , α and β are respectively the velocities of longitudinal and transverse waves, v x and v y are the velocity components of particle vibration, τ 12 , τ 111 , τ 111x and τ 111y are stress components of particle vibration, the partial derivative with respect to time t;

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

[0105]

[0106] where p represents a wavefield variable, and Δt is a time discretization step;

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

[0108]

[0109] where h φ represents a spatial discretization step in the φ direction, represents a spatial discretization step in the r direction;

[0110] (4) The stability condition of the difference format is as follows:

[0111]

[0112] where

[0113] (5) The numerical simulation of the entire wavefield p(x, y, t) is calculated using the following spatiotemporal high-order difference operator n

[0114]

[0115] where φ represents an orthogonal coordinate axis other than the principal axis of derivative r, r, φ ∈ (x, z), r ≠ φ, h r represents a spatial discretization step in the r direction, and N r represents half of 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 the reverse wavefield data based on the velocity model and the actual observation data includes:

[0118] S331, using formula (2) to calculate the reverse wavefield data d n (x, y, t) by taking the actual observation data as a boundary value condition.

[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] In the formula, d n (x,y,t) is based on v n Inverse wavefield data for (x,y),

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

[0122] (x1,y1;x2,y2;…x M+N ,y M+n ) represents the position parameter of the detector, N+M represents the total number of detectors, t represents time, and S represents the number of seismic sources.

[0123] In step S33 above, the expression for the objective equation is:

[0124]

[0125] In the formula, F[v n [x,y)] represents the objective function value after the nth iteration, v n (x,y) represents the velocity model after the nth iteration, p n (x,y,t) is based on v n Forward propagation wavefield data for (x,y), d n (x,y,t) is based on v n The inverse wavefield data of (x,y), where ||·|| is the L2 norm, ζ is the regularization scale factor, W is the regularization weight operator, and v apr This is the approximate solution to the expectation.

[0126] In step S34 above, the calculated forward propagation wavefield data p is used. n (x,y,t) and the inverse wavefield data d n Calculate the gradient matrix of (x, y, t) with respect to the wave velocity. The expression for the gradient matrix is:

[0127]

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

[0129] In step S35 above, based on the gradient matrix g[vn (x,y)] to calculate the Hessian matrix, the expression of which is:

[0130]

[0131] where k n is the update amount of the gradient, k n = g[v n+1 (x,y)] - g[v n (x,y)], q n is the update amount with respect to the velocity model v n (x,y), q n = v n+1 (x,y) - v n (x,y), T represents the matrix transpose, is the initial positive definite Hessian matrix of the n-th iteration, is the approximate Hessian matrix after the j-th recursive iteration based on , and represents the first-order partial derivative with respect to v n (x,y), i.e., the second-order partial derivative of F[v n (x,y)] with respect to v n (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 length includes:

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

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

[0135] where a n is the iteration step length, v n (x,y) is the n-th iteration velocity model, v n+1 (x,y) is the velocity model after the n+1-th iteration, p n is the iteration direction, and the expression of p n is:

[0136]

[0137] where γ is a damping coefficient, and I is a unit matrix.

[0138] In step S361, the initial velocity model v0(x, y) is substituted into formula (6) to update, where n is 0, to obtain the updated velocity model v1(x, y), which is iteratively updated, v n (x, y) is obtained by iteratively updating v n-1 (x, y) is obtained by iteratively updating v n+1 (x, y) is obtained by iteratively updating v n (x, y) is obtained by iteratively updating v

[0139] After the initial velocity model is updated in step S37, a new velocity model is obtained. After the initial velocity model is replaced by the new velocity model, steps S31 to S36 are repeatedly executed until a preset convergence criterion is met to obtain a final velocity model. The preset convergence criterion can be, for example, that the objective function value F[v(x, y)] is less than or equal to a preset threshold, which can be set to a minimum value ε0, ε0=0.000001. For example, after the initial velocity model v0(x, y) is updated to obtain a new velocity model v1(x, y), the new velocity model v1(x, y) is used as the initial velocity model, and the objective function value F[v(x, y)] is obtained through steps S31 and S32. It is determined whether the objective function value F[v(x, y)] is less than or equal to ε0. If yes, the velocity model v1(x, y) at this time is used as the final velocity model of the full waveform inversion and output. Otherwise, steps S33 to S36 are continuously executed to update the velocity model v1(x, y) to obtain a new velocity model v2(x, y), and the iteration is repeatedly performed until the objective function value F[v(x, y)] is less than or equal to ε0.

[0140] Embodiment 2

[0141] A tunnel seismic advanced detection cave system, comprising:

[0142] A detection system erection unit is configured to arrange multiple seismic sources at equal intervals on a side wall of a tunnel and multiple geophones at equal intervals on the side wall of the tunnel, the multiple geophones being located between the seismic sources and a tunnel face, and multiple geophones being arranged at equal intervals on the tunnel face, each geophone and each seismic source being located at the same height.

[0143] An observation data acquisition unit is configured to connect each seismic source and all geophones through equipment to form a seismic advanced detection system, sequentially excite the seismic sources, and acquire actual observation data through the seismic advanced detection system.

[0144] A velocity model updating unit is configured to establish an initial velocity model, update the velocity model by using a Newton method, and obtain a final velocity model.

[0145] The step of updating the velocity model by using the Newton method to obtain the final velocity model comprises:

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

[0147] A backward wavefield data acquisition module is configured to obtain backward wavefield data according to the velocity model and actual observation data;

[0148] A target function construction module is configured to determine a target equation of a residual according to the forward wavefield data and the backward wavefield data;

[0149] A gradient matrix calculation module is configured to calculate a gradient matrix of wave velocity by using the forward wavefield data and the backward wavefield data;

[0150] A Hessian matrix calculation module is configured to calculate a Hessian matrix based on the gradient matrix;

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

[0152] An iteration judgment module is configured to repeatedly execute the forward wavefield data acquisition module to the model updating module until a preset convergence criterion is met, and obtain the final velocity model.

[0153] Embodiment 3

[0154] A storage medium storing a program, the program being executed by a processor to implement a tunnel seismic advanced detection cave method according to Embodiment 1, and the method comprises the following steps of:

[0155] S1. A plurality of seismic sources are arranged at equal intervals on one side wall of a tunnel, and a plurality of geophones are arranged at equal intervals between the seismic sources and a working face. A plurality of geophones are arranged at equal intervals on the working face, and each geophone and each seismic source are located at the same height.

[0156] S2. Each seismic source and all geophones are connected by a device to form a seismic advanced detection system, the seismic sources are sequentially excited, and actual observation data are acquired by the seismic advanced detection system.

[0157] S3. An initial velocity model is established, the velocity model is updated by using the Newton method, and a final velocity model is obtained.

[0158] The step of updating the velocity model by using the Newton method to obtain the final velocity model comprises the following steps of:

[0159] S31. Synthetic observation data and forward wavefield data are determined according to the velocity model.

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

[0161] S33, determining a target equation of residual according to the forward-propagation wave field data and the back-propagation wave field data;

[0162] S34, calculating a gradient matrix of wave velocity by using the forward-propagation wave field data and the back-propagation 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 a preset iteration step;

[0165] S37, repeatedly performing 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 Embodiment 1, which will not be repeated here.

[0167] In this embodiment, the storage medium can be a disk, an optical disk, a computer memory, a read-only memory, a random access memory, a U disk, a mobile hard disk, etc.

[0168] Embodiment 4

[0169] A computing device includes a processor and a memory for storing a processor-executable program, and the processor implements a tunnel seismic advanced detection cave method according to Embodiment 1 when executing the program stored in the memory, as follows:

[0170] S1, a plurality of sources and a plurality of receivers are arranged at equal intervals on a side wall of a tunnel, the plurality of receivers are located between the sources and a working face, a plurality of receivers are arranged at equal intervals on the working face, and each receiver and each source are located at the same height;

[0171] S2, each source and all receivers are connected by a device to form a seismic advanced detection system, the sources are sequentially excited, and actual observation data are obtained by the seismic advanced detection system;

[0172] S3, an initial velocity model is established, the velocity model is updated by using a Newton method, and a final velocity model is obtained;

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

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

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

[0176] S33, determining a target equation of residual according to the forward-propagation wave field data and the back-propagation wave field data;

[0177] S34, calculating a gradient matrix of wave velocity by using the forward-propagation wave field data and the back-propagation 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 a preset iteration step length;

[0180] S37, repeatedly performing steps S31 to S36 until a preset convergence criterion is met, and obtaining a final velocity model.

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

[0182] In this embodiment, the computing device can be a desktop computer, a notebook computer, a PDA handheld terminal, a tablet computer or the like terminal device.

[0183] The above embodiments are the preferred embodiments of the present application, and cannot limit the present application, and any changes or other equivalent replacement manners without departing from the technical scheme of the present application are all included in the protection scope of the present application.

Claims

1. A method of tunnel seismic advanced detection of cavities, characterized in that, The method comprises the following steps: S1, arranging multiple seismic sources and multiple geophones at equal intervals on a sidewall of a tunnel, the multiple geophones being located between the seismic sources and a tunnel face, multiple geophones being arranged at equal intervals on the tunnel face, each geophone and each seismic source being located at the same height; S2, connecting each seismic source and all geophones through a device to form a seismic advanced detection system, sequentially exciting the seismic sources, and obtaining actual observation data through the seismic advanced detection system; S3, establishing an initial velocity model, updating the velocity model by using a Newton method, and obtaining a final velocity model; The step of updating the velocity model by using the Newton method and obtaining the final velocity model comprises: S31, determining synthetic observation data and forward wave field data according to the velocity model; S32, obtaining backward wave field data according to the velocity model and actual observation data; S33, determining a target equation of a residual error according to the forward wave field data and the backward wave field data; In the step of determining the target equation of the residual error according to the forward wave field data and the backward wave field data, the expression of the target equation is: where 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 the forward wavefield data based on v n (x,y), d n (x,y,t) is the backward wavefield data based on v n (x,y), ||·|| is the two-norm, ζ is the regularization scaling factor, W is the regularization weight operator, and v apr is the desired approximate solution. S34, calculating a gradient matrix of wave velocity by using the forward wave field data and the backward wave field data; S35, calculating a Hessian matrix based on the gradient matrix; In the step of calculating the Hessian matrix based on the gradient matrix, the expression of the Hessian matrix is: where k n is the update amount of the gradient, k n = g[v n+1 (x, y)] - g[v n (x, y)], q n is the update amount with respect to the velocity model v n (x, y), q n = v n+1 (x, y) - v n (x, y), and T represents matrix transposition, is the initial positive definite Hessian matrix of the nth iteration, is the approximate Hessian matrix after the jth recursive iteration is performed on the basis of , and represents the first-order partial derivative with respect to v n (x, y), and I is the unit matrix; S36, updating the velocity model according to the gradient matrix, the Hessian matrix, and a preset iteration step length; S37, repeatedly performing steps S31 to S36 until a preset convergence criterion is met, and obtaining the final velocity model.

2. The method of tunnel seismic advanced detection of cavities according to claim 1, characterized in that, The step of determining the synthetic observation data and the forward wave field data according to the velocity model comprises: determining an observation system parameter according to the seismic advanced detection system, calculating the forward wave field data by using a time-space high-order finite difference algorithm based on the observation system parameter, and determining the synthetic observation data based on the forward wave field data. L[v n (x,y)]p n (x,y,t)=S(x S ,y S ,t) (1); where L[·] is the acoustic wave time-space high-order finite difference operator, v n (x,y) is the velocity model after the nth iteration, p n (x,y,t) is the forward wavefield data based on v n (x,y), S(x S ,y S ,t) is the source wavelet, (x S ,y S ) is the source position, and t is the time parameter. The step of obtaining the backward wave field data according to the velocity model and the actual observation data comprises:

3. The method of claim 1, wherein, taking the actual observation data as a boundary value condition, and calculating the backward wave field data by using formula (2): In the step of calculating the gradient matrix of wave velocity by using the forward wave field data and the backward wave field data, the expression of the gradient matrix is: L[v n (x,y)]d n (x,y,t)=D(S;x1,y1;x2,y2;…x M+N ,y M+N ;t) (2); where v n (x, y) is the velocity model after the nth iteration, d n (x, y, t) is the reverse wavefield data based on v n (x, y), 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 position parameter of the receiver, N+M represents the total number of receivers, and t is the time.

4. The method of claim 1, wherein, The step of updating the velocity model according to the gradient matrix, the Hessian matrix, and the preset iteration step length comprises: wherein F[v n (x,y)] is the first order partial derivative of F[v n (x,y) with respect to v n (x,y) is the velocity model after the nth iteration.

5. The method of claim 1, wherein, updating the velocity model by using formula (6), and the expression of formula (6) is: wherein γ is a damping coefficient, and I is a unit matrix. v n+1 (x,y) = v n (x,y) + a n p n (6); where a n is the iteration step, v n is the velocity model after the nth iteration, v n+1 is the velocity model after the (n+1)th iteration, p n is the iteration direction, p n is given by The method comprises the following steps:

6. A system for tunnel seismic advanced detection of cavities, characterized in that, The detection system setting unit is configured to arrange multiple seismic sources and multiple geophones at equal intervals on a sidewall of a tunnel, the multiple geophones being located between the seismic sources and a tunnel face, multiple geophones being arranged at equal intervals on the tunnel face, each geophone and each seismic source being located at the same height; ​ An observation data acquisition unit is configured to connect each of the seismic sources and all the geophones through devices to form a seismic advanced detection system, sequentially excite the seismic sources, and acquire actual observation data through the seismic advanced detection system; A velocity model updating unit is configured to establish an initial velocity model, update the velocity model by using a Newton method, and obtain a final velocity model; The updating of the velocity model by using the Newton method and the obtaining of the final velocity model include: A forward wave field data acquisition module is configured to determine synthetic observation data and forward wave field data according to the velocity model; A backward wave field data acquisition module is configured to obtain backward wave field data according to the velocity model and the actual observation data; A target function construction module is configured to determine a target equation of a residual according to the forward wave field data and the backward wave field data; The expression of the target equation is: where 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 the forward wavefield data based on v n (x,y), d n (x,y,t) is the backward wavefield data based on v n (x,y), and ||·|| is the two-norm, ζ is a regularization scaling factor, and W is a regularization weight operator, v apr is the desired approximate solution. A gradient matrix calculation module is configured to calculate a gradient matrix of wave velocity by using the forward wave field data and the backward wave field data; A Hessian matrix calculation module is configured to calculate a Hessian matrix based on the gradient matrix; The expression of the Hessian matrix is: where k n is the update amount of the gradient, k n = g[v n+1 (x, y)] - g[v n (x, y)], q n is the update amount with respect to the velocity model v n (x, y), q n = v n+1 (x, y) - v n (x, y), and T represents matrix transposition, is the initial positive definite Hessian matrix of the nth iteration, is the approximate Hessian matrix after the jth recursive iteration based on , and represents the first-order partial derivative with respect to v n (x, y), and I is the unit matrix; A model updating module is configured to update the velocity model according to the gradient matrix, the Hessian matrix, and a preset iteration step length; An iteration judgment module is configured to repeatedly execute the forward wave field data acquisition module to the model updating module until a preset convergence criterion is met, and obtain a final velocity model.

7. A storage medium, characterized by A program is stored, and the program is executed by a processor to implement the tunnel seismic advanced detection cave method of any one of claims 1-5.

8. A computing device, comprising: A processor and a memory for storing a processor executable program are included, and the processor executes the program stored in the memory to implement the tunnel seismic advanced detection cave method of any one of claims 1-5.

Images