Tunnel seismic data linear moveout correction and phase delay focused reverse time migration imaging method

Abstract

This invention discloses a linear dynamic correction and phase delay focusing reverse time migration imaging method for tunnel seismic data, comprising the following steps: Step 1: First, the linear dynamic correction method is used to eliminate the influence of offset, converting single-shot tunnel seismic data into self-excited and self-received data, and the corresponding shot point positions are also converted to the new positions; Step 2: During forward and reverse time extension, the wavelet and self-excited and self-received data are sequentially given different phase delays according to the location of the record trace to perform reverse time migration imaging; Step 3: These reverse time migration imaging results are weighted and superimposed using signal-to-noise ratio statistics to obtain the final focused reverse time migration imaging result. This invention aims to overcome the influence of small offsets, better identify and distinguish the types and spatial structural characteristics of adverse geological bodies, and improve the accuracy of reverse time migration imaging.

Description

Technical Field

[0001] This invention relates to the field of tunnel seismic detection and engineering survey technology, and in particular to a method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data. Background Technology

[0002] Tunnel seismic advance detection is one of the most effective geophysical methods for predicting long-distance geological hazards within tunnels. It can accurately predict geological anomalies 100m to 150m in the advance detection direction, including fault fracture zones, karst caves, weak rock layers, and water-rich zones (Andisheh et al., 2008; Shi et al., 2014). Migration imaging is one of the core technologies for tunnel seismic data processing. The narrow tunnel space limits the acquisition of tunnel seismic data, resulting in small offset distances, multiple components, and small data volume. Many scholars have studied tunnel seismic migration imaging methods. Ashida et al. (2001) proposed the equal-travel-time surface migration method based on the scattering superposition theory, which can achieve rapid imaging, but the imaging accuracy is low. Luth et al. (2005) proposed restricting the Kirchhoff migration operator within the Fresnel zone volume of the reflection path to improve the imaging accuracy of tunnel seismic data with small migration apertures. Compared to the Kirchhoff migration method (Lüthet et al., 2008; Rechlin et al., 2009; Tzavaras et al., 2012; Bellino et al., 2013), the reverse time migration (RTM) method utilizes the cross-correlation principle of forward and reverse time continuation for imaging. The imaging area is not limited by the migration aperture, making it the migration method with the highest imaging accuracy to date, especially suitable for migration imaging of short-offset tunnel seismic data. However, there are relatively few cases of RTM application to tunnel seismic data; it is mostly used in surface seismic data processing.

[0003] Many scholars have studied reverse time migration (RTM) methods for the ground, achieving breakthroughs in imaging theory, imaging conditions, and computational speed, and obtaining good results in practical applications. Chang et al. (1986) proposed excitation time imaging conditions for RTM imaging. Mulder et al. (2004) found that iterative migration can reduce migration artifacts in RTM. Yan et al. (2008) extended the cross-correlation imaging conditions and proposed elastic wave images for angular domain analysis. Fletcher et al. (2009) realized acoustic RTM migration of tilted transverse isotropy (TTI) media. Liu et al. (2015) improved the accuracy of RTM imaging with multiple waves, which can enhance the illumination of the target and reduce the influence of false anomalies. Silvestrov et al. (2016) realized post-stack RTM imaging. High-performance computing methods such as GPU (Liu et al., 2010) and FPGA (Medeiros et al., 2013) have also been used to improve the computational efficiency of RTM. Zha et al. (2018) achieved reverse time migration imaging of the two-dimensional elastic wave equation of the tunnel velocity model. Zou et al. (2021) proposed a passive source reverse time migration method for imaging three-dimensional complex structures.

[0004] Reverse-time migration methods have evolved from two-dimensional to three-dimensional, but the problem of small offset distances remains unsolved. Combining the principles of seismic wave propagation in tunnel excavation sites and within tunnel spaces, this invention proposes a linear dynamic correction and phase delay focusing reverse-time migration method for tunnel seismic data. This method aims to focus the reverse-time migration imaging range, suppress low-frequency noise interference, improve migration imaging accuracy, and ultimately enhance the accuracy of tunnel seismic advance prediction.

[0005] References

