A tunnel engineering detection method and device based on viscoelastic wave field simulation

Abstract

The present invention provides a tunnel engineering detection method and device based on viscoelastic wave field simulation in the field of tunnel engineering geological detection technology. The method includes the following steps: completing the source and sensor settings according to the specific conditions of the construction environment in the tunnel; after each source is excited, the sensor automatically collects data and stores the collected data, and performs denoising on the collected data; designing a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation; using the stored data to perform full waveform inversion and reverse time migration imaging. It can better fit the real physical process of wave propagation, thereby improving the inversion imaging accuracy in geological detection, and achieves the purpose of improving the accuracy of seismic wave inversion imaging by using a forward modeling method that uses viscoelastic waves as a physical model of wave propagation in a medium and a full waveform inversion method based on the viscoelastic wave forward modeling method.

Description

Technical Field

[0001] The present invention relates to the technical field of tunnel engineering geological detection, and in particular to a tunnel engineering geological detection method and device based on viscoelastic wave field simulation. Background Art

[0002] The main difficulty in rapid tunnel excavation is to understand the geological conditions and rock mechanics parameters ahead. The geological interface of the tunnel axis may cause serious problems during construction and excavation. In tunnel construction, if you encounter bad geological structures, it is easy to cause disasters such as water and mud inrush, especially when these disasters occur at the same time, the problem will be more serious, seriously threatening the lives of construction workers and causing damage to construction equipment, resulting in huge property losses.

[0003] In order to minimize the occurrence of accidents during tunnel construction and ensure the safety of construction, tunnel advance geological prediction is an effective way to solve this problem. Among them, the geophysical advance prediction method in the tunnel advance geological prediction method is more and more widely used in tunnels. The geophysical advance prediction method is to use geophysical exploration methods in the tunnel to detect the tunnel construction face and the surrounding adjacent areas, and to find out the nature, location and scale of the unfavorable geological body according to the difference in physical properties between the surrounding rock and the unfavorable geological body. At present, the geophysical detection methods for tunnel advance geological prediction mainly include high-density electricity, ground penetrating radar, seismic wave method, infrared detection method and beam (bore-tunnenng elec-trical ahead monitoring) method. Among them, the seismic wave method has become one of the important methods for tunnel advance prediction due to its advantages of high interface imaging accuracy and long detection distance.

[0004] However, the traditional seismic wave method is based on the classical wave equation, but a lot of actual data evidence shows that the propagation of waves in actual media will experience waveform attenuation, and the classical non-attenuated wave equation cannot explain the phenomenon of waveform attenuation. Summary of the invention

[0005] The traditional seismic wave method is based on the classical wave equation, but a lot of actual data evidence shows that the propagation of waves in actual media will experience waveform attenuation, and the classical non-attenuated wave equation cannot explain the phenomenon of waveform attenuation. The present invention provides a tunnel engineering geological detection method and device based on viscoelastic wave field simulation.

[0006] In a first aspect, the technical solution of the present invention provides a tunnel engineering detection method based on viscoelastic wave field simulation, comprising the following steps:

[0007] Complete the setting of seismic sources and sensors according to the specific conditions of the construction environment in the tunnel;

[0008] After each source is excited, the sensor automatically collects and stores data, and performs denoising on the collected data;

[0009] Design a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation;

[0010] Full waveform inversion and reverse time migration imaging are performed using stored data.

[0011] In this application, a viscoelastic wave forward modeling method is designed to solve the viscoelastic seismic wave equation and resolve the bottleneck problem of non-local wave equation calculation.

[0012] As a preferred embodiment of the technical solution of the present invention, the step of denoising the collected data uses a wavelet threshold denoising method to perform denoising, specifically comprising:

[0013] Use the wavelet decomposition formula to perform wavelet decomposition;

[0014] The threshold λ is calculated using the formula: N is the number of decomposition levels;

[0015] The soft threshold denoising formula is used for denoising; the formula is as follows:

[0016]

[0017] Among them, w is the wavelet coefficient;

[0018] Design a wavelet reconstruction formula and use the wavelet reconstruction formula to obtain the denoised signal;

[0019] Wavelet reconstruction formula:

[0020]

[0021]

[0022] As a preferred embodiment of the technical solution of the present invention, the steps of designing a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation include:

[0023] The viscoelastic seismic wave equations are constructed by integrating the detected area, the locations of the seismic source and the sensors into the numerical model;

