A tunnel advance detection imaging method and device

Abstract

The present invention belongs to the field of tunnel advance detection, and specifically, is a tunnel advance detection imaging method and device, which can effectively process tunnel advance detection seismic data in the harsh environment of tunnel construction. The source wave field is extended to obtain the source extension wave field; seismic data on both sides of the tunnel and the face are collected; the seismic data is used as the boundary condition to extend the detector wave field to obtain the detector extension wave field; imaging condition calculation is performed on the source extension wave field and the detector extension wave field. The calculation efficiency is improved, the space occupation is reduced, and the availability and real-time performance of tunnel advance detection are guaranteed.

Description

Technical Field

[0001] The present invention belongs to the field of tunnel seismic prediction (TSP), and in particular to a tunnel seismic prediction imaging method and device. Background Art

[0002] Tunnel advance detection technology is an indispensable technology used in tunnel construction. It uses the reflection of seismic waves to image the unexcavated area behind the tunnel face. The imaging results can be used to guide the construction party to formulate a suitable construction plan to avoid accidents.

[0003] Due to the limitation of computing resources and real-time requirements, the TSP systems commonly used by construction parties currently mostly use the depth migration imaging technology based on the Kirchhoff method, such as the method proposed by Zha Xinjie, Wang Wei, Gao Xing, etc. in the paper "Application of Quasi-VSP and Kirchhoff Migration Method in Tunnel Advance Forecasting" (Geophysical and Geochemical Exploration, 2016, 40(01): 214-219.). This method has fast computing speed and small memory usage, but has poor energy focusing for inclined interfaces and the imaging results have obvious "arcing" phenomenon. The reverse time migration technology proposed by Baysal E and Kosloff D in the paper Reverse-Time Migration (Geophysics, 1983, 48: 1514-1524) is a time migration technology that has been widely used as computer computing efficiency improves. This method can effectively solve the problems of the Kirchhoff method and improve the imaging effect.

[0004] CN206833004U discloses a three-channel wireless seismograph data acquisition system, which can collect elastic wave signals and generally upload the collected data to a server or a large workstation for processing. However, in actual tunnel construction, the construction environment is humid and the space is small, so high-performance computing equipment cannot work, and the computing program can only be run in a well-sealed small-volume embedded device. However, the computing rate and computing resources of embedded devices are limited, and the existing reverse time migration imaging method is not well applicable to the field of tunnel advance detection. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a tunnel advance detection imaging method and device. The present invention is implemented as follows:

[0006] A tunnel advance detection imaging method, the method comprising:

[0007] The source wave field is extended to obtain the source extended wave field;

[0008] Collect seismic data on both sides of the tunnel and on the face;

[0009] Taking the seismic data as the boundary condition to carry out the detector wave field extension, the detector extension wave field is obtained;

[0010] Perform imaging condition calculations on the source extended wavefield and the receiver extended wavefield.

[0011] Furthermore, the source wave field is extended to obtain the source extended wave field including:

[0012] The two-dimensional first-order velocity stress elastic wave equation is selected as the basic equation for extension. The two-dimensional first-order velocity stress elastic wave equation is:

[0013]

[0014] Where x is the direction parallel to the tunnel face, z is the direction perpendicular to the tunnel face, t is the running time, σ xx , σ yy , σ zz is the normal stress on the microelement, v x and v z are the velocity components of the microelement in the x and z directions respectively, σ zx , σ xy , σ yz is the shear stress on the microelement, the Lame coefficients λ and μ and the propagation speed of elastic waves in the medium, v p and v s And the medium density ρ satisfies the following relationship:

[0015]

[0016]

[0017] Select the main frequency f m =500Hz Ricker wavelet as the source when the source wave field is extended;

[0018] The high-order staggered grid finite difference method is used as the numerical solution method to solve the two-dimensional first-order velocity stress elastic wave equation to obtain the differential iterative equation of any point in the simulation area, where the differential formula is:

[0019]

[0020] Where u is the velocity or stress term to be determined, A0 is the differential point to be determined, ΔA is the grid step size of the staggered grid, 2N is the difference order, and a is m is the differential coefficient, and A is x or z.

[0021] Further, the differential coefficient a is solved using the following equation: m :

[0022]

[0023] Where 2N is the difference order, b1=1, b n Satisfy the following formula:

[0024]

[0025] The parameters satisfy the following formula:

[0026] c n =(cosθ) 2n-1 β n , d n =(sinθ) 2n-1 β n , θ=π / 8.

[0027] 4. The method according to claim 2 is characterized in that an absorption function is selected to be introduced into the Ricker wavelet as a boundary condition, and the absorption function expression is as follows:

[0028]

[0029] Where L is the total number of layers in the absorption region, k is the number of layers from the forward modeling region, and R = 0.0001;

[0030] The Ricker wavelet expression is:

[0031]

[0032] Furthermore, the differential iterative equations at any point include:

[0033] The differential iteration equation on the x component is:

[0034]

[0035] The differential iteration equation on the z component is:

[0036]

[0037] Where: d(x0,z0) is the absorption matrix, which is 0 in the simulation area, now represents the stress or velocity value at the point (x0,z0) at the current moment, next represents the stress or velocity value to be calculated at the point (x0,z0) at the next moment, and dif x () indicates to find the difference of the variable in the x direction, dif z () indicates taking the difference of the variable in the z direction.

[0038] Furthermore, performing detector wave field extension includes:

[0039] The two-dimensional first-order velocity stress elastic wave equation is selected as the basic equation for extension. The two-dimensional first-order velocity stress elastic wave equation is:

[0040]

[0041] Where x is the direction parallel to the tunnel face, z is the direction perpendicular to the tunnel face, t is the running time, σ xx , σ yy , σ zz is the normal stress on the microelement, v x and v z are the velocity components of the microelement in the x and z directions respectively, σ zx , σ xy , σ yz is the shear stress on the microelement, the Lame coefficients λ and μ and the propagation speed of elastic waves in the medium, v p and v s And the medium density ρ satisfies the following relationship:

[0042]

[0043]

[0044] Select the main frequency f m =500Hz Ricker wavelet as the source when the source wave field is extended;

[0045] The high-order staggered grid finite difference method is used as the numerical solution method to solve the two-dimensional first-order velocity stress elastic wave equation to obtain the differential iterative equation of any point in the simulation area, where the differential formula is:

[0046]

[0047] Where u is the velocity or stress term to be determined, A0 is the differential point to be determined, ΔA is the grid step size of the staggered grid, 2N is the difference order, and a is m is the differential coefficient, and A is x or z.

[0048] Furthermore, the differential coefficient:

[0049] Use the following equation to solve for the differential coefficient a m .

[0050]

[0051] Where 2N is the difference order, b1=1, b n Satisfy the following formula

[0052]

[0053] The parameters in the formula satisfy the following formula:

[0054] c n =(cosθ) 2n-1 β n , d n =(sinθ) 2n-1 β n .

[0055] Furthermore, the imaging condition calculation of the source extended wave field and the detector extended wave field includes:

[0056] Calculate travel time to filter direct waves and remove direct waves;

[0057] Using the cross-correlation imaging condition, the results of the source wavefield extension and the receiver wavefield extension are cross-correlated. The cross-correlation imaging condition is as follows:

[0058]

[0059] Use normalization operation: Use fast Fourier transform to transform the source extension wave field S(x,z,t) and the detector extension wave field R(x,z,t) into the frequency domain to obtain S(x,z,ω) and R(x,z,ω), and then calculate the conjugate S of the frequency domain source extension wave field S(x,z,ω) * (x,z,ω). The operation expression of the normalized cross-correlation imaging condition is as follows:

[0060]

[0061] Use the prewitt operator to extract the boundary of the imaging result image(x,z).

[0062] A tunnel advance detection imaging device, based on an embedded computer device of Stm32, comprising:

[0063] The multi-channel collector is connected to the network module based on the ksz8463 chip through the SPI interface, which synchronizes the multi-channel collector clock and transmits and receives data;

[0064] The switch module based on triode circuit is used to control the triggering of the explosion source;