[0006] Andisheh,A.,Ali,M.,Reza,N,Mojtaba,ZS,Afshin,E,2008,Prediction ofgeological hazardous zones in front of a tunnel face using TSP-203 andartificial neural networks.Tunnelling and Underground Space Technology,23:711-717.

[0007] Ashida Y.2001.Seismic imaging ahead of a tunnel face with three-component geophones.International Journal of Rock Mechanics and MiningSciences,38:823-831.

[0008] Bellino,A.,Garibaldi,L.,Godio,A.,2013.An automatic method for dataprocessing of seismic in tunneling.Journal of Applied Geophysics,98:243-253.

[0009] Chang,W.F.,McMechan,G.A.,1986.Reverse-time migration of offsetvertical seismic profiling data using the excitation-time imagingcondition.Geophysics,51(1):67-84.

[0010] Fletcher,R.P.,Du,X.,Fowler,P.J.,2009.Reverse time migration in tiltedtransversely isotropic(TTI)media.Geophysics,74(6):Wca179-Wca187.

[0011] Liu,H.W.,Li,B.,Liu,H.,Tong,X.L.,Liu,Q.,2010.The algorithm of highorder finite difference pre-stack reverse time migration and GPUimplementation.Chin.J.Geophys,53(7):1725-1733.

[0012] Liu,Y.K.,Hu,H.,Xie,X.B.,Zheng,Y.C.,Li,P.,2015.Reverse time migrationof internal multiples for subsalt imaging.Geophysics,80(5):S175-S185.

[0013] Luth S.,Buske S.,Giese R.,Goertz A.2005.Fresnel volume migration ofmulticomponent data.Geophysics,70:121-129.

[0014] Lüth,S.,Giese,R.,Otto,P.,Krüger,K.,Mielitz,S.,Bohlen,T.,Dickmann,T.,2008.Seismic investigations of the Piora Basin using S-wave conversions atthe tunnel face of the Piora adit(Gotthard Base Tunnel).Int.J.RockMech.Min.Sci.45(1):86-93.

[0015] Medeiros,V.,Barros,A.,Silva-Filho,A.,de Lima,M.E.2013.Highperformance implementation of RTM seismic modeling on FPGAs:architecture,arithmetic and power issues.In:High-performance Computing UsingFPGAs.Springer,pp.305-334.

[0016] Mulder,W.A.,Plessix,R.E.,2004.A comparison between one-way and two-way wave equation migration.Geophysics,69(6):1491-1504.

[0017] Rechlin,A.J.,Lüth,S.,Giese,R.,2009.OnSITE:integrated seismic imagingand interpretation for tunnel excavation.Proceedings of the InternationalConference on Rock Joints and Jointed Rock Masses,pp.1-7.

[0018] Shi,S.S.,Li,S.C.,Li,L.P.,Zhou,Z.Q.,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.

[0019] Silvestrov,I.,Baina,R.,Landa,E.,2016.Poststack diffraction imagingusing reverse-time migration.Geophysical Prospecting,64:129-142.

[0020] Tzavaras,J.,Buske,S.,Groβ,K.,Shapiro,S.,2012.Three-dimensionalseismic imaging of tunnels.Int.J.Rock Mech.Min.Sci.49:12-20.

[0021] Yan,J.,Sava,P.,2008.Isotropic angle-domain elastic reverse-timemigration.Geophysics,73(6):S229-S239.

[0022] Zha,XJ,Gao,X.,Wang,W.,et al.,2018.Advanced prediction migrationmethod research in tunnel engineering investigation.Chinese Journal ofGeophysics(in Chinese),61(3):1150-1156.

[0023] Zou,P.,Cheng,JB,2021.Three-dimensional passive-source anisotropicreverse time migration for imaging lithospheric discontinuities.Geophys.J.Int.226:2103-2115. Summary of the Invention

[0024] In view of at least one deficiency of the prior art, the purpose of the present invention is to provide a method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data, which aims to overcome the influence of small offset distance, better identify and distinguish the types and spatial structural features of adverse geological bodies, and improve the accuracy of reverse time migration imaging.

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

[0026] A linear dynamic correction and phase delay focusing reverse time migration imaging method for tunnel seismic data, the key of which includes the following steps:

[0027] Step 1: Use the Tunnel Seismic Prediction (TSP) observation system to collect seismic data, pick out the first arrival and travel time of the seismic data, and calculate the average velocity v of the surrounding rock using the least squares fitting method;