[0024] The constructed viscoelastic wave equations are transformed into an equivalent problem that is local in time and non-local in space using the short memory splitting format;

[0025] The transformed problem is numerically solved using the spline method;

[0026] The boundary is treated using the auxiliary differential equation perfectly matched layer method to obtain the absorbing boundary region equation set.

[0027] As a preferred embodiment of the technical solution of the present invention, the viscoelastic seismic wave equations are:

[0028]

[0029]

[0030] And the time fractional stress-strain relationship:

[0031]

[0032] in, u1 and u3 are deformation components, ε 11 , ε 13 , ε 33 is the strain component, σ 11 , σ 13 , σ 33 is the stress component, f1 and f3 are source terms;

[0033] By transforming the flat singular kernel t -2γ Transformed into a relatively localized exotic nucleus y 4γ-1 , that is,

[0034]

[0035] Convert the equation into an equivalent problem that is local in time and non-local in space.

[0036] As a preferred embodiment of the technical solution of the present invention, the step of numerically solving the transformed problem using the spline method includes:

[0037] The spline simulation formula is:

[0038]

[0039] make:

[0040] get:

[0041] Additional boundary conditions:

[0042]

[0043] Obtain γ ν (y j );

[0044] In the actual calculation, f(x, y) is the actual wave field value, s(x, y) is the simulated wave field value, and η ν,βis the coefficient of the conventional two-dimensional calculation, γ ν (y j ) is the coefficient after conversion to a one-dimensional problem, B ν (x) is the x-direction spline basis function, B β (y) is the spline basis function in the y direction.

[0045] As a preferred embodiment of the technical solution of the present invention, the step of using the auxiliary differential equation perfectly matched layer method to process the boundary to obtain the absorbing boundary region equation group includes:

[0046] Set the relevant calculation formulas of the stretching coordinate system;

[0047] According to the differential theorem of Fourier transform, the auxiliary differential equation is obtained;

[0048] The auxiliary differential equations are solved using an integral solver, and finally the absorbing boundary region equations are obtained.

[0049] As a preferred embodiment of the technical solution of the present invention, in the step of performing full waveform inversion and reverse time migration imaging using stored data, the step of performing full waveform inversion includes:

[0050] Calculate forward wavefield data using the initial model;

[0051] Record simulated observation data and actual observation data;

[0052] Calculate the residual data by subtracting the actual observation data from the simulated observation data;

[0053] Using residual data to obtain reverse propagation wavefield data;

[0054] The forward wave field data and the reverse wave field data are cross-correlated to obtain the gradient matrix of the model error with respect to the wave velocity;

[0055] The gradient matrix is ​​imported into the model updating program, and the velocity model is updated using the quasi-Newton method. After the velocity model is updated, the error function is calculated until the calculation result of the error function is less than or equal to the allowable error, and the inversion velocity model is obtained.

[0056] As a preferred embodiment of the technical solution of the present invention, the step of performing cross-correlation calculation on the forward wavefield data and the reverse wavefield data to obtain a gradient matrix of the model error with respect to the wave velocity includes:

[0057] The forward wave field data d 正 and reverse propagation wavefield data d 逆 At each time point, a corresponding multiplication operation is performed in two-dimensional space, and then each time node is added, that is, a cross-correlation operation is performed. The result of the cross-correlation calculation is the gradient matrix J of the full waveform inversion:

[0058]

[0059] Among them, d 正 is the forward wave field data, d 逆 It is the reverse propagation wavefield data.

[0060] As a preferred embodiment of the technical solution of the present invention, in the step of performing full waveform inversion and reverse time migration imaging using stored data, the step of performing reverse time migration imaging includes:

[0061] The inversion model is used to perform reverse time migration imaging, and the wave field forward propagation is performed in the inversion velocity model;

[0062] The collected data is reversely transmitted in the wave field;

[0063] The forward wavefield data and the reverse wavefield data are cross-correlated to obtain the reverse time migration imaging result.

[0064] The technical solution of the present invention provides a tunnel engineering detection device based on viscoelastic wave field simulation, including a data acquisition and preprocessing module, a forward modeling module, and an inversion imaging module;

[0065] The seismic source module is set in the construction tunnel according to the specific conditions of the construction environment in the tunnel, and is used to output simulated seismic signals after excitation;

[0066] The data acquisition and preprocessing module is used to automatically acquire data from the simulated seismic signal after each source excitation, collect and store the data, and perform denoising on the collected data;