[0065] The SD card storage module is used to store the wave field data with a large amount of data generated during the imaging calculation process and the seismic data received through the network module;

[0066] LCD touch screen module is used to interact with the user;

[0067] The source wave field extension module extends the source wave field to obtain the source extension wave field;

[0068] A detector wave field extension module performs detector wave field extension to obtain a detector extension wave field;

[0069] The cross-correlation imaging module performs imaging condition calculations on the extended wavefield of the source and the extended wavefield of the detector.

[0070] Furthermore, multiple collectors are evenly distributed on both sides of the tunnel and on the face. The collectors on both sides of the tunnel are symmetrically positioned, and one collector is placed at each of the three points of the face. The seismic sources are set at the two corners of the tunnel face.

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

[0072] 1. The reverse time migration method is used as the imaging technology for tunnel advance detection. Compared with the widely used Kirchhoff imaging technology, it has better imaging effect on large-angle interfaces and less "arcing" phenomenon.

[0073] 2. By using staggered grids, variable step size grids, optimized PML absorbing boundaries and time-space domain difference coefficients, the time complexity and space complexity of finite difference calculations are effectively reduced, making it more suitable for lightweight embedded systems and easy to use in harsh environments.

[0074] 3. Use normalized cross-correlation imaging conditions to reduce the ray noise and dispersion noise caused by the imaging process, so that the signal-to-noise ratio of the tunnel detection imaging result is higher and the imaging effect is better. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 It is a basic structural diagram of the device of the present invention;

[0076] Figure 2 It is an observation system arrangement suitable for tunnel advance detection;

[0077] Figure 3 It is a staggered grid diagram of the first-order velocity stress elastic wave equation;

[0078] Figure 4 It is the dispersion curve diagram of different order staggered grid finite difference method;

[0079] Figure 5 It is a schematic diagram of the PML absorption boundary;

[0080] Figure 6 This is a schematic diagram of the reverse time migration imaging principle;

[0081] Figure 7 It is a schematic diagram of the inclined layered tunnel model, where the black part is the excavated tunnel, the gray part is the unexcavated rock mass, and the dotted lines represent the medium stratification;

[0082] Figure 8 yes Figure 7Under the model, the collected seismic data are common shot point gathers, where (a) is the z component and (b) is the x component;

[0083] Fig. 9 yes Figure 7 Reverse time migration imaging results of the model;

[0084] Fig.10 yes Fig. 9 Binary image of the model’s boundary extraction results. DETAILED DESCRIPTION

[0085] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0086] See also Figure 1 As shown, the present invention provides a tunnel advance detection imaging device, which is an embedded computer device based on Stm32, and the device structure is as follows Figure 1 As shown, the main control chip is linked to the network module 1 based on the ksz8463 chip through the SPI interface, and synchronizes the collector clock through it and performs data transmission and reception; the switch module 2 based on the triode circuit is used to control the triggering of the explosion source; the SD card storage module 3 is used to store the wave field data with a large amount of data generated during the imaging calculation process and the seismic data received through the network module 1; the LCD touch screen module 4 is used to interact with the user; the power supply module 5 is used for voltage conversion and power supply to the main control chip and other modules. The collected data will be processed by the background program in the MCU, including the initial parameter setting module, the source wave field extension module, the detector wave field extension module, and the cross-correlation imaging module. The source wave field extension module extends the source wave field to obtain the source extension wave field; the detector wave field extension module extends the detector wave field to obtain the detector extension wave field; the cross-correlation imaging module performs imaging condition calculations on the source extension wave field and the detector extension wave field.

[0087] A tunnel advance detection imaging method, the method comprising:

[0088] The source wave field is extended to obtain the source extended wave field;

[0089] Collect seismic data on both sides of the tunnel and on the face;

[0090] Taking the seismic data as the boundary condition to carry out the detector wave field extension, the detector extension wave field is obtained;

[0091] Perform imaging condition calculations on the source extended wavefield and the receiver extended wavefield.

[0092] In one embodiment, the specific process includes:

[0093] The first step is initial setup.

[0094] according to Figure 2 As shown in the figure, the observation system is set up, which is characterized by the even distribution of acquisition stations on both sides of the tunnel and the face. Six acquisition stations are placed at symmetrical positions on both sides of the tunnel, with a spacing of 12m. The first station is 8m away from the face, and the other two are placed at the three points of the face. The source is set at the two corners of the tunnel face, and two shots are performed respectively. The purpose of this is to obtain more comprehensive seismic data.

[0095] Subsequently, the initial position of the source and the detector, and the initial velocity model are input through the touch screen module 4 and acquired by the background calculation program. The network module 1 is used to send the initial setting parameters such as the collector sampling frequency and sampling duration to all collectors through UDP data broadcast, and the working status of the collector is checked.

[0096] The second step is to extend the source wave field.

[0097] Continuation is the computer simulation of the propagation of seismic waves in the medium, which is the basis of reverse time migration technology. Its principle is to solve the wave equation using numerical simulation methods. In order to better apply it to the calculation of the above-mentioned embedded control equipment, the following specific scheme is used for continuation.

[0098] 1) The two-dimensional first-order velocity stress elastic wave equation is selected as the basic equation for extension.

[0099] The elastic wave equation was chosen because it can better reflect the actual propagation of seismic waves in the medium. Therefore, the use of the elastic wave equation can simulate and emulate the underground medium more closely to the actual situation, and can also obtain seismic data containing more information.

[0100] Two dimensions are chosen because construction parties are often concerned about medium changes in the direction parallel to the tunnel face, but not changes in the vertical direction. Therefore, reducing the dimension and focusing on key information can more effectively reduce computing costs and improve computing efficiency.

[0101] The first-order velocity stress form is chosen because it can be combined with the staggered grid method in 3) below to solve the problem. This method can effectively reduce the amount of calculation and numerical dispersion. The first-order velocity stress elastic wave equation is shown as follows:

[0102]

[0103] In the formula, x is the direction parallel to the tunnel face, z is the direction perpendicular to the tunnel face, and t is the running time σ xx , σ yy , σ zz is the normal stress on the microelement, v x and v zis the velocity component of the microelement in the x and z directions, σ zx , σ xy , σ yz is the shear stress on the microelement, the Lame coefficients λ and μ and the propagation speed v of the elastic wave in the medium p and v s And the medium density ρ satisfies the following relationship:

[0104]

[0105]

[0106] 2) Select the main frequency f m =500Hz Ricker wavelet is used as the earthquake source when the earthquake source wave field is extended forward.

[0107] Ricker wavelet is a commonly used wavelet in seismic exploration, which can better reflect the basic characteristics of seismic wave propagation. 500Hz is the middle value of the commonly used frequency of explosion sources. Its formula is as follows:

[0108]

[0109] 3) Select the high-order staggered grid finite difference method as the numerical solution method.

[0110] The present invention uses a 6th order precision variable step size staggered grid finite difference method to solve the wave equation in 1). The principle is to use the difference in the computer to replace the differential of the formula in 1) and perform iterative numerical solution. The schematic diagram of the staggered grid is shown in Figure 3 Without loss of generality, taking the x direction as an example, the differential formula is as follows:

[0111]

[0112] Where u is the velocity or stress term to be determined, x0 is the differential point to be determined, Δx is the grid step of the staggered grid, the general accuracy requirement is 1m, 2N is the difference order, and a m is the difference coefficient.

[0113] In order to determine the appropriate calculation parameters, the present invention performs dispersion calculations on staggered grids of different orders and conventional grid finite difference methods. The dispersion curves are compared. Figure 4 As shown. Since the main frequency of the sub-wave is 500Hz, at the longitudinal wave velocity v pAt 4000m / s, Δx / λ=0.125. It can be seen that under the staggered grid operation of the 6th order progress, the dispersion improvement brought by further increasing the order is very small, but it will bring huge computational consumption, so the 6th order spatial difference is selected for numerical calculation. In addition, the present invention also proposes a variable step size idea to reduce computational consumption. We change the spatial step size of the density change part in the extension area to 1m, the spatial step size of the uniform medium area to 2m, and the spatial step size of the PML absorption area to 4m to reduce the number of computational grid points, thereby reducing the complexity of the calculation.

[0114] 4) Use finite difference coefficients in time and space domains.

[0115] For the differential coefficient a in 3) m , the present invention uses a time-space domain differential coefficient to solve, the differential coefficient can effectively reduce numerical dispersion and improve imaging accuracy. We use the following equation to solve the differential coefficient a m .

[0116]

[0117] Where 2N is the difference order, b1=1, b n Satisfy the following formula

[0118]

[0119] The parameters in the formula satisfy the following formula

[0120] c n =(cosθ) 2n-1 β n , d n =(sinθ) 2n-1 β n

[0121] When the differential order 2N = 6, the differential coefficient a is m =[1.167,-0.0632,0.0045].

[0122] 5) A PML absorbing boundary using an efficient absorption function.

[0123] Artificial boundaries are boundary problems that should not occur in computer numerical simulations due to the limited calculation area. PML (Perfectly Matched Layer) is a method to deal with false reflections caused by artificial boundaries in computer numerical simulations. The principle is to construct a suitable absorption function to reduce the false reflections in a limited area to a number and size that will not affect the normal numerical simulation. The present invention proposes an absorption function with higher absorption efficiency, and the expression is as follows:

[0124]

[0125] Where L is the total number of layers in the absorption region, k is the number of layers from the forward modeling region, and R = 0.0001. This absorption function can balance the reflections caused by the traditional absorption function, including the head reflection of the quadratic function and the tail reflection of the exponential function. Figure 5 Schematic diagram of the PML absorbing boundary.

[0126] After integrating all steps 1)-5), the differential iteration equation for any point in the simulation area is as follows:

[0127] The differential iteration equation on the x component is:

[0128]

[0129] The difference iteration equation on the z component is

[0130]

[0131] Where: d(x0,z0) is the absorption matrix, which is 0 in the simulation area and is determined by the absorption function in 5) in the absorption area. now represents the stress or velocity value at the point (x0,z0) at the current moment, and next represents the stress or velocity value to be calculated at the point (x0,z0) at the next moment. dif x () indicates to find the difference of the variable in the x direction, dif z () indicates differentiating the variable in the z direction, and the difference method is given in 3).

[0132] According to the second step specific plan 5), the neutron wave function in 2) is introduced as the boundary condition, and the differential formula in 3) and the differential coefficient in 4) are used to solve the two-dimensional first-order velocity stress elastic wave equation. After the solution is completed, the source wave field extension result S (x, z, t) is saved in the SD card. The advantage of this is that the extension process of the source wave field is completed before acquisition, which can greatly reduce the calculation time during imaging.

[0133] The third step is seismic data collection and acquisition.

[0134] The user turns on the transistor switch of the control module 2 through the touch screen module 4, triggering the circuit to conduct and stimulating the explosion source. Prior to this, the network module 1 uses UDP data broadcast to send out a synchronous acquisition signal in advance, so that all acquisition stations can collect at the sampling frequency set in the first step.

[0135] After the sampling time set in the first step, the network module 1 automatically uses the FTP protocol to download the seismic data collected by the collector back to the control device. After the transmission is completed, the data is saved in the SD card and the waveform is displayed on the touch screen module 4. The network communication tasks involved in the first and third steps are programmed based on the lightweight network communication protocol LWIP, which has fast communication speed, occupies less resources, and is a bottom-level encapsulation protocol suitable for embedded systems.

[0136] The fourth step is to extend the detector wave field.

[0137] The detector wavefield extension method is the same as the source wavefield extension method. It still uses the equation in the second step specific solution 5) to substitute the differential formula in 3) and the differential coefficient in 4) for solution, but the boundary condition is no longer the source function, but the multi-channel seismic signal collected in the third step. Finally, the obtained detector wavefield extension result R (x, z, t) is also saved in the SD card.

[0138] For the extension of the earthquake source wavefield, two boundary conditions are known. The first is the wavefield value at time t=0 and the wavelet function at the earthquake source (x0,z0), namely U(x,z,t=0) and U(x=x0,z=z0,t). Based on these, the extension calculation of the earthquake source wavefield can be performed to obtain the earthquake source extension wavefield S(x,z,t).

[0139] For the detector wave field extension, it is known that the detector (x d ,z d ) received at the earthquake signal data U(x=x d ,z=z d ,t), assuming t>t end When the wave field value is 0, that is, U(x,z,t>t end )=0, based on this, the detector wave field can be extended inversely to obtain the detector extended wave field R(x,z,t). The calculation of reverse extension is similar to that of forward extension. Forward extension is to use the finite difference method to calculate the detector wave field when U is known. t In the case of U t+1 , and the reverse extension is to use the same formula in the known U t+1 In the case of U t .

[0140] Step 5: Imaging.

[0141] Reverse time migration imaging is based on the principle of time consistency imaging, that is, if the forward extension and reverse extension wave peaks reach the same position at the same time, there must be a wave impedance difference interface at this place. Therefore, as long as certain operations are performed on the source wave field S and the detector wave field R, the imaging of the underground structure can be obtained. Such operations are called imaging conditions in reverse time migration imaging technology. The principle of reverse time migration imaging is as follows Figure 6 As shown:

[0142] In order to achieve good imaging effects and reduce the amount of calculation to be applicable to embedded systems, the present invention proposes the following imaging conditions:

[0143] Calculate travel time to filter direct waves. Since the observation system is fixed, the travel time of direct waves is easy to calculate. By calculating the travel time, direct waves can be filtered to reduce their impact on the imaging results, especially near-surface noise. Therefore, after obtaining seismic data in the third step, direct waves must be eliminated.

[0144] Use cross-correlation imaging conditions. Compared with traditional excitation time imaging conditions and amplitude ratio imaging conditions, this imaging condition does not need to traverse the continuation results to find the excitation time of each grid point, but directly cross-correlate the results of the forward continuation of the source wave field and the reverse continuation of the detector wave field to obtain the imaging result. The advantages of this are: first, the calculation speed is fast without traversal; second, there will be no ratio distortion due to the minimum value. However, this will take up a large storage space, and an external SD card is required to store the calculation results. The cross-correlation imaging condition is shown in the following formula:

[0145]

[0146] Use normalization operation. The purpose of this is to eliminate the ray noise generated by the calculation. This method requires the use of fast Fourier transform to first transform the source extension wave field S(x,z,t) and the detector extension wave field R(x,z,t) into the frequency domain to obtain S(x,z,ω) and R(x,z,ω), and then calculate the conjugate S of the frequency domain source extension wave field S(x,z,ω) * (x,z,ω). The operation expression of the normalized cross-correlation imaging condition is as follows:

[0147]

[0148] Use the prewitt operator to extract the boundary of the imaging result image(x,z).

[0149] image_grey(x,z)=prewitt(imagey(x,z))

[0150] The prewitt operator is a first-order differential operator with a relatively fast calculation speed. It can effectively extract the medium stratification and fracture boundary in the imaging result, and effectively display the result to the construction personnel in the form of a binary image. It is suitable for embedded systems with high real-time requirements. Finally, the obtained imaging result image_grey(x,z) is output on the touch screen module 4, and the result is saved in the SD card.

[0151] Figure 7-10It is used to illustrate the working process of the present invention. Figure 7 is a possible two-layer inclined tunnel model. Figure 8 is the seismic data collected by the observation system. Fig. 9 is the imaging result image(x,z) obtained by normalized cross-correlation operation, Fig.10 It is the imaging result image_grey(x,z) obtained after boundary extraction.

[0152] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A tunnel advance detection imaging method, characterized in that: The method includes: The source wave field is extended to obtain the source extended wave field; Collect seismic data on both sides of the tunnel and on the face; Taking the seismic data as the boundary condition to carry out the detector wave field extension, the detector extension wave field is obtained; Perform imaging condition calculation on the extended wave field of the source and the extended wave field of the receiver; The extended source wave field includes: The two-dimensional first-order velocity stress elastic wave equation is selected as the basic equation for extension. The two-dimensional first-order velocity stress elastic wave equation is: Where x is the direction parallel to the tunnel face, z is the direction perpendicular to the tunnel face, t is the running time, σ xx , σ yy , σ zz is the normal stress on the microelement, v x and v z are the velocity components of the microelement in the x and z directions respectively, σ zx , σ xy , σ yz is the shear stress on the microelement, the Lame coefficients λ and μ and the propagation speed of elastic waves in the medium, v p and v s And the medium density Satisfy the following relationship: Select the main frequency f m =500Hz Ricker wavelet as the source when the source wave field is extended; The high-order staggered grid finite difference method is used as the numerical solution method to solve the two-dimensional first-order velocity stress elastic wave equation to obtain the differential iterative equation of any point in the simulation area, where the differential formula is: Where u is the velocity or stress term to be determined, A0 is the differential point to be determined, ΔA is the grid step size of the staggered grid, 2N is the difference order, and a is m is the difference coefficient, A is x or z; Use the following equation to solve for the differential coefficient a m : Where 2N is the difference order, b1=1, b n Satisfy the following formula: The parameters in the formula satisfy: c n =(cosθ) 2n-1 β n , d n =(sinθ) 2n-1 β n , θ=π8; select the absorption function to bring into the Ricker wavelet as the boundary condition, the absorption function expression is as follows: Where L is the total number of layers in the absorption region, k is the number of layers from the forward modeling region, and R = 0.0001; The Ricker wavelet expression is: The difference iteration equations at any point include: The differential iteration equation on the x component is: The differential iteration equation on the z component is: Where: d(x0,z0) is the absorption matrix, which is 0 in the simulation area, now represents the stress or velocity value at the point (x0,z0) at the current moment, next represents the stress or velocity value to be calculated at the point (x0,z0) at the next moment, and dif x () indicates to find the difference of the variable in the x direction, dif z () means taking the difference of the variable in the z direction; The multi-channel collectors are evenly distributed on both sides of the tunnel and on the face. The collectors on both sides of the tunnel are symmetrically positioned, and one collector is placed at each of the three points of the face. The seismic sources are set at the two corners of the tunnel face.

2. The method according to claim 1, characterized in that Performing detector wavefield extension includes: The two-dimensional first-order velocity stress elastic wave equation is selected as the basic equation for extension. The two-dimensional first-order velocity stress elastic wave equation is: Where x is the direction parallel to the tunnel face, z is the direction perpendicular to the tunnel face, t is the running time, σ xx , σ yy , σ zz is the normal stress on the microelement, v x and v z are the velocity components of the microelement in the x and z directions respectively, σ zx , σ xy , σ yz is the shear stress on the microelement, the Lame coefficients λ and μ and the propagation speed of elastic waves in the medium, v p and v s And the medium density Satisfy the following relationship: Select the main frequency f m =500Hz Ricker wavelet as the source when the source wave field is extended; The high-order staggered grid finite difference method is used as the numerical solution method to solve the two-dimensional first-order velocity stress elastic wave equation to obtain the differential iterative equation of any point in the simulation area, where the differential formula is: Where u is the velocity or stress term to be determined, A0 is the differential point to be determined, ΔA is the grid step size of the staggered grid, 2N is the difference order, and a is m is the differential coefficient, and A is x or z.

3. The method according to claim 2, characterized in that The differential coefficient: Use the following equation to solve for the differential coefficient a m , Where 2N is the difference order, b1=1, b n Satisfy the following formula The parameters in the formula satisfy the following formula:

4. The method according to claim 2, characterized in that The imaging condition calculations for the source extended wave field and the receiver extended wave field include: Calculate travel time to filter direct waves and remove direct waves; Using the cross-correlation imaging condition, the results of the source wavefield extension and the receiver wavefield extension are cross-correlated. The cross-correlation imaging condition is as follows: Use normalization operation: Use fast Fourier transform to transform the source extension wave field S(x,z,t) and the detector extension wave field R(x,z,t) into the frequency domain to obtain S(x,z,ω) and R(x,z,ω), and then calculate the conjugate S of the frequency domain source extension wave field S(x,z,ω) * (x, z, ω), the operation expression of the normalized cross-correlation imaging condition is as follows: Use the prewitt operator to extract the boundary of the imaging result image(x,z).

Images