[0028] Step 2: Preprocess tunnel seismic data and extract reflected wave records (RC) in the leading direction;

[0029] Step 3: Based on Step 2, using velocity v as the linear dynamic correction velocity, linear dynamic correction is applied to the reflected wave record RC to eliminate the influence of offset, and the single shot record is corrected to a self-excited and self-receiving seismic record RE with zero offset.

[0030] Step 4: Based on Step 3, convert the self-excited and self-received seismic records (RE) into seismic records (REA) with different phase delays;

[0031] Step 5: Based on Step 4, perform phase-delay focused reverse-time migration imaging. During the forward continuation calculation, the Ricker wavelet phase delay corresponds to a given delay. During the reverse-time continuation calculation, the input wavefield is the phase-delay corrected seismic record REA. The focused reverse-time migration imaging results are obtained by applying cross-correlation imaging conditions.

[0032] Step Six: Repeat Step Five, successively giving phase delays Δt at different intervals. i (i = 1, 2, ..., m, where m can be considered as the number of scan stackings) Focusing reverse time migration is performed separately, and these migration imaging results are weighted and superimposed with signal-to-noise ratio statistics to obtain the final phase-delayed focusing reverse time migration imaging result.

[0033] In step two, bandpass filtering, automatic gain control, geometric diffusion correction, FK filtering, and adaptive polarization filtering are used sequentially to preprocess the tunnel seismic data.

[0034] In step three, the formula for calculating the time of linear dynamic correction is:

[0035]

[0036] This is the time after linear dynamic correction, t0 is the time of the original tunnel seismic record, and offset is... k is the distance from the detector to the shot point, and v is the velocity of the linear dynamic correction, which is also the velocity obtained by fitting the first arrival wave travel time in step one. After calculating the linear dynamic correction time, the single-shot record RC is corrected to a self-excited and self-received wavefield record RE. The calculation formula is:

[0037]

[0038] x1 k It records the actual position of each track of the RC gun, x k These are the actual locations of each trace in the self-excited and self-received wavefield record (RE), k = 1, 2, ..., N, where N is the trace number in the tunnel seismic record. The purpose of linear dynamic correction is to eliminate the influence of offset.

[0039] Step four includes assuming there are N tunnel seismic records, and giving equally spaced delays Δt, with the first record having a delay of 0, the second record having a delay of Δt, and the Nth record having a delay of (N-1)×Δt, and giving different delay values ​​Δt for each. i (i = 1, 2, ..., m), corrected to the corresponding tunnel seismic record REA i (i = 1, 2, ..., m); m is the number of scan stacking times.

[0040] In step four, the derivation formula for the phase delay method is as follows:

[0041] REA(x 2k ,t+(k-1)Δt i )=RE(x k ,t) (3);

[0042] x 2k The actual position of the gather after linear dynamic correction, where t is time, and Δt i This is for phase delay.

[0043] In step five, the correction formula for the wavelet phase delay during the forward extension process in the phase delay focusing reverse time offset is:

[0044] CRW(x 2k ,t+(k-1)Δt i )=RW(t) (4);

[0045] RW stands for Reck wavelet, and CRW stands for wavelet after phase delay correction.

[0046] In the forward and reverse time continuation calculations of reverse time migration, the corrected wavelet CRW and the tunnel seismic record REA were used as input wavefields to begin the iterative calculations. The reverse time migration imaging results were obtained by applying cross-correlation imaging conditions, and the calculation formula is as follows:

[0047]

[0048] Scontinuation(x,z,t) is the forward continuation wave field at spatial location (x,z) at different times t. Rcontinuation(x,z,t) is the reverse continuation wave field at spatial location (x,z) at different times t. RESULT(x,z) is the reverse time migration result.

[0049] In step six, the formula for calculating the weighted superposition of the focusing reverse-time offset based on the signal-to-noise ratio statistics is as follows:

[0050]

[0051] S is the number of samples per channel in the reverse time offset data, w i (j,k) represents the weights, and their values ​​are positively correlated with the signal-to-noise ratio (SNR) statistics of the reverse-time offset data. Regions with high SNR are assigned larger weights, while regions with low SNR are assigned smaller weights. The sum of the weights is m, and the specific calculation formula is as follows:

[0052]

[0053] Significant effects: This invention provides a method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data, which can better overcome the influence of small offset distances, better identify and distinguish the types and spatial structural features of adverse geological bodies, and improve the accuracy of reverse time migration imaging. Attached Figure Description

[0054] Figure 1 Cave model;

[0055] Figure 2 (a) Full-wave field record; (b) TSP field record; (c) TSP field record with linear dynamic correction.

[0056] Figure 3 Wavefield snapshot at 7ms: (a) Forward continuation profile of conventional reverse time migration method; (b) Forward continuation profile of focused reverse time migration method; (c) Reverse time continuation profile of conventional reverse time migration method; (d) Reverse time continuation profile of focused reverse time migration method;

[0057] Figure 4 Comparison of reverse time migration results for cave models: (a) conventional reverse time migration method; (b) focused reverse time migration method;

[0058] Figure 5 This is a flowchart of the method of the present invention. Detailed Implementation

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

[0060] This invention relates to a linear dynamic correction and phase delay focusing reverse time migration imaging method for tunnel seismic data, which can better overcome the influence of small offset distances, better identify and distinguish the types and spatial structural features of adverse geological bodies, and improve the accuracy of reverse time migration imaging.

[0061] like Figures 1-5 As shown, the present invention adopts the following technical solution:

[0062] 1. A linear dynamic correction and phase delay focusing reverse time migration imaging method for tunnel seismic data, specifically, firstly, the linear dynamic correction method is used to eliminate the influence of the offset, converting the single-shot tunnel seismic data into self-excited and self-received data, and the corresponding shot point positions are also converted to the new positions; then, the wavelet and self-excited and self-received data are given different phase delays according to the position of the record trace in order to perform reverse time migration imaging; finally, these reverse time migration imaging results are weighted and superimposed with signal-to-noise ratio statistics to obtain the final focused reverse time migration imaging result.

[0063] The formula for calculating the time of linear dynamic correction is:

[0064]

[0065] This is the time after linear dynamic correction, t0 is the time of the original tunnel seismic record, and offset is... k is the distance from the detector to the shot point, and v is the velocity of the linear dynamic correction, which is also the velocity obtained by fitting the first arrival wave travel time in step one. After calculating the linear dynamic correction time, the single-shot record RC is corrected to a self-excited and self-received wavefield record RE. The calculation formula is:

[0066]

[0067] x1 k It records the actual position of each track of the RC gun, x k These are the actual locations of each trace in the self-excited and self-received wavefield record (RE), k = 1, 2, ..., N, where N is the trace number in the tunnel seismic record. The purpose of linear dynamic correction is to eliminate the influence of offset.

[0068] The derivation formula for the phase delay method is as follows:

[0069] REA(x 2k ,t+(k-1)Δt i )=RE(x k ,t) (3);

[0070] x 2k The actual position of the gather after linear dynamic correction, where t is time, and Δt i For phase delay, i = 1, 2, ..., m, where m is the number of scan stacking.

[0071] The formula for correcting the wavelet phase delay during the forward extension process in phase-delay focusing reverse-time migration is:

[0072] CRW(x 2k ,t+(k-1)Δt i )=RW(t) (4);

[0073] RW stands for Reck wavelet, and CRW stands for wavelet after phase delay correction.

[0074] In the forward and reverse time continuation calculations of reverse time migration, the corrected wavelet CRW and the tunnel seismic record REA were used as input wavefields to begin the iterative calculations. The reverse time migration imaging results were obtained by applying cross-correlation imaging conditions, and the calculation formula is as follows:

[0075]

[0076] Scontinuation(x,z,t) is the forward continuation wave field at spatial location (x,z) at different times t. Rcontinuation(x,z,t) is the reverse continuation wave field at spatial location (x,z) at different times t. RESULT(x,z) is the reverse time migration result.

[0077] The formula for calculating the weighted superposition of focusing reverse-time offset based on signal-to-noise ratio statistics is as follows:

[0078]