[0067] Forward modeling module, used to design viscoelastic wave forward modeling method to solve viscoelastic seismic wave equations;

[0068] Inversion imaging module for full waveform inversion and reverse time migration imaging using stored data.

[0069] Based on the viscoelastic wave equation, compared with the classical wave equation, when studying problems such as wave propagation in viscoelastic and anelastic media, it can better fit the real physical process of wave propagation, thereby improving the inversion imaging accuracy in geological exploration. The purpose of improving the accuracy of seismic wave inversion imaging is achieved by using the forward modeling method of viscoelastic waves as the physical model of wave propagation in the medium and the full waveform inversion method based on the viscoelastic wave forward modeling method.

[0070] It can be seen from the above technical solutions that the present invention has the following advantages: 1. In terms of the wave equation model, the viscoelastic wave model term used in the present invention can be more in line with the wave propagation process in the real rock and soil medium for the classical wave equation model, so it is more reasonable in the physical model construction; 2. In terms of solving the fractional-order wave equation, the short memory operator splitting (SMOS) format designed by the present invention shows in actual calculations that when the fractional order is very small (≤0.06), SMOS can use 32 layers of memory variables to achieve the accuracy of the original method of 500 layers of memory variables (L^2 error is about 2%), and reduces the storage and calculation time, greatly reducing the calculation cost; 3. In terms of inversion imaging, since the physical model used is more practical, the imaging resolution will be higher than using the classical wave equation as the forward physical model.

[0071] In addition, the invention has a reliable design principle, a simple structure and a very broad application prospect.

[0072] It can be seen that compared with the prior art, the present invention has outstanding substantive features and significant progress, and the beneficial effects of its implementation are also obvious. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0074] Figure 1 is a schematic flow chart of a method according to an embodiment of the present invention.

[0075] Figure 2 It is a schematic flow chart of performing reverse time migration imaging in a method according to an embodiment of the present invention.

[0076] Figure 3 It is a schematic flow chart of performing full waveform inversion in a method according to an embodiment of the present invention.

[0077] Figure 4 is a schematic block diagram of an apparatus according to an embodiment of the present invention.

[0078] Figure 5 Schematic diagram of a tunnel example according to an embodiment of the present invention. DETAILED DESCRIPTION

[0079] The traditional seismic wave method is based on the classical wave equation, but a lot of actual data evidence shows that the propagation of waves in actual media will experience waveform attenuation, and the classical wave equation without attenuation cannot explain this phenomenon. The present application provides a viscoelastic wave equation, which can better fit the real physical process of wave propagation when studying problems such as wave propagation in viscoelastic and anelastic media, thereby improving the inversion imaging accuracy in geological exploration. In order to enable those skilled in the art to better understand the technical solutions in the present invention, 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 only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of the present invention.

[0080] like Figure 1 As shown, an embodiment of the present invention provides a tunnel engineering detection method based on viscoelastic wave field simulation, comprising the following steps:

[0081] Step 1: Complete the setting of the seismic source and sensors according to the specific conditions of the construction environment in the tunnel;

[0082] Step 2: After each source is excited, the sensor automatically collects and stores data, and performs denoising on the collected data;

[0083] Step 3: Design a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation;

[0084] Step 4: Perform full waveform inversion and reverse time migration imaging using the stored data.

[0085] In this application, a viscoelastic wave forward modeling method is designed to solve the viscoelastic seismic wave equation and resolve the bottleneck problem of non-local wave equation calculation.

[0086] In the embodiment of the present invention, the step of denoising the collected data uses a wavelet threshold denoising method to perform denoising, which specifically includes:

[0087] Step 21: Use the wavelet decomposition formula to perform wavelet decomposition; first approximate the original signal with a step signal: f j Decomposed into: f j =ω j-1 +f j-1

[0088] in:

[0089]

[0090]

[0091] Finally, we can get: f j =ω j-1 +ω j-2 +…+ω0+f0.

[0092] Step 22: Calculate the threshold λ using the formula: N is the number of decomposition levels;

[0093] Step 23: Use the soft threshold denoising formula to perform denoising; the formula is as follows:

[0094]

[0095] Among them, w is the wavelet coefficient;

[0096] Step 24: Design a wavelet reconstruction formula and use the wavelet reconstruction formula to obtain the denoised signal;

[0097] Wavelet reconstruction formula:

[0098]

[0099]

[0100] In the example of the present invention, the steps of designing a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation include:

[0101] Step 31: constructing a set of viscoelastic seismic wave equations by integrating the detected area, the locations of the seismic source and the sensors into the numerical model;

[0102] Step 32: Use the short memory splitting format to transform the constructed viscoelastic wave equations into an equivalent problem that is local in time and non-local in space;

[0103] The Short-Memory Operator Splitting (SMOS) format for the two-dimensional fractional wave equation is designed; the viscoelastic seismic wave equations are:

[0104]

[0105]

[0106] And the time fractional stress-strain relationship:

[0107]

[0108] in, u1 and u3 are deformation components, ε 11 , ε 13 , ε 33 is the strain component, σ 11 , σ 13 , σ 33 is the stress component, f1 and f3 are source terms;

[0109] By transforming the flat singular kernel t -2γ Transformed into a relatively localized exotic nucleus y 4γ-1 , that is,

[0110]

[0111] Convert the equation into an equivalent problem that is local in time and non-local in space.

[0112] Step 33: Use the spline method to numerically solve the transformed problem;

[0113] Use cubic B-spline basis functions:

[0114]

[0115] This formula can be obtained from the Cox-de Boor recursive formula:

[0116]

[0117]

[0118] For one-dimensional problems:

[0119] Among them, i=-1 and n+1 are boundary conditions; when n=1,2,…,n-1

[0120]

[0121] Add the following two boundary conditions:

[0122] The matrix calculation of the system of equations is in the form:

[0123]

[0124] Where A is:

[0125] C and F are:

[0126] C=[c -1 ,c0,c1,…,c n-1 ,c n ,c n+1 ] T

[0127] F=[f'(u0),f(u0),f(u1),…,f(u n-1 ),f(u n ),f'(u n+1 )] T

[0128] For two-dimensional problems, in order to speed up the calculation, the original two-dimensional spline simulation method is no longer used. Instead, the two-dimensional problem is converted into a combination of one-dimensional problems for solution.

[0129] The spline simulation formula is:

[0130]

[0131] make:

[0132] get:

[0133]

[0134] Additional boundary conditions:

[0135]

[0136] Obtain γ ν (y j ). All coefficients can be obtained by performing a second solution using similar steps. The entire coefficient calculation process only needs to be solved once each time. After the coefficients are solved, the approximate spatial derivatives are obtained by multiplying the coefficients with the derivatives of the basis functions.

[0137] Step 34: Use the auxiliary differential equation perfectly matched layer method to process the boundary to obtain the absorbing boundary region equation group, which specifically includes:

[0138] Step 341: Set the relevant calculation formula of the stretching coordinate system:

[0139]

[0140]

[0141]

[0142] in, represents the stretched coordinate system, and in addition:

[0143]

[0144] k i =1+(k max -1) m

[0145]

[0146] Where L is the thickness of the PML layer, N = 2, σ0 = -(N+1)v pmax log(R) / 2L,v pmax is the velocity of the pressure wave, R is the target theoretical reflection coefficient, which is selected as 0.1% here. max =πf0, where f0 is the main frequency of the earthquake source. max Usually between 1 and 20, we take m = 2. The first-order acoustic wave equation in the frequency domain under stretched coordinates is:

[0147]

[0148]

[0149]

[0150] In the frequency domain we have:

[0151]

[0152] make: get:

[0153] Step 342: Obtain an auxiliary differential equation according to the differential theorem of Fourier transform; the auxiliary differential equation is:

[0154] The acoustic wave equations for the entire PML region are:

[0155]

[0156] Step 343: Use an integral solver to solve the auxiliary differential equation, and finally obtain the absorbing boundary region equation group:

[0157]

[0158] In the embodiment of the present invention, in the step of using the stored data to perform full waveform inversion and reverse time migration imaging, as Figure 2 As shown, the steps of full waveform inversion include:

[0159] Step 41: Calculate the forward wave field data using the initial model, record the simulated observation data and the actual observation data; and calculate the residual data by subtracting the actual observation data from the simulated observation data;

[0160] Calculate the forward wavefield data d using the initial model 正 , record the simulated observation data ob 模拟 , the actual observed data is denoted as ob实际 , use the formula to calculate the residual data: ob 残差 =ob 实际 -ob 模拟 ;

[0161] Step 42: Using the residual data ob 残差 Get reverse propagation wavefield data d 逆 ;

[0162] Step 43: performing cross-correlation calculation on the forward wavefield data and the reverse wavefield data to obtain a gradient matrix of the model error with respect to the wave velocity;

[0163] The forward wave field data d 正 and reverse propagation wavefield data d 逆 At each time point, a corresponding multiplication operation is performed in two-dimensional space, and then each time node is added, that is, a cross-correlation operation is performed. The result of the cross-correlation calculation is the gradient matrix J of the full waveform inversion:

[0164]

[0165] Among them, d 正 is the forward wave field data, d 逆 It is the reverse propagation wavefield data.

[0166] Step 44: importing the gradient matrix into the model updating program, updating the velocity model using the quasi-Newton method, calculating the error function after updating the velocity model, until the calculation result of the error function is less than or equal to the allowable error, and obtaining the inversion velocity model; specifically including:

[0167] Import the gradient matrix into the model update program and use the quasi-Newton method to calculate the update matrix dm of the velocity model. First, calculate the search direction b: Then calculate the speed update matrix dm: λ = arg minf(m+λb)

[0168] dm=λb

[0169] Where D is the approximate Hessian matrix, and its update method is: dJ = J (t+1) -J (t)

[0170]

[0171] The updated velocity model matrix m1 is obtained by adding the original velocity model matrix and the updated matrix. The model error is calculated using the least square method. The updated model parameters are imported into the forward modeling program. The simulated seismic data are recorded. The second norm of the difference between the seismic data of the updated model and the seismic data before the update is calculated, which can be expressed as: Φ = ||f(m1)-f(m0)||2;

[0172] Where Φ is the error of the calculated velocity model, and the allowable error ε of the model is set to 0.01. If Φ>ε, the next inversion calculation is performed, and the updated velocity model is subjected to the above-mentioned gradient calculation process and model update process again until Φ≤ε, and the inversion velocity model m is obtained.

[0173] In the embodiment of the present invention, in the step of using the stored data to perform full waveform inversion and reverse time migration imaging, as follows Figure 3 As shown, the steps of performing reverse time migration imaging include:

[0174] Step 45: Perform reverse time migration imaging using the inversion model, and perform wave field forward propagation in the inversion velocity model;

[0175] Step 46: Reverse transmission of the collected data;

[0176] Step 47: Perform cross-correlation calculation on the forward wavefield data and the reverse wavefield data to obtain a reverse time migration imaging result.

[0177] like Figure 4 As shown, an embodiment of the present invention provides a tunnel engineering detection device based on viscoelastic wave field simulation, including a seismic source module, a data acquisition and preprocessing module, a forward modeling module, and an inversion imaging module;

[0178] The seismic source module is set in the construction tunnel according to the specific conditions of the construction environment in the tunnel, and is used to output simulated seismic signals after excitation;

[0179] The data acquisition and preprocessing module is used to automatically acquire data from the simulated seismic signal after each source excitation, collect and store the data, and perform denoising on the collected data;

[0180] Forward modeling module, used to design viscoelastic wave forward modeling method to solve viscoelastic seismic wave equations;

[0181] The inversion imaging module is used to perform full waveform inversion and reverse time migration imaging using the stored data. Figure 5 shown.

[0182] The first step is to generate active source seismic wave signals through the source module. The source can be a hammer source or a pneumatic source, etc. The source excitation frequency is about 200 to 400 Hz. In this example, two hammer sources are arranged on each side wall, and the interval between the sources is 2m.

[0183] The second step is to receive, transmit and denoise the seismic signals through the data acquisition and preprocessing module. In this example, the sensors are arranged about 5m away from the tunnel face on the front side wall of the tunnel, and holes are drilled every 2m with an electric drill. Every 2m, a hole is drilled with a diameter of about 50mm and a depth of about 1.60m. The total length is 10m, and there are 5 sensors on each side wall; the sensors can also be detectors here.

[0184] After each source is excited, the data acquisition module automatically collects data once. Each signal reception time is 0.1s, the sampling frequency is 8000 times, and each collected data is automatically stored until all sources are excited. After the data is stored, the collected data is denoised. In the preprocessing module, the wavelet threshold denoising method is used for denoising. The wavelet decomposition formula in the above text is used for wavelet decomposition. Four wavelet functions, Haar, Daubechies, Coiflets, and Symlets, are selected to perform 32-layer wavelet decomposition of the seismic data, and then the threshold is calculated. The calculation formula is:

[0185]

[0186] Where N is the number of decomposition layers, and λ is calculated to be 1.735.

[0187] Then use the soft threshold denoising formula for denoising:

[0188]

[0189] Among them, w is the wavelet coefficient;

[0190] Finally, the denoised signal is obtained using the wavelet reconstruction formula mentioned above.

[0191] The third step is to construct a forward model through a forward model construction module. In this embodiment, the detected area and the locations of the source and sensor are constructed in the numerical model. The source is set to the Ricker wavelet, and the formula is:

[0192] Where fp is the main frequency of the earthquake source, which is 300Hz. According to the sampling frequency, the time interval dt is:

[0193]

[0194] The constructed viscoelastic wave equations are solved using the short-memory operator splitting (SMOS) format. During the numerical simulation, the spline method is used to calculate the spatial derivative, the difference method is used to calculate the time derivative, and the auxiliary differential equation perfectly matched layer (ADE-Pml) is used for the absorbing boundary.

[0195] The fourth step is to perform inversion calculation. First, give the initial velocity model m0; import the initial velocity model matrix, perform a forward simulation in the initial velocity model, the simulation duration is 0.1s, the simulation time nodes are 0s to 0.1s, record the simulated wave field data and simulated seismic records, and then subtract the simulated seismic records from the previously stored sensor data to obtain residual data; use the residual data to perform reverse propagation wave field simulation in the wave field, the simulation duration is 0.1s, the simulation time nodes are 0.1s to 0s, and also record the simulated wave field data of the reverse propagation wave field; perform corresponding multiplication operations on the obtained forward wave field data and reverse wave field data in two-dimensional space at each time point, and then perform addition operations at each time node, that is, perform cross-correlation operations:

[0196]

[0197] The result of the cross-correlation calculation is the gradient matrix J of the full waveform inversion (the gradient of the objective function to the velocity). The gradient matrix is ​​imported into the model update program, and the update matrix dm of the velocity model is calculated using the quasi-Newton method. First, the search direction b is calculated:

[0198]

[0199] Then calculate the velocity update matrix dm:

[0200] λ=arg minf(m+λb)

[0201] dm=λb

[0202] Where D is the approximate Hessian matrix, and its update method is: dJ = J (t+1) -J (t)

[0203]

[0204] The updated velocity model matrix m1 is obtained by adding the original velocity model matrix and the updated matrix. The model error is calculated using the least squares method. The updated model parameters are imported into the forward modeling program. The simulated seismic data are recorded and the second norm of the difference between the seismic data of the updated model and the seismic data before the update is calculated. It can be expressed as:

[0205] Φ=||f(m1)-f(m0)||2

[0206] Where Φ is the error of the calculated velocity model, and the allowable error of the model is set to 0.01. If the calculated model error is greater than the allowable error, the next inversion calculation is performed, and the updated velocity model is subjected to the above-mentioned gradient calculation process and model update process again until the model error is less than the allowable error. The inversion velocity model m is obtained, and the wave field is forward propagated in the velocity model m, and then the collected data is reversely propagated. Finally, the forward propagation wave field data and the reverse propagation wave field data are cross-correlated to obtain the reverse time migration imaging result.

[0207] Although the present invention has been described in detail by referring to the accompanying drawings and in combination with the preferred embodiments, the present invention is not limited thereto. Without departing from the spirit and essence of the present invention, a person of ordinary skill in the art may make various equivalent modifications or substitutions to the embodiments of the present invention, and these modifications or substitutions shall be within the scope of the present invention. Any person of ordinary skill in the art may easily think of changes or substitutions within the technical scope disclosed by the present invention, and these shall be within the scope of protection of the present invention. Therefore, the scope of protection of the present invention shall be subject to the scope of protection of the claims.

Claims

1. A tunnel engineering detection method based on viscoelastic wave field simulation, characterized in that: The steps include: Complete the setting of seismic sources and sensors according to the specific conditions of the construction environment in the tunnel; After each source is excited, the sensor automatically collects and stores data, and performs denoising on the collected data; Design a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation; Full waveform inversion and reverse time migration imaging using stored data; Viscoelastic seismic wave equations: And the time fractional stress-strain relationship: in, , and is the deformation component, is the strain component, is the stress component, is the source term; By placing the flat singular core Transformed into relatively localized exotic nuclei , that is, Convert the equation into an equivalent problem that is local in time and non-local in space.

2. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 1 is characterized in that: The step of denoising the collected data uses a wavelet threshold denoising method to denoise the data, specifically including: Use the wavelet decomposition formula to perform wavelet decomposition; Use the formula to calculate the threshold; Use the soft threshold denoising formula to perform denoising; Design a wavelet reconstruction formula and use it to get the denoised signal.

3. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 2 is characterized in that: The steps of designing a viscoelastic wave forward modeling method to solve the viscoelastic seismic wave equation include: The viscoelastic seismic wave equations are constructed in the numerical model by taking the detected area and the locations of the seismic source and the sensors; The constructed viscoelastic wave equations are transformed into an equivalent problem that is local in time and non-local in space using the short memory splitting format; The transformed problem is numerically solved using the spline method; The boundary is treated using the auxiliary differential equation perfectly matched layer method to obtain the absorbing boundary region equation set.

4. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 3 is characterized in that: The steps to numerically solve the transformed problem using the spline method include: The spline simulation formula is: make: get: Additional boundary conditions: Obtain ; Among them, in actual calculation is the actual wave field value, is the simulated wave field value, are the coefficients for conventional two-dimensional calculations, are the coefficients after conversion to a one-dimensional problem, for Directional spline basis functions, for Directional spline basis functions.

5. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 4 is characterized in that: The steps of using the auxiliary differential equation perfectly matched layer method to process the boundary and obtain the absorbing boundary region equations include: Set the relevant calculation formulas of the stretching coordinate system; According to the differential theorem of Fourier transform, the auxiliary differential equation is obtained; The auxiliary differential equations are solved using an integral solver, and finally the absorbing boundary region equations are obtained.

6. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 5 is characterized in that: In the step of performing full waveform inversion and reverse time migration imaging using stored data, the steps of performing full waveform inversion include: Calculate forward wavefield data using the initial model; Record simulated observation data and actual observation data; Calculate the residual data by subtracting the actual observation data from the simulated observation data; Using residual data to obtain reverse propagation wavefield data; The forward wave field data and the reverse wave field data are cross-correlated to obtain the gradient matrix of the model error with respect to the wave velocity; The gradient matrix is ​​imported into the model updating program, and the velocity model is updated using the quasi-Newton method. After the velocity model is updated, the error function is calculated until the calculation result of the error function is less than or equal to the allowable error, and the inversion velocity model is obtained.

7. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 6 is characterized in that: The steps of performing cross-correlation calculation on the forward wavefield data and the reverse wavefield data to obtain the gradient matrix of the model error with respect to the wave velocity include: The forward wave field data obtained and reverse propagation wavefield data At each time point, a corresponding multiplication operation is performed in two-dimensional space, and then each time node is added, that is, a cross-correlation operation is performed. The result of the cross-correlation calculation is the gradient matrix J of the full waveform inversion: in, It is the positive transmission wave field data. It is the reverse propagation wavefield data.

8. The tunnel engineering detection method based on viscoelastic wave field simulation according to claim 7 is characterized in that: In the step of performing full waveform inversion and reverse time migration imaging using stored data, the step of performing reverse time migration imaging includes: The inversion model is used to perform reverse time migration imaging, and the wave field forward propagation is performed in the inversion velocity model; The collected data is reversely transmitted in the wave field; The forward wavefield data and the reverse wavefield data are cross-correlated to obtain the reverse time migration imaging result.

9. A tunnel engineering detection device based on viscoelastic wave field simulation, characterized in that: Source module, data acquisition and preprocessing module, forward modeling module, inversion imaging module; The seismic source module is set in the construction tunnel according to the specific conditions of the construction environment in the tunnel, and is used to output simulated seismic signals after excitation; The data acquisition and preprocessing module is used to automatically acquire data from the simulated seismic signal after each source excitation, collect and store the data, and perform denoising on the collected data; Forward modeling module, used to design viscoelastic wave forward modeling method to solve viscoelastic seismic wave equations; An inversion imaging module for performing full waveform inversion and reverse time migration imaging using stored data; Viscoelastic seismic wave equations: And the time fractional stress-strain relationship: in, , and is the deformation component, is the strain component, is the stress component, is the source term; By placing the flat singular core Transformed into relatively localized exotic nuclei , that is, Convert the equation into an equivalent problem that is local in time and non-local in space.

Images