Three-dimensional reverse time extension method of seismic waves and high-resolution tunnel earthquake advance prediction method

Abstract

This invention discloses a three-dimensional reverse time-delay topology method for seismic waves and a high-resolution tunnel earthquake early prediction method, comprising the following steps: Step 1: Acquire the first arrival travel time information of seismic data excited and received by the tunnel wall, and calculate the P-wave velocity Vp and S-wave velocity Vs of the surrounding rock by combining the spacing of the shot points; Step 2: Using the calculated P-wave velocity Vp and S-wave velocity Vs as a velocity model, perform three-dimensional reverse time-delay topology calculation from the maximum time to the minimum time, during which virtual detectors are deployed at the tunnel face to receive the three-dimensional reverse time-delay topology wavefield; Step 3: Perform autocorrelation imaging on the three-dimensional reverse time-delay topology wavefield, further converting the seismic wavefield excited by the tunnel wall and received by the tunnel face into a seismic wavefield excited by the tunnel face and received by the tunnel face. This invention aims to suppress strong surface wave interference in the seismic data received by the tunnel wall, improve the resolution of reflected waves, and thus improve the accuracy of tunnel earthquake early prediction.

Description

Technical Field

[0001] This invention relates to the field of tunnel seismic detection and engineering survey technology, and in particular to a three-dimensional reverse time extension method for seismic waves and a high-resolution tunnel seismic advance prediction method. Background Technology

[0002] Tunnel earthquake advance prediction methods can effectively predict adverse geological bodies such as fault fracture zones and karst caves in front of the tunnel face, reducing the risk of tunnel excavation (Dickmann and Sander, 1996; Shi et al., 2014). Existing advanced geological prediction methods for tunnel seismic events mainly include 2-D tunnel seismic prediction (TSP) (Lin et al., 2006; Shi et al., 2014), integrated seismic imaging system (Borm et al., 2003), horizontal seismic profiling (HSP) (Inazaki et al., 1999), sonic soft ground probing (Kneib et al., 2000), tunnel seismic tomography (Zhao et al., 2006), tunnel seismic-while-drilling method (Petronio & Poletto 2002), and 3D true reflection tomography (3D TRT) (Shang et al., 2012). The methods described above all involve exciting seismic wave signals at the tunnel wall and receiving them at the tunnel wall. This is mainly because it is very difficult to deploy geophones at the tunnel face due to risks such as support and rockfall. Furthermore, due to the influence of the tunnel excavation damage zone (Feng et al., 2016; Ma Chunde et al., 2018), the seismic wave signals received by the tunnel wall suffer from severe surface wave interference. High-energy surface waves mask reflected waves, reducing the signal-to-noise ratio of the reflected wave data.

[0003] Current methods for processing seismic wave signals received through tunnel walls mainly consist of three modules: preprocessing, data imaging, and geological interpretation (Lu et al., 2020). The preprocessing module includes data editing, spectral analysis, bandpass filtering, automatic gain control, inverse Q-filtering, predictive deconvolution, and FK filtering. The data imaging module includes full waveform inversion and 2D / 3D migration imaging. The geological interpretation module includes data imaging interpretation and geological data analysis. Because strong-energy surface waves in the raw seismic data received through tunnel walls can mask reflected wave data, subsequent data processing methods are significantly less accurate in extracting reflected wave signals, leading to insufficient accuracy in geological prediction. Current methods for processing seismic wave signals received through tunnel walls still have room for improvement and refinement.

[0004] Previous tunnel earthquake early prediction methods, considering safety and construction efficiency, failed to deploy geophones at the tunnel face to receive seismic wave signals. Subsequent data processing methods also consistently processed seismic wave signals received from the tunnel walls. Combining the principles of seismic wave propagation at the tunnel excavation site and within the tunnel space, this invention proposes a high-resolution tunnel earthquake prediction method based on three-dimensional reverse time-delay topology. This method aims to suppress strong surface wave interference in the seismic data received from the tunnel walls, improve the resolution of reflected waves, and thus enhance the accuracy of tunnel earthquake early prediction.

[0005] References

[0006] Borm,G.,Giese,R.,Klose,C.,Mielitz,S.,Otto,P.&Bohlen,T.,2003.ISIS-integrated seismic imaging system for the geological prediction ahead ofunderground construction,65 th Conference&Exhibition,EAGE,Stavanger,Norway.

[0007] Dickmann,T.,Sander,BK,1996.Drivage concurrent tunnel seismicprediction.Felsbau,14,406-411.

[0008] Feng,X.T.,Zhang,C.Q.,Qiu,S.L.,Zhou,H.,Jiang,Q.,Li,S.J.,2016.Dynamicdesign method for deep hard rock tunnels and its application.Journal of RockMechanics and Geotechnical Engineering,8,443-461.

[0009] Inazaki,T.,Isahai,H.,Kawamura,S.,Kuruhashi,T.,Hayashi,H.,1999.Stepwise application of horizontal seismic profiling for tunnelprediction ahead of the face.The Leading Edge,18,1429-1431.

[0010] Kneib,G.,Kassel,A.&Lorenz,K.,2000.Automatic seismic prediction aheadof the tunnel bore machine.First Break,18,295-302.

[0011] Lin,C.N.,Jiao,Y.Y.,Liu,Q.S.,2006.Site experiment for predictinghazardous geological formations ahead of tunnel face.Key EngineeringMaterials,326,461-464.

[0012] Lu,X.L.,Liao,X.,Wang,Y.,Wang,G.M.,Fu,Z.H.,Tai,H.M.,2020.The tunnelseismic advance prediction method with wide illumination and a high signal-to-noise ratio.Geophysical Prospecting,68(8),2444-2458.

[0013] Petronio,L.&Poletto,F.,2002.Seismic-While-drilling by using tunnelboring-machine noise.Geophysics,67,1798-1809.

[0014] Shang,JL, Luo,

[0015] Shi,SS,Li,SC,Li,LP,Zhou,ZQ,Wang,J.,2014.Advance optimizedclassification and application of surrounding rock based on fuzzy analytichierarchy process and Tunnel Seismic Prediction.Automation in Construction,37,217-222.

[0016] Zhao,YG,Jiang,H.,Zhao,XP,2006.Tunnel seismic tomography method for geological prediction and its application.Applied Geophysics,3(2),69-74.

[0017] Ma Chunde, Guo Chunzhi, Fu Wei, Zhou Yanan. 2018. Study on the correspondence between rock wave velocity, stress and energy storage during loading and unloading within the elastic range. Mining and Metallurgical Engineering. Vol. 38, No. 1. Summary of the Invention

[0018] In view of at least one deficiency of the prior art, the purpose of this invention is to provide a three-dimensional reverse time extension method for seismic waves, which aims to suppress strong surface wave interference in seismic data received by tunnel walls, improve the resolution of reflected waves, and thus improve the accuracy of tunnel earthquake advance prediction.

[0019] To achieve the above objectives, the present invention adopts the following technical solution:

[0020] A three-dimensional reverse time-delay topology method for seismic waves, the key of which includes the following steps:

[0021] Step 1: Pick up the first arrival travel time information of the seismic data excited by the tunnel wall and received by the tunnel wall, and calculate the P-wave velocity Vp and S-wave velocity Vs of the surrounding rock by combining the spacing of the shot points;

[0022] Step 2: Using the calculated longitudinal wave velocity Vp and transverse wave velocity Vs as the velocity model, perform three-dimensional reverse time extension calculation from the maximum time to the minimum time, converting the wave field received by the tunnel wall into the wave field received by the virtual detector at the tunnel face. During the calculation process, virtual detectors are deployed at the tunnel face to receive the three-dimensional reverse time extension wave field.

[0023] Step 3: Perform autocorrelation imaging on the three-dimensional inverse time-delayed topology wavefield to further convert the seismic wavefield excited by the tunnel wall and received by the tunnel face into a seismic wavefield excited by the tunnel face and received by the tunnel face.

[0024] The specific derivation formula for the three-dimensional inverse time delay topology in step two is as follows:

[0025] The formula for calculating the decoupled first-order partial differential equation in the three-dimensional elastic wave equation is as follows:

[0026]

[0027] Differential equations (1) to (9) are discretized to obtain the difference formulas for the three-dimensional inverse time-delay topology:

[0028]

[0029]

[0030] For velocity components:

[0031]

[0032] In the above formula, i, j, k are the subscripts of the spatial directions X, Y, and Z. x, y, and z are the subscripts of the three components of the spatial directions X, Y, and Z, respectively. n is the time superscript. Δx, Δy, and Δz are the grid spacings of the X, Y, and Z directions, respectively, with a sampling interval of Δt. X, Y, and Z refer to three directions in the tunnel space. X is parallel to the tunnel face and points horizontally, Y is the direction of tunnel forward excavation, and Z is parallel to the tunnel face and points vertically downward. Formulas (11) to (18) have second-order precision in both space and time. Similarly, they can be derived to second-order precision in space, where L is an integer from 1 to 6. In equations (1) to (18), v x v y v z Represents the velocity components in the X, Y, and Z directions. (τ)xx ,τ yy ,τ zz ,τ xy ,τ yz ,τ xz ) is the stress component, τ xx Let τ be the stress component in the X direction. yy Let τ be the stress component in the Y direction. zz Let τ be the stress component in the Z direction. xy τ is the stress component of the resultant force in the X and Y directions. yz τ is the stress component of the resultant force in the Y and Z directions. xz Let f be the stress component of the resultant force in the X and Z directions. x ,f y ,f z ) represents the body force source components in the X, Y, and Z directions. In the three-dimensional reverse time-delay topology, this represents the three-component seismic data received by the tunnel wall. ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the P-wave velocity Vp and S-wave velocity Vs is as follows:

[0033] λ=ρ(Vp 2 -2Vs 2 (19);

[0034] μ=ρVs 2 (20);

[0035] Among them, the calculation formulas (9) to (18) are used to calculate the reverse time extension from the maximum time to the minimum time. During the calculation process, virtual detectors are set up on the working face to receive the wave field of the reverse time extension.

[0036] Data processing techniques such as bandpass filtering, automatic gain control, predictive deconvolution, full waveform inversion, and reverse time offset are implemented.

[0037] The key to the three-dimensional reverse time extension method for seismic waves is that it also includes step four, which involves performing bandpass filtering, automatic gain control, predictive deconvolution, full waveform inversion, and reverse time migration data processing based on the results of step three.

[0038] A high-resolution tunnel earthquake early prediction method, including the aforementioned three-dimensional reverse time-delay topology method for seismic waves, is characterized by the following key features:

[0039] A virtual detector is set up at the tunnel face. The observation system that is excited by the tunnel wall and received by the tunnel wall is converted into an observation system that is excited by the tunnel face and received by the tunnel face through three-dimensional inverse time extension calculation. The wave field received by the tunnel wall is then recursively pushed to the virtual detector position at the tunnel face, and finally converted into an observation system that is excited by the tunnel face and received by the tunnel face.

[0040] The specific data processing steps are as follows:

[0041] Step A: Data preprocessing: including data editing, bandpass filtering, FK filtering, and automatic gain control;

[0042] Step B: Data Imaging: 3D inverse time extension topology, predictive deconvolution, full waveform inversion, and 2D and 3D inverse time migration imaging, wherein the 3D inverse time extension topology uses the aforementioned seismic wave 3D inverse time extension topology method to process the data;

[0043] Step C: Geological Interpretation: Combining the full waveform inversion and reverse time migration imaging results, extracting adverse geological bodies, classifying surrounding rock grades, and providing support recommendations.

[0044] Significant effects: This invention provides a three-dimensional reverse time extension method for seismic waves and a high-resolution tunnel earthquake advance prediction method, which can suppress strong surface wave interference in seismic data received by the tunnel wall, improve the resolution of reflected waves, and thus improve the accuracy of tunnel earthquake advance prediction. Attached Figure Description

[0045] Figure 1 Schematic diagram of the three-dimensional reverse time extension principle of tunnels;

[0046] Figure 2 Flowchart of data processing for a high-resolution tunnel earthquake early prediction method based on 3D inverse time-delay topology;

[0047] Figure 3 : Three-dimensional model of tunnel karst caves and weak interlayers;

[0048] Figure 4 Three-component seismic records received by tunnel wall excitation-tunnel wall detector G1: (a) X component; (b) Y component; (c) Z component;

[0049] Figure 5 Three-component seismic records received by tunnel wall excitation-tunnel wall detector G2: (a) X component; (b) Y component; (c) Z component;

[0050] Figure 6 Three-component seismic records excited by the tunnel wall and received at the tunnel face: (a) X component; (b) Y component; (c) Z component;

[0051] Figure 7 Three-dimensional reverse time-delayed seismic records from the tunnel wall excitation-face reception: (a) X component; (b) Y component; (c) Z component;

[0052] Figure 8 Three-component seismic autocorrelation records excited by the tunnel wall and received at the working face: (a) X component; (b) Y component; (c) Z component;

[0053] Figure 9Three-component seismic autocorrelation records of the tunnel, generated by three-dimensional reverse time delay and transmitted to the tunnel wall and received at the tunnel face: (a) X component; (b) Y component; (c) Z component;

[0054] Figure 10 Flowchart of the three-dimensional reverse time extension method for seismic waves. Detailed Implementation

[0055] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0056] This invention relates to a three-dimensional reverse time-delay topology method for seismic waves and a high-resolution tunnel earthquake advance prediction method, which aims to suppress strong surface wave interference in seismic data received by the tunnel wall, improve the resolution of reflected waves, and thus improve the accuracy of tunnel earthquake advance prediction.

[0057] The present invention adopts the following technical solution:

[0058] 1. A three-dimensional reverse time-delayed seismic wave method and a high-resolution tunnel earthquake advance prediction method, the key of which is, taking a conventional tunnel earthquake advance prediction observation system as an example, such as... Figure 1 and Figure 3 As shown, the specific data acquisition method involves firing 24 shots and using 2 geophones to receive seismic wave signals from the tunnel wall. The specific process of the tunnel seismic 3D reverse time-delay topology is as follows: a virtual geophone is set up at the tunnel face to receive seismic wave signals. The seismic wave signals received from the tunnel wall are then reconstructed to the tunnel face position through 3D reverse time-delay topology calculations. By changing the receiving position of the seismic wave signals, the observation system of tunnel wall firing and receiving is transformed into an observation system of tunnel wall firing and receiving at the tunnel face. The specific derivation formula for the 3D reverse time-delay topology is as follows:

[0059] The formula for calculating the decoupled first-order partial differential equation in the three-dimensional elastic wave equation is as follows:

[0060]

[0061] Equations (1) to (9) are the three-dimensional numerical simulation equations for elastic waves. In the forward modeling, the wave field at time T+1 is recursively calculated using the wave field at time T. Unlike the forward modeling calculation process, in the three-dimensional reverse time extension, the wave field at time T is recursively calculated using the wave field at time T-1. The difference discretization equations are significantly different from the forward modeling equations. By differentially discretizing the differential equations (1) to (9), the difference formulas for the three-dimensional reverse time extension are obtained:

[0062]

[0063]

[0064] For velocity components:

[0065]

[0066] In the above formula, i, j, k are the subscripts of the X, Y, and Z spatial directions. x, y, and z are the subscripts of the three components of the X, Y, and Z spatial directions, respectively. n is the time superscript. Δx, Δy, and Δz are the grid spacings in the X, Y, and Z directions, respectively, with a sampling interval of Δt. X, Y, and Z refer to three directions in the tunnel space: X is parallel to the tunnel face and points horizontally; Y is the direction of forward tunnel excavation; and Z is parallel to the tunnel face and points vertically downward.

[0067] Equations (11) to (18) have second-order precision in both space and time. Similarly, they can be derived to second-order precision in space, where L is an integer from 1 to 6. In equations (1) to (18), v x v y v z Represents the velocity components in the X, Y, and Z directions. (τ) xx ,τ yy ,τ zz ,τ xy ,τ yz ,τ xz ) is the stress component, where τ xx Let τ be the stress component in the X direction. yy Let τ be the stress component in the Y direction. zz Let τ be the stress component in the Z direction. xy τ is the stress component of the resultant force in the X and Y directions. yz τx represents the stress components of the resultant force in the Y and Z directions, and τxz represents the stress components of the resultant force in the X and Z directions. (f) x ,f y ,f z ) represents the body force source components in the X, Y, and Z directions. The three-dimensional reverse time-delay topology represents 24 three-component seismic data received from the tunnel wall. ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the P-wave velocity Vp and S-wave velocity Vs is as follows:

[0068] λ=ρ(Vp 2 -2Vs 2 (19);

[0069] μ=ρVs 2 (20).

[0070] The present invention will now be described in further detail with reference to the accompanying drawings.

[0071] like Figures 1-10 As shown, this invention relates to a three-dimensional reverse time-delay topology method for seismic waves and a high-resolution method for predicting tunnel earthquakes in advance. Figure 1As shown, a virtual detector is set up at the tunnel face. The observation system, which is based on tunnel wall excitation and tunnel wall reception, is transformed into a tunnel face excitation and tunnel face reception observation system through three-dimensional reverse time extension calculation. The wavefield received by the tunnel wall is then recursively extrapolated to the virtual detector location at the tunnel face, ultimately resulting in a tunnel face excitation and tunnel face reception observation system. The specific advantages of this invention are:

[0072] (1) An observation system that can be implemented by tunnel wall excitation and face reception can be realized without increasing the data acquisition cost at the tunnel construction site;

[0073] (2) By using three-dimensional inverse time-delay topology calculation to suppress strong surface wave interference in the seismic data received by the tunnel wall, the resolution of the forward direction reflected wave is improved, thereby improving the prediction accuracy.

[0074] like Figure 2 As shown, the specific data processing steps are as follows:

[0075] (1) Preprocessing section: including data editing, bandpass filtering, FK filtering, automatic gain control, etc.;

[0076] (2) Data imaging section: 3D inverse time delay topology, predictive deconvolution, full waveform inversion and 2D and 3D offset imaging, etc.

[0077] like Figure 1 , Figure 3 , Figure 10 As shown, the detailed processing part of the three-dimensional inverse time-delay topology is as follows:

[0078] 1) Pick up the first arrival travel time information of the seismic data excited by the tunnel wall and received by the tunnel wall, and calculate the P-wave velocity Vp and S-wave velocity Vs of the surrounding rock by combining the spacing of the shot points;

[0079] 2) Using the calculated longitudinal wave velocity Vp and transverse wave velocity Vs as the velocity model, the calculation formulas (9) to (18) are used to calculate the reverse time extension from the maximum time to the minimum time. During the calculation process, virtual detectors are set up on the working face to receive the reverse time extension wave field.

[0080] 3) Autocorrelation imaging is performed on the inverse time-delayed topology wavefield to further convert the seismic wavefield excited by the tunnel wall and received by the tunnel face into the seismic wavefield excited by the tunnel face and received by the tunnel face.

[0081] 4) Based on the results of step 3), perform data processing such as bandpass filtering, automatic gain, predictive deconvolution, full waveform inversion, and reverse time offset.

[0082] (3) Geological interpretation: Combining the full waveform inversion and reverse time migration imaging results, we will extract adverse geological bodies, classify surrounding rock grades, provide support suggestions, and continuously monitor the excavation results to verify the prediction conclusions.

[0083] This invention analyzes the propagation process and wavefield characteristics of elastic waves in a three-dimensional tunnel space, and creatively proposes a three-dimensional reverse time extension method. This method transforms the original observation system, which is based on tunnel wall excitation and tunnel wall reception, into an observation system based on tunnel wall excitation and face reception. This suppresses strong surface wave interference in the tunnel wall received wavefield and improves the resolution of reflected waves in the leading direction.

[0084] The three-dimensional reverse time extension method is a key step in the processing of high-resolution tunnel seismic data and is within the scope of patent protection. The three-dimensional reverse time extension method is not limited to the observation system with 24 shot points and 2 detectors as described in the patent. The three-dimensional reverse time extension method is also within the scope of patent protection for data collected by other tunnel seismic observation systems.

[0085] Model Experiment

[0086] 1. Tunnel seismic data acquisition

[0087] To verify the method of this invention, a three-dimensional model of a karst cave and a weak interlayer was designed, such as... Figure 3 As shown. G1 and G2 are two geophone locations, and S1, S2, to S24 are 24 shot points. The P-wave velocity of the surrounding rock is 3000 m / s, the S-wave velocity is 1730 m / s, and the density is 2.3 g / cm³. 3 The longitudinal wave velocity in the tunnel space is 340 m / s, the transverse wave velocity is 0 m / s, and the density is 1.0 g / cm³. 3 The cave has a diameter of 10m, a longitudinal wave velocity of 800m / s, a transverse wave velocity of 375m / s, and a density of 1.4g / cm³. 3 The weak interlayer has a thickness of 8 m, a P-wave velocity of 1700 m / s, a S-wave velocity of 1000 m / s, and a density of 1.8 g / cm³. 3 The data sampling interval is 40 μs, and the number of sampling points per channel is 4000. Simulated 24-shot data is used, with three-component seismic wave signals received at positions G1 and G2. The three-component seismic wave signal received by the G1 detector is as follows: Figure 4 As shown. The three-component seismic wave signal received by the G2 detector is as follows. Figure 5 As shown. The symbols D represent direct wave, S represent surface wave, A represent acoustic wave, and R represent reflected wave in the leading direction. When elastic waves propagate through the tunnel wall, they generate strong surface wave interference, affecting the resolution of reflected waves in the leading direction. This makes it difficult to identify and accurately extract reflected wave signals, severely impacting the prediction accuracy of adverse geological formations.

[0088] To facilitate comparison of the three-dimensional reverse-time continuation results, geophones were simultaneously deployed at the tunnel face during the numerical simulation, with the geophone height parallel to the positions of the shot point and the receiver point. The three-component wavefield records received by the geophones at the tunnel face are as follows: Figure 6 As shown. The wavefield record received by the G2 detector in the tunnel wall was reconstructed to the tunnel face location using the three-dimensional reverse time delay topology method mentioned in this invention. The result is as follows. Figure 7 As shown. For easier comparison, [the text is incomplete]. Figure 6 and Figure 7 Autocorrelation imaging was performed on the wavefield records in the tunnel to convert the wavefield records excited at the tunnel wall and received at the tunnel face into wavefield records excited at the tunnel face and received at the tunnel face. The results are as follows: Figure 8 and Figure 9 As shown. Figure 8 and Figure 9 There is a clear similarity, especially in the Y component. Compared to Figure 5 , Figure 9 The resolution of reflected waves is significantly improved, and surface wave interference is reduced. Numerical simulation results demonstrate that the resolution of the wavefield record after three-dimensional reverse time-delay topology is close to that of the wavefield record received at the tunnel face, and superior to that of the wavefield record received at the tunnel wall.

[0089] The above invention is currently primarily applied in the field of tunnel earthquake early warning.

[0090] Finally, it should be noted that the above are only specific embodiments of the present invention. Of course, those skilled in the art can make modifications and variations to the present invention. If these modifications and variations fall within the scope of the claims of the present invention and their equivalents, they should be considered as being within the protection scope of the present invention.

Claims

1. A three-dimensional reverse time-delay topology method for seismic waves, characterized in that, Includes the following steps: Step 1: Pick up the first arrival travel time information of the seismic data excited by the tunnel wall and received by the tunnel wall, and calculate the P-wave velocity Vp and S-wave velocity Vs of the surrounding rock by combining the spacing of the shot points; Step 2: Using the calculated longitudinal wave velocity Vp and transverse wave velocity Vs as velocity models, perform three-dimensional reverse time extension calculations from the maximum time to the minimum time. During the calculation process, virtual detectors are deployed on the tunnel face to receive the three-dimensional reverse time extension wave field. Step 3: Perform autocorrelation imaging on the three-dimensional inverse time-delayed topology wavefield to further convert the seismic wavefield excited by the tunnel wall and received by the tunnel face into a seismic wavefield excited by the tunnel face and received by the tunnel face.

2. The three-dimensional reverse time-delay topology method for seismic waves according to claim 1, characterized in that: In step two, the derivation formula for the three-dimensional inverse time delay topology is as follows: The formula for calculating the decoupled first-order partial differential equation in the three-dimensional elastic wave equation is as follows: Differential equations (1) to (9) are discretized to obtain the difference formulas for the three-dimensional inverse time-delay topology: For velocity components: In the above formulas, i, j, k are the subscripts of the spatial directions X, Y, and Z; n is the time superscript; Δx, Δy, and Δz are the grid spacings of the X, Y, and Z directions, respectively, and the sampling interval is Δt; X, Y, and Z refer to three directions in the tunnel space, X is parallel to the tunnel face and points horizontally, Y is the tunnel forward excavation direction, and Z is parallel to the tunnel face and points vertically downward; Formulas (11) to (18) have second-order precision in both space and time, and can be similarly derived to second-order precision in space, where L is an integer from 1 to 6; in equations (1) to (18), v x v y v z Represents the velocity components in the X, Y, and Z directions; (τ) xx ,τ yy ,τ zz ,τ xy ,τ yz ,τ xz ) is the stress component, where τ xx Let τ be the stress component in the X direction. yy Let τ be the stress component in the Y direction. zz Let τ be the stress component in the Z direction. xy τ is the stress component of the resultant force in the X and Y directions. yz τ is the stress component of the resultant force in the Y and Z directions. xz The stress components of the resultant force in the X and Z directions; (f x ,f y ,f z ) represents the body force source components in the X, Y, and Z directions. In the three-dimensional reverse time-delay topology, this represents the three-component seismic data received by the tunnel wall. ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the P-wave velocity Vp and S-wave velocity Vs is as follows: λ=ρ(Vp 2 -2Vs 2 ) (19); μ=ρVs 2 (20); Among them, the calculation formulas (9) to (18) are used to calculate the reverse time extension from the maximum time to the minimum time. During the calculation process, virtual detectors are set up on the working face to receive the wave field of the reverse time extension.

3. The three-dimensional reverse time-delay topology method for seismic waves according to claim 1, characterized in that: It also includes step four, which involves performing bandpass filtering, automatic gain control, predictive deconvolution, full waveform inversion, and inverse time offset data processing based on the results of step three.

4. A high-resolution tunnel earthquake advance prediction method, comprising the three-dimensional reverse time delay extrapolation method for seismic waves as described in any one of claims 1-3, characterized in that: A virtual detector is set up at the tunnel face. The observation system that is excited by the tunnel wall and received by the tunnel wall is converted into an observation system that is excited by the tunnel face and received by the tunnel face through three-dimensional inverse time extension calculation. The wave field received by the tunnel wall is then recursively pushed to the virtual detector position at the tunnel face, and finally converted into an observation system that is excited by the tunnel face and received by the tunnel face. The specific data processing steps are as follows: Step A: Data preprocessing: including data editing, bandpass filtering, FK filtering, and automatic gain control; Step B: Data Imaging: 3D inverse time extension topology, predictive deconvolution, full waveform inversion and 2D and 3D inverse time migration imaging, wherein the 3D inverse time extension topology uses the seismic wave 3D inverse time extension topology method described in any one of claims 1-3 to process the data; Step C: Geological Interpretation: Combining the full waveform inversion and reverse time migration imaging results, extracting adverse geological bodies, classifying surrounding rock grades, and providing support recommendations.

Images