[0079] m is the number of scan stacks, N is the number of channels, and S is the number of samples per channel in the reverse time offset data. In formula (6), This is equivalent to normalizing to between 0 and 1, w i (j,k) represents the weights, whose values ​​are positively correlated with the signal-to-noise ratio (SNR) statistics of the reverse-time offset data. Regions with high SNR are assigned larger weights, while regions with low SNR are assigned smaller weights. The sum of the weights is m, and the specific calculation formula is as follows:

[0080]

[0081] The specific advantages of this invention are:

[0082] (1) The velocity accuracy used in linear dynamic correction is high, and linear dynamic correction processing of tunnel seismic data can completely eliminate the influence of offset.

[0083] (2) Focused reverse time migration imaging of adverse geological bodies can be achieved without increasing the data acquisition cost at the tunnel construction site;

[0084] (3) The focused reverse time migration method can better overcome the influence of small offset distance, better identify and distinguish the types and spatial structure characteristics of undesirable geological bodies, and improve the detection accuracy of undesirable geological bodies.

[0085] This invention simplifies the propagation process of elastic waves in three-dimensional tunnel space into a two-dimensional model, and simplifies elastic waves into sound waves. It creatively proposes a phase delay focusing reverse time migration imaging method, which improves the detection accuracy of adverse geological bodies in the leading direction by changing the delay of the wave field during forward and reverse time extension and focusing the reverse time migration imaging results.

[0086] The phase-delayed focusing reverse time-shifting imaging method is within the scope of patent protection. The phase-delayed focusing reverse time-shifting imaging method is not limited to the two-dimensional application described in the patent specification; the three-dimensional phase-delayed focusing reverse time-shifting imaging method is also within the scope of patent protection.

[0087] Model Experiment

[0088] 1. Tunnel seismic data acquisition

[0089] To verify the method of this invention, a two-dimensional cave model was designed, such as... Figure 1 As shown. Numerical simulation and reverse time migration are achieved using two-dimensional acoustic wave equations. G1 and G2 are two detector positions, and shots are 24 shot points. The P-wave velocity Vp of the surrounding rock is 5400 m / s, and the density Rho is 2.12 g / cm³. 3 The cave has a diameter of 36m, a longitudinal wave velocity of 2500m / s, and a density of 1.4g / cm³. 3 The data sampling interval is 35 μs, and the number of sampling points per channel is 4000. Simulated 24-shot data is used, with acoustic signals received at positions G1 and G2. The acoustic signal received by the G1 detector is as follows... Figure 2 As shown. Symbol D1 represents the direct wave, R1 and R2 represent the reflected waves from interfaces I and II, and M1 to M5 represent multiple waves. Figure 2 (a) in the figure is a full-range wavefield record with a lateral range of 0m to 200m. The Y position of the receiver is the same as that of G1. Figure 2 (b) shows the extracted 24-shot TSP seismic record. Linear dynamic correction was applied at a velocity of 5400 m / s, resulting in the following... Figure 2 (c) The corrected record.

[0090] To visually demonstrate the effect of the focused reverse-time offset, wavefield snapshots of the forward and reverse-time continuations at 7ms are extracted, as follows: Figure 3 As shown. Figure 3 (a) and Figure 3 (b) shows the wavefield snapshots at 7ms forward extension using the conventional reverse time migration and focused reverse time migration methods, respectively. P1 to P3 are the wavefronts of the reflected waves. P1 is approximately semi-circular, while the energy of P2 and P3 begins to weaken, and P4 is strengthened. Figure 3 (c) and Figure 3 (d) represents wavefield snapshots taken at 7 ms in both conventional and focused reverse-time migration methods, respectively. P5 to P6 are the wavefronts of the reflected waves, compared to... Figure 3 (c) in the middle Figure 3 In (d), the energies of P5 and P6 are stronger and more pronounced.

[0091] The conventional reverse time migration method and the focused reverse time migration method were used to obtain the following results: Figure 4 The offset imaging results show that I-1 and I-2 are strongly reflective interfaces of interface I, and II-1 is a strongly reflective interface of interface II. N1 to N4 represent artifact noise. Compared to Figure 4 (a) in the middle, Figure 4 In (b), the imaging of cave interface I is more focused, making it easier to identify the structural features of the cave. In particular, I-3 and I-4 can clearly indicate the interface outline of the cave. Figure 4In (b), the artifact noise energy is also weaker, especially in N4.

[0092] Numerical simulation results demonstrate that focused reverse time migration imaging can more clearly display the structural features of caves, suppress artifact noise interference, and improve the resolution of cave anomalies.

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

[0094] 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 method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data, characterized in that, Includes the following steps: Step 1: Use the tunnel earthquake prediction observation system to collect earthquake data, pick up the first arrival and travel time of the earthquake data, and calculate the average velocity v of the surrounding rock using the least squares fitting method; Step 2: Preprocess tunnel seismic data and extract reflected wave records (RC) in the leading direction; Step 3: Using velocity v as the linear dynamic correction velocity, perform linear dynamic correction on the reflected wave record RC to eliminate the influence of offset, and correct the single shot record to a self-excited and self-receiving seismic record RE with zero offset; Step 4: Convert the self-excited and self-received seismic records (RE) into seismic records (REA) with different phase delays; Step 5: Perform phase-delay focused reverse-time migration imaging. In the forward continuation calculation, the phase delay of the Ricker wavelet corresponds to the given delay. In the reverse-time continuation calculation, the input wavefield is the phase-delay corrected seismic record REA. The focused reverse-time migration imaging results are obtained by applying the cross-correlation imaging conditions. Step Six: Repeat Step Five, successively giving phase delays Δt at different intervals. i Each phase-delayed focusing reverse time migration is performed separately, and the migration imaging results are weighted and superimposed with signal-to-noise ratio statistics to obtain the final phase-delayed focusing reverse time migration imaging result.

2. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 1, characterized in that: In step two, bandpass filtering, automatic gain control, geometric diffusion correction, FK filtering, and adaptive polarization filtering are used sequentially to preprocess the tunnel seismic data.

3. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 1, characterized in that: In step three, the formula for calculating the time of linear dynamic correction is: This is the time after linear dynamic correction, t0 is the time of the original tunnel seismic record, and offset is... k is the distance from the detector to the shot point, v is the velocity of the linear dynamic correction, which is also the velocity obtained by fitting the first arrival wave travel time in step one; after calculating the linear dynamic correction time, the single-shot record RC is corrected to a self-excited and self-received wavefield record RE, and the calculation formula is: x1 k It is a single-shot record of the actual position of each track of the RC, x k It represents the actual location of each trace in the self-excited and self-received wavefield record (RE), k = 1, 2, ..., N, where N is the trace number of the tunnel seismic record.

4. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 3, characterized in that: Step four includes assuming there are N tunnel seismic records, and giving equally spaced delays Δt, with the first record having a delay of 0, the second record having a delay of Δt, and the Nth record having a delay of (N-1)×Δt, and giving different delay values ​​Δt for each. i (i = 1, 2, ..., m, where m is the number of scan stacks), corrected to the corresponding tunnel seismic record REA i (i = 1, 2, ..., m).

5. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 4, characterized in that: In step four, the derivation formula for the phase delay method is as follows: REA(x 2k ,t+(k-1)Δt i )=RE(x k ,t) (3); x 2k t represents the actual position of the gather after linear dynamic correction, and t represents time.

6. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 5, characterized in that: In step five, the correction formula for the wavelet phase delay during the forward extension process in the phase delay focusing reverse time offset is: CRW (x 2k ,t+(k-1)Δt i )=RW(t) (4); RW represents the Ricker wavelet, and CRW represents the phase-delay-corrected wavelet. In the forward and reverse-time extrapolation calculations of reverse-time migration, the corrected wavelet CRW and the tunnel seismic record REA are used as the input wavefields to begin iterative calculations. The reverse-time migration imaging results are obtained by applying cross-correlation imaging conditions, and the calculation formula is as follows: Scontinuation(x,z,t) is the forward continuation wave field at spatial location (x,z) at different times t; Rcontinuation(x,z,t) is the reverse continuation wave field at spatial location (x,z) at different times t; RESULT(x,z) is the reverse time migration result.

7. The method for linear dynamic correction and phase delay focusing reverse time migration imaging of tunnel seismic data according to claim 6, characterized in that: In step six, the formula for calculating the weighted superposition of the focusing reverse-time offset based on the signal-to-noise ratio statistics is as follows: S is the number of samples per channel in the reverse time offset data, w i (j,k) represents the weights, and the sum of the weights is m. The specific calculation formula is as follows:

Images