Seismic wave three-dimensional inverse time continuation method and high-resolution tunnel earthquake advanced forecasting method

Abstract

The invention discloses a seismic wave three-dimensional reverse time continuation method and a high-resolution tunnel earthquake advanced forecasting method, and the method comprises the following steps: 1, picking up the first arrival travel time information of tunnel wall excitation-tunnel wall receiving earthquake data, and calculating the longitudinal wave speed Vp and the transverse wave speed Vs of surrounding rock through combining the distance between shot points; 2, taking the calculated longitudinal wave velocity Vp and transverse wave velocity Vs as velocity models, starting three-dimensional inverse time continuation calculation from the maximum moment to the minimum moment, and arranging a virtual detector on a tunnel face to receive a three-dimensional inverse time continuation wave field in the calculation process; and step 3, implementing autocorrelation imaging on the wave field after three-dimensional reverse time continuation, and further converting a seismic wave field excited by the tunnel wall and received by the tunnel face into a seismic wave field excited by the tunnel face and received by the tunnel face. The invention aims to suppress strong surface wave interference in seismic data received by a tunnel wall, improve the resolution of reflected waves and further improve the advanced prediction precision of tunnel earthquakes.

Description

Technical Field

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

[0002] The tunnel earthquake advance prediction method can effectively predict adverse geological bodies such as fault fracture zones and karst caves in front of the tunnel face, thereby reducing the risk of tunnel excavation construction (Dickmann and Sander, 1996; Shi et al., 2014). The existing tunnel seismic geological prediction mainly includes 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 (Tunnel seismic tomography) (Zhao et al., 2006), tunnel seismic-while drilling method (Petronio & Poletto 2002) and 3D True Reflection Tomography (True Reflection Tomography - 3D TRT) (Shang et al., 2012). The above methods basically excite the tunnel wall to receive seismic wave signals, mainly because it is very difficult to deploy geophones at the tunnel face due to risks such as support and block falling. 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 are seriously interfered by surface waves. The strong energy surface waves cover the reflected waves, reducing the signal-to-noise ratio of the reflected wave data.

[0003] At present, the processing methods of seismic wave signals received by tunnel walls are mainly divided into three modules: preprocessing, data imaging, and geological interpretation (Lu et al., 2020). The preprocessing part includes: data editing, spectrum analysis, bandpass filtering, automatic gain, inverse Q filtering, predictive deconvolution, FK filtering, etc.; the data imaging part includes: full waveform inversion and two-dimensional and three-dimensional offset imaging; the geological interpretation module includes: data imaging interpretation and geological data analysis. Because the strong energy surface waves in the original seismic data received by the tunnel wall will cover the reflected wave data, the later data processing methods are obviously insufficient in accurately extracting the reflected wave signals, resulting in insufficient geological prediction accuracy. The current processing methods of seismic wave signals received by tunnel walls still have room for improvement.

[0004] In the early stage, the tunnel earthquake advance prediction method failed to deploy detectors at the tunnel face to receive seismic wave signals due to safety and construction efficiency. The later data processing method has been processing the seismic wave signals received by the tunnel wall. Combining the seismic wave propagation principle at the tunnel excavation site and in the tunnel space, the present invention proposes a high-resolution tunnel earthquake prediction method based on three-dimensional reverse time extension, which aims to suppress the strong surface wave interference in the seismic data received by the tunnel wall, improve the resolution of the reflected wave, and thus improve the accuracy of tunnel earthquake advance 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 defect of the prior art, the object of the present invention is to provide a three-dimensional reverse time extension method for seismic waves, aiming to suppress strong surface wave interference in the seismic data received by the tunnel wall, improve the resolution of the reflected wave, and thereby improve the accuracy of advance prediction of tunnel earthquakes.

[0019] In order to achieve the above object, the present invention adopts the following technical solution:

[0020] A three-dimensional inverse time extension method for seismic waves, the key of which is that it includes the following steps:

[0021] Step 1: Pick up the first arrival travel time information of the seismic data of tunnel wall excitation and tunnel wall reception, and calculate the longitudinal wave velocity Vp and shear wave velocity Vs of the surrounding rock in combination with the spacing of the shot points;

[0022] Step 2: Using the calculated longitudinal wave velocity Vp and shear wave velocity Vs as the velocity model, three-dimensional reverse time extension calculation is performed from the maximum time to the minimum time, and the wave field received by the tunnel wall is converted into the wave field received by the virtual geophone on the tunnel face. During the calculation process, virtual geophones are arranged on the tunnel face to receive the wave field of three-dimensional reverse time extension;

[0023] Step 3: Perform autocorrelation imaging on the wave field after three-dimensional reverse time extension, and further convert the seismic wave field of tunnel wall excitation and tunnel face reception into the seismic wave field of tunnel face excitation and tunnel face reception.

[0024] The specific derivation formula of three-dimensional inverse time extension in step 2 is:

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

[0026]

[0027] By discretizing the differential equations (1) to (9), we obtain the differential formula for three-dimensional inverse time extension:

[0028]

[0029]

[0030] For velocity components:

[0031]

[0032] 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 superscript of time. Δx, Δy and Δz are the grid spacings in the X, Y and Z directions respectively, and the sampling interval is Δt. X, Y and Z refer to the three directions in the tunnel space, X is parallel to the face and points to the horizontal direction, Y is the forward excavation direction of the tunnel, and Z is parallel to the face and points to the vertical downward direction. Formulas (11) to (18) are second-order accuracy in space and time, and can be derived to 2L-order in space by the same logic, 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 is the stress component in the X direction, τ yy is the stress component in the Y direction, τ zz is 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 is the stress component of the resultant force in the X and Z directions. (f x ,f y ,f z ) is the body force source component in the X, Y and Z directions, is the three-component seismic data received by the tunnel wall during three-dimensional reverse time extension, ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the longitudinal wave velocity Vp and the shear wave velocity Vs is:

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

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

[0035] The calculation formulas (9) to (18) are used to perform reverse time delay calculation from the maximum time to the minimum time. During the calculation process, virtual detectors are arranged on the tunnel face to receive the reverse time delayed wave field.

[0036] Implement data processing such as bandpass filtering, automatic gain, predictive deconvolution, full waveform inversion, and reverse time migration.

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

[0038] A high-resolution tunnel earthquake advance prediction method includes the three-dimensional inverse time extension method of seismic waves, the key of which is:

[0039] A virtual geophone is set at the tunnel face position, and the tunnel wall excitation-tunnel wall reception observation system is calculated through three-dimensional reverse time extension, and the wave field received by the tunnel wall is reversely time-recursed to the position of the virtual geophone at the tunnel face, and finally converted into a tunnel face excitation-tunnel face reception observation system;

[0040] The specific data processing steps are:

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

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

[0043] Step C: Geological interpretation: Combine the full waveform inversion and reverse time migration imaging results to extract unfavorable geological bodies, classify surrounding rock grades, and provide support recommendations.

[0044] Significant effect: The present invention provides a three-dimensional reverse time-delay 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. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0046] Figure 2 :Data processing flow chart of high-resolution tunnel earthquake advance prediction method based on three-dimensional reverse time extension;

[0047] Figure 3 : Tunnel three-dimensional cave and weak interlayer model;

[0048] Figure 4 : Tunnel wall excitation - three-component seismic record received by tunnel wall geophone G1: (a) X component; (b) Y component; (c) Z component;

[0049] Figure 5 : Tunnel wall excitation - three-component seismic record received by tunnel wall geophone G2: (a) X component; (b) Y component; (c) Z component;

[0050] Figure 6 : Three-component seismic records of tunnel wall excitation and face reception: (a) X component; (b) Y component; (c) Z component;

[0051] Figure 7 : Three-component seismic records of tunnel 3D reverse time extension to tunnel wall excitation and face reception: (a) X component; (b) Y component; (c) Z component;

[0052] Figure 8 : Three-component seismic autocorrelation records of tunnel wall excitation and face reception: (a) X component; (b) Y component; (c) Z component;

[0053] Fig. 9: Three-component seismic autocorrelation records of tunnel 3D reverse time extension to tunnel wall excitation and face reception: (a) X component; (b) Y component; (c) Z component;

[0054] Fig.10 :Flowchart of the three-dimensional reverse time extension method of seismic waves. DETAILED DESCRIPTION

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

[0056] The present invention relates to a three-dimensional reverse time extension method for seismic waves and a high-resolution tunnel earthquake advance prediction method, which are intended 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 extension method for seismic waves and a high-resolution tunnel earthquake advance prediction method, the key lies in taking a conventional tunnel earthquake advance prediction observation system as an example, such as Figure 1 and Figure 3 As shown in the figure, the specific data acquisition method is to excite 24 shots on the tunnel wall and 2 geophones to receive seismic wave signals. The specific process of three-dimensional reverse time extension of tunnel seismic is: set up a virtual geophone at the tunnel face to receive seismic wave signals, restore the seismic wave signals received by the tunnel wall to the tunnel face through three-dimensional reverse time extension calculation, change the receiving position of the seismic wave signal, and convert the observation system of tunnel wall excitation-tunnel wall reception into an observation system of tunnel wall excitation-tunnel face reception. The specific derivation formula of three-dimensional reverse time extension is:

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

[0060]

[0061] Formulas (1) to (9) are three-dimensional numerical simulation equations for elastic waves. In forward simulation, the wave field at time T is used to recursively calculate the wave field at time T+1. Different from the forward simulation calculation process, the wave field at time T-1 is used to recursively calculate the wave field at time T in three-dimensional reverse time extension. The differential discretization equation is significantly different from the forward simulation equation. The differential equations (1) to (9) are discretized to obtain the differential formula for three-dimensional reverse time extension:

[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, and the sampling interval is Δt. X, Y and Z refer to the three directions in the tunnel space, X is parallel to the face and points to the horizontal direction, Y refers to the forward excavation direction of the tunnel, and Z is parallel to the face and points to the vertical downward direction.

[0067] Equations (11) to (18) are second-order in space and time, and can be derived to 2L-order in space by the same logic, 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 is the stress component in the X direction, τ yy is the stress component in the Y direction, τ zz is 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, and τxz is the stress component of the resultant force in the X and Z directions. (f x ,f y ,f z ) is the body force source component in the X, Y and Z directions, 24 three-component seismic data received by the tunnel wall during three-dimensional reverse time extension, ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the longitudinal wave velocity Vp and the shear wave velocity Vs is as follows:

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

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

[0070] The present invention will be further described in detail below in conjunction with the accompanying drawings.

[0071] like Figure 1-Figure 10 As shown, the present invention relates to a three-dimensional reverse time extension method for seismic waves and a high-resolution tunnel earthquake advance prediction method. Figure 1As shown, a virtual detector is set at the tunnel face position, and the tunnel wall excitation-tunnel wall reception observation system is calculated by three-dimensional reverse time extension, and the wave field received by the tunnel wall is reversely recursively transferred to the tunnel face virtual detector position, and finally converted into a tunnel face excitation-tunnel face reception observation system. The specific advantages of the present invention are:

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

[0073] (2) The strong surface wave interference in the seismic data received by the tunnel wall is suppressed through three-dimensional reverse time extension calculation, and the resolution of the reflected waves in the leading direction is improved, thereby improving the prediction accuracy.

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

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

[0076] (2) Data imaging: 3D inverse time extension, predictive deconvolution, full waveform inversion, and 2D and 3D migration imaging;

[0077] like Figure 1 , Figure 3 , Fig.10 As shown, the detailed processing part of the three-dimensional inverse time extension part is:

[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 longitudinal wave velocity Vp and shear wave velocity Vs of the surrounding rock in combination with the spacing of the shot points;

[0079] 2) Taking the calculated P-wave velocity Vp and S-wave velocity Vs as the velocity model, the calculation formulas (9) to (18) are used to perform reverse time extension calculation from the maximum time to the minimum time. During the calculation process, virtual geophones are arranged on the tunnel face to receive the reverse time extended wave field.

[0080] 3) Autocorrelation imaging is performed on the wave field after reverse time extension, and the seismic wave field of tunnel wall excitation and tunnel face reception is further converted into the seismic wave field of tunnel face excitation and tunnel face reception;

[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 migration.

[0082] (3) Geological interpretation: Combined with the full waveform inversion and reverse time migration imaging results, unfavorable geological bodies are extracted, surrounding rock grades are classified, support suggestions are given, and the excavation results are continuously monitored to verify the prediction conclusions.

[0083] The present invention analyzes the propagation process and wave field characteristics of elastic waves in three-dimensional tunnel space, and creatively proposes a three-dimensional reverse time-delay method, which converts the original tunnel wall excitation-tunnel wall reception observation system into a tunnel wall excitation-face reception observation system, suppresses the strong surface wave interference in the tunnel wall receiving wave field, and improves the resolution of the reflected wave in the leading direction.

[0084] The three-dimensional reverse time extension method, as a key step in high-resolution tunnel seismic data processing, falls within the scope of patent protection. The three-dimensional reverse time extension method is not limited to the observation system of 24 shot point excitation and 2 detector reception in the patent book. The three-dimensional reverse time extension method in data collected by other tunnel seismic observation systems also falls within the scope of patent protection.

[0085] Model experiment

[0086] 1. Tunnel seismic data acquisition

[0087] In order to verify the method of the present invention, a three-dimensional cave and weak interlayer model was designed, such as Figure 3 As shown. G1 and G2 are two geophone positions, S1, S2, to S24 are 24 shot points. The longitudinal wave velocity of the surrounding rock is 3000m / s, the shear wave velocity is 1730m / s, and the density is 2.3g / cm 3 The longitudinal wave velocity in the tunnel space is 340m / s, the transverse wave velocity is 0m / s, and the density is 1.0g / cm 3 The diameter of the cave is 10m, the longitudinal wave velocity is 800m / s, the transverse wave velocity is 375m / s, and the density is 1.4g / cm 3 The thickness of the weak interlayer is 8m, the longitudinal wave velocity is 1700m / s, the shear wave velocity is 1000m / s, and the density is 1.8g / cm 3 The data sampling interval is 40μs, and the number of sampling points per channel is 4000. Simulate 24-shot data and receive three-component seismic wave signals at positions G1 and G2. The three-component seismic wave signals received by the G1 detector are as follows: Figure 4 The three-component seismic wave signal received by the G2 detector is shown in Figure 5 As shown in the figure, the symbol D is the direct wave, S is the surface wave, A is the acoustic wave, and R is the reflected wave in the leading direction. When the elastic wave propagates on the tunnel wall, it will produce strong surface wave interference, which will affect the resolution of the reflected wave in the leading direction, making it difficult to identify and accurately extract the reflected wave signal, which will seriously affect the prediction accuracy of the poor geological body.

[0088] In order to facilitate the comparison of 3D reverse time extension results, geophones were deployed at the tunnel face during numerical simulation, and the height of the geophones was parallel to the positions of the shot points and the geophones. The three-component wave field records received by the geophones at the tunnel face are shown in Figure 2. Figure 6 The wave field record received by the tunnel wall G2 detector is restored to the tunnel face position using the three-dimensional reverse time extension method mentioned in the present invention. The result is as follows: Figure 7 For the convenience of comparison, Figure 6 and Figure 7 The wave field records in the tunnel are used for autocorrelation imaging, and the wave field records of tunnel wall excitation and tunnel face reception are converted into wave field records of tunnel face excitation and tunnel face reception. The results are as follows: Figure 8 and Fig. 9 shown. Figure 8 and Fig. 9 There are obvious similarities, especially the Y component. Figure 5 , Fig. 9 The resolution of the reflected waves in the tunnel is significantly improved, and the surface wave interference is small. The numerical simulation results show that the resolution of the wave field record after three-dimensional reverse time extension is close to that of the wave field record received at the tunnel face, and is better than the resolution of the wave field record received at the tunnel wall.

[0089] The above invention content is currently mainly used in the field of tunnel earthquake advance prediction.

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

Claims

1. A three-dimensional inverse time extension method for seismic waves, characterized in that: The steps include: Step 1: Pick up the first arrival travel time information of the seismic data of tunnel wall excitation and tunnel wall reception, and calculate the longitudinal wave velocity Vp and shear wave velocity Vs of the surrounding rock in combination with the spacing of the shot points; Step 2: Using the calculated P-wave velocity Vp and S-wave velocity Vs as the velocity model, three-dimensional reverse time extension calculation is performed from the maximum time to the minimum time. During the calculation process, virtual geophones are arranged on the tunnel face to receive the wave field of three-dimensional reverse time extension. Step 3: Perform autocorrelation imaging on the wave field after three-dimensional reverse time extension, and further convert the seismic wave field of tunnel wall excitation and tunnel face reception into the seismic wave field of tunnel face excitation and tunnel face reception.

2. The three-dimensional inverse time extension method of seismic waves according to claim 1, characterized in that: In step 2, the derivation formula of three-dimensional inverse time extension is: The calculation formula of the decoupled first-order partial differential equation in the three-dimensional elastic wave equation is as follows: By discretizing the differential equations (1) to (9), we obtain the differential formula for three-dimensional inverse time extension: For velocity components: In the above formulas, i, j, k are subscripts of the X, Y and Z spatial directions; n is the superscript of time; Δx, Δy and Δz are the grid spacings in 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 to the horizontal direction, Y is the forward excavation direction of the tunnel, and Z is parallel to the tunnel face and points to the vertical downward direction; Formulas (11) to (18) are second-order accuracy in space and time, and can be derived to 2L-order in space by the same logic, 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 is the stress component in the X direction, τ yy is the stress component in the Y direction, τ zz is 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 is the stress component of the resultant force in the X and Z directions; (f x ,f y ,f z ) is the body force source component in the X, Y and Z directions, is the three-component seismic data received by the tunnel wall during three-dimensional reverse time extension, ρ is the density of the medium; λ and μ are the first and second parameters of the Lamé constant, respectively, and their relationship with the longitudinal wave velocity Vp and the shear wave velocity Vs is: λ=ρ(Vp 2 -2Vs 2 ) (19); μ=ρVs 2 (20); The calculation formulas (9) to (18) are used to perform reverse time delay calculation from the maximum time to the minimum time. During the calculation process, virtual detectors are arranged on the tunnel face to receive the reverse time delayed wave field.

3. The three-dimensional inverse time extension method of seismic waves according to claim 1, characterized in that: The method further includes step 4, which implements bandpass filtering, automatic gain, predictive deconvolution, full waveform inversion, and reverse time migration data processing based on the processing result of step 3.

4. A high-resolution tunnel earthquake advance prediction method, comprising the three-dimensional reverse time extension method of seismic waves according to any one of claims 1 to 3, characterized in that: A virtual geophone is set at the tunnel face position, and the tunnel wall excitation-tunnel wall reception observation system is calculated through three-dimensional reverse time extension, and the wave field received by the tunnel wall is reversely time-recursed to the position of the virtual geophone at the tunnel face, and finally converted into a tunnel face excitation-tunnel face reception observation system; The specific data processing steps are: Step A: Data preprocessing: including data editing, bandpass filtering, FK filtering, and automatic gain; Step B: Data imaging: three-dimensional reverse time extension, predictive deconvolution, full waveform inversion and two-dimensional and three-dimensional reverse time migration imaging, wherein the three-dimensional reverse time extension adopts the three-dimensional reverse time extension method of seismic waves described in any one of claims 1 to 3 to process the data; Step C: Geological interpretation: Combine the full waveform inversion and reverse time migration imaging results to extract unfavorable geological bodies, classify surrounding rock grades, and provide support recommendations.

Images