Pre-stack elastic reverse time migration in tunnels based on cylindrical coordinates

Yuxio Ren,Jinsen Wng,Zhicho Yng,Xinji Xu,*,Lei Chen,**

a Geotechnical and Structural Engineering Research Center,Shandong University,Jinan,250061,China

b School of Qilu Transportation,Shandong University,Jinan,250061,China

Keywords:Tunnel geological forward-prospecting Elastic reverse time migration Cylindrical coordinates Numerical singularity Field application

ABSTRACT Seismic forward-prospecting in tunnels is an important step to ensure excavation safety.Nowadays,most advanced imaging techniques in seismic exploration involve calculating the solution of elastic wave equation in a certain coordinate system.However,considering the cylindrical geometry of common tunnel body,Cartesian coordinate system seemingly has limited applicability in tunnel seismic forwardprospecting.To accurately simulate the seismic signal received in tunnels,previous imaging method using decoupled non-conversion elastic wave equation is extended from Cartesian coordinates to cylindrical coordinates.The proposed method preserves the general finite-difference time-domain (FDTD)scheme in Cartesian coordinates,except for a novel wavefield calculation strategy addressing the singularity issue inherited at the cylindrical axis.Moreover,the procedure of cylindrical elastic reverse time migration(CERTM)in tunnels is introduced based on the decoupled non-conversion elastic wavefield.Its imaging effect is further validated via numerical experiments on typical tunnel models.As indicated in the synthetic examples,both the PP-and SS-images could clearly show the geological structure in front of the tunnel face without obvious crosstalk artifacts.Migration imaging using PP-waves can present satisfactory results with higher resolution information supplemented by the SS-images.The potential of applying the proposed method in real-world cases is demonstrated in a water diversion tunnel.In the end,we share our insights regarding the singularity treatment and further improvement of the proposed method.

1.Introduction

The prospecting of the unknown geology prior to excavation is of vital importance for tunneling projects,especially for mechanized tunneling using tunnel boring machines (TBMs).The lack of knowledge regarding adverse geology like faults,fractured zones or other lithologic and structural heterogeneities beforehand may bring danger to the tunneling activities as well as the construction personnel.However,common geological surveys on ground surface can only provide a coarse exploration of the subsurface geology.To meet the requirement of accurate detection in tunnels,seismic methods have been widely used due to its long detection distance and the interface recognition effect (Li et al.,2017).

With the increase of tunneling difficulty,researchers and engineers have developed various seismic prospecting systems,including tunnel seismic prediction (Sattel et al.,1992; Alimoradi et al.,2008),true reflection tomography (Otto et al.,2002;Yamamoto et al.,2011),and integrated seismic imaging system(Bohlen et al.,2003).These methods mainly consider the seismic wave reflected from the geological structures ahead and achieve the locating and imaging of the adverse geology based on ray theory.Besides reflection waves,seismic prospecting methods using shear and surface waves have been developed to provide an improved detection image(Bohlen et al.,2007;Jetschny et al.,2011;Bharadwaj et al.,2017).Moreover,with the advance in computation powers,migration techniques used in tunnels have been promoted from ray-based methods like Kirchhoff method (Bleistein,1987;Ashida,2001;Lüth et al.,2008;Tzavaras et al.,2012)to wave-based methods like reverse time migration(RTM),which is considered to be a more accurate tool to image the complex subsurface reflectors(Baysal et al.,1983;McMechan,1983;Whitmore,1983;Zhou et al.,2018).

Specifically for RTM in tunnels,most researchers simplify the prospecting situation to an acoustic model and migrate the seismic wavefield based on the velocity of surrounding rock.Chang et al.(2006) made a preliminary attempt in numerical experiments where the tunnel body is neglected given a velocity the same as the surrounding rock.Similar strategy was also verified by Gong et al.(2010) in field data and the advantages of using RTM over Kirchhoff integral migration are demonstrated in various aspects.Recently,by taking the effect of tunnel body in consideration,Li et al.(2019) incorporated the beamforming method into RTM imaging condition and achieved an improved imaging effect in both numerical and field experiments.With the elastic RTM (ERTM) showing superior imaging results in surface seismic survey (Sun and McMechan,2001; Zhang et al.,2007;Zhang and McMechan,2010; Tang and McMechan,2018),its application in tunnels begin to capture the interest of scholars.One major obstacle on the way is the crosstalk artifacts caused by converted waves on tunnel surfaces.To mitigate this issue,Cheng et al.(2014) proposed a two-dimensional (2D) ERTM method on the basis of the non-conversion wave equation and achieved a good application effect in tunnel data.Moreover,by carefully analyzing the practical influence of the tunnel body,Ren et al.(2020) proposed a novel ERTM method based on the decoupled non-conversion elastic equation.In this way,P-and S-wavefield can be propagated independently without generating converted waves and satisfactory imaging results can be obtained in practical tunneling projects.

However,current study of tunnel based RTM is mostly in Cartesian coordinate system,while in practical cases,tunnel cross-section is usually considered to be circular and the seismic wavefield studied in Cartesian coordinates may fail to fully reveal the true feature of practical seismic wavefield.For example,tunnel body is usually approximately simplified as a cylinder in Cartesian coordinates,the corresponding wavefield would be seriously influenced by the edges and corners of the cubic grid,leading to diffraction waves on tunnel sidewalls.Thus,in order to avoid the diffraction interference,cylinder is regarded as a better approximation of the tunnel body and cylindrical coordinate system is recommended in the study of seismic prospecting in tunnels.For example,Liu et al.(2017) has carried out forward modeling of tunnel environment using polar coordinates.Jetschny et al.(2010) simulated the surface wavefield on the side wall of the tunnel through cylindrical coordinates and observed the propagation characteristics of tunnel surface wave.However,there is little research on using cylindrical coordinates to simulate the full elastic wavefield in tunnel geological forwardprospecting.

As for the elastic modeling in cylindrical coordinates,Kessler and Kosloff (1991) proposed a pseudo-spectral method to address the problem of 2D elastic wave propagation in polar coordinates.Liu and Sinha (2003) derived the three-dimensional(3D) elastic wave equation in cylindrical coordinates and achieved good modeling results in a pressurized borehole using an improved cylindrical perfectly matched layer as the absorbing boundary condition.To mitigate the numerical instability caused by the increasing polar radius and azimuth inflation,Pissarenko et al.(2010) developed a finite-difference method with periodical refinement of the grid step in the azimuth direction.Another numerical instability issue of seismic wave simulation in cylindrical coordinates lies on the vicinity of the cylindrical axis.The term 1/r in the finite-difference scheme will become a singularity at the radius r=0 (Chen et al.,1998; Cheng and Blanch,2008;He et al.,2013).Unfortunately,most seismic modeling methods in cylindrical coordinates are developed for borehole situation,where the vicinity of the cylindrical axis can be neglected in the seismic wavefield computation.While in tunnels,geological condition along the tunnel axis is the major target area for seismic prospecting and it will greatly affect the tunnel excavation process.Thus,in order to take full advantage of the cylindrical coordinate system in tunnel geological forwardprospecting,it is necessary to study the method of accurate seismic wavefield simulation and imaging in cylindrical coordinates,especially near the axis of cylinder.

In this work,we extend our previous study in Cartesian coordinate system (Ren et al.,2020) to cylindrical coordinate system.After the derivation of the decoupled non-conversion elastic wave equation in cylindrical coordinates,the proposed strategy is followed to deal with the singularity at the axis of the cylinder.Then,procedure of cylindrical ERTM (CERTM) in tunnels is proposed and its application effect is verified via numerical examples and field cases.Insights on the proposed approximate strategy handling the wavefield singularity issue are further discussed in the end.

2.Theory

2.1.Decoupled non-conversion elastic wave equation

As analyzed in our previous study (Ren et al.,2020),a good ERTM image in tunnel geological forward-prospecting usually relies on the quality and independence of reflective signal in seismic record.Specifically,to mitigate the crosstalk artifacts in ERTM,wavefield conversion at the tunnel surface should be avoided and P-and S-wavefield separation via decoupled wave equation or field data separation methods is required.Thus,in this work,we continue focusing on the decoupled nonconversion elastic wave equation and extend its application in cylindrical ERTM to improve the imaging effect of forwardprospecting in tunnels.

Starting from the non-conversion elastic wave equation (Yao,1996) in elastic isotropic medium with constant density,we have

where u denotes the particle displacement of elastic wavefield;and Vpand Vsare the P-and S-wave velocities,respectively.According to the Helmholtz theory,this equation can be rewritten as the summation of P-wave wavefield upand S-wave wavefield us:

Moreover,by introducing the stress field τp=V2p?·u and τs=-V2s?×u and defining the particle vibration velocity v=?u/?t,we can derive the first-order velocity-stress equation for the decoupled non-conversion elastic wave in cylinder coordinates(r,θ,z) as

where vi(i={r,θ,z}) represents the particle velocity in the i direction,τpis the compressional stress,and τsij(i,j={r,θ,z}) denotes the shear stress.Moreover,similarly to the Cartesian coordinates,the main difficulty of separating P-and S-waves in cylindrical coordinates is similar to calculate the divergence and curl of the wavefield,which involves some extra calculation of rθ-related terms.

To numerically solve Eq.(3),a staggered-grid finite-difference time-domain (FDTD) method (Yee,1966; Chen et al.,1998;Pitarka,1999) is applied in this work.As shown in Fig.1,the compressional stress is placed at the center of the staggered grid cell,shear stresses are located at the cell edges and the particle velocity components are sampled at the center of the cell faces.To mitigate the numerical dispersion resulting from the increasing polar radius and azimuth inflation,the grid refinement strategy proposed by Pissarenko et al.(2010) is applied.The detailed variable interpolation and update formulae are exemplified in Appendix A.In addition,perfectly matched layers (Liu and Sinha,2003) are applied as absorbing boundary condition to prevent artificial reflections from the exterior boundary of the cylindrical girds.

Fig.1.A unit staggered grid cell with discretized velocity and stress components.

2.2.Wavefield calculation at the cylindrical axis

As mentioned above,when numerically calculating the wavefield propagation in cylindrical coordinates using the FDTD method,the presence of the term 1/r in the wave equation will cause a mathematical singularity at the cylindrical axis,i.e.the radius r=0.However,this singularity issue arises only because of the choice of coordinate system and the physical wavefield itself does not possess any singular point.Considering the different physical natures of P-and S-wavefields,different singularity treatments are proposed for the decoupled P-and Swaves,respectively.

According to Eq.(2),the non-conversion elastic wavefield can be separated into a scalar P-wavefield and a vectorial S-wavefield,which is in fact a vortex field.Accordingly,the wave equation variables placed in the first cell from the center can be separated as Fig.2.One can easily see that the singularity issue will only affect the variable vprin the P-wavefield and variables vsr,τsrθand τszrin the S-wavefield.Thus,we will address the strategy of calculating these variables with respect to different wavefields.

Fig.2.A staggered grid cell at the cylindrical axis defined for (a) the decoupled Pwavefield and (b) the decoupled S-wavefield.The singularity variables on the cylindrical axis are marked in red.

For P-waves,the finite-difference scheme for calculating vprat r=0 should be

However,there is no τp(-1/2,j+1/2,k+1/2;t) term in the cylindrical coordinate system.Inspired by the study in fluid mechanics (Fukagata and Kasagi,2002; Zhuang et al.,2018),considering the fact that there is no such singularity in the Cartesian coordinate system,transformation from the cylindrical coordinate system to the Cartesian coordinate system could provide an interpolation solution to approximate the variable value of vprat the cylindrical axis.Thus,as shown in Fig.3,the value of vprat r=0 can be approximated by all of the nearby velocity vpθ via two intermediate variables vxand vyat r=0.Note that,a vector in the cylindrical coordinates (r,θ,z) can be expressed in the Cartesian coordinates as er=excosθ+eysinθ,one can obtain the procedure of handling the singularity in P-wavefield as follows.

The intermediate variables vxand vyare calculated from vpθ:

where N denotes the number ofin the circle near cylindrical axis and θj=jΔθ.

The vpris calculated from the intermediate variables vxand vy:

As for S-waves,similar observation can be made about the missing variables when calculating vsr,τsrθand τszrat the cylindrical axis using FDTD method.Different from the P-wavefield,S-wavefield is defined by curl operation as

Fig.3.Cross-section diagram of P-wavefield near cylindrical axis at z=k+ 1/ 2.

Thus,by Stokes’ theorem (Stewart,2012),we can calculate the singularity variable using the surrounding variables on the plane perpendicular to the singularity variable.For example,as shown in Fig.4,the vector vsrat r=0 can be calculated using the vectors τsrθand τszron the circle adjacent to vsr.Similarly,vector τrθat the cylindrical axis can be calculated from vsθat the circle of radius r=Δr/2 and vector τs2rat r=0 can be obtained based on the vsrand vsz.The detailed wavefield calculation formulae are summarized as follows:

To demonstrate the validity of the proposed singularity treatment,we propagate wavefield in a homogeneous velocity model.As shown in Fig.5,by loading the excitation source at the cylindrical axis,a stable wavefield can be numerically calculated.That is,both the P-and S-waves can propagate through the cylindrical axis without obvious diffraction,which lays the groundwork for accurate imaging in tunnels.

2.3.Procedure of the cylindrical ERTM in tunnels

Similarly with the normal ERTM procedure in tunnels,decoupled non-conversion ERTM using cylindrical coordinates also requires forward and backward wavefield extrapolation and proper imaging conditions(Wapenaar et al.,1987).To ensure a satisfactory migration result,separated P-and S-waves are usually used to mitigate crosstalk artifacts,and its good application effect in tunnel forward prospecting has been verified under the 2D Cartesian coordinate system (Ren et al.,2020).In this work,we extend the similar philosophy into 3D cylindrical coordinates.After addressing the key issues in wavefield extrapolations and imaging conditions,an effective seismic imaging procedure based on the decoupled non-conversion elastic wave equation suitable for geological forward-prospecting in tunnels is proposed,whose application effect will be further illustrated via numerical examples and field cases.The detailed CERTM imaging procedure shown in Fig.6 can be described in the following aspects.

Fig.4.S-wavefield singularity calculation strategy for (a) vsr,(b) τsrθ and (c) τszr.

Fig.5.Wavefield snapshot of vr component for (a) the original elastic wavefield,(b) the decoupled P-wavefield,and (c) the decoupled S-wavefield.

Firstly,in tunnel seismic detection,hammer sources are usually excited on tunnel sidewalls.Cylindrical coordinates could provide an accurate simulation of the hammer source in terms of the source excitation position and direction.Thus,when conducting forward wavefield extrapolation in CERTM,we prefer loading the source wavelet (e.g.Ricker wavelet) onto the variable vrto simulate the particle movement caused by hammering.As for the backward wavefield extrapolation,data preprocessing like τ-p filtering (e.g.Russell et al.,1990)is required to separate the practical observation data into P-and S-waves and the different components of the separated waves will be injected into the corresponding variables vr,vθand vz,respectively,as the initial condition.

As for the imaging condition,considering that the wavefield obtained via Helmholtz decomposition includes a scalar stress field τPand a vectorial stress field τS,a scalar cross-correlation imaging condition(Du et al.,2012)is adopted between PP-waves or SS-waves.Moreover,as shown in the roz planes in Fig.5,there is a wavefield polarity reversal phenomenon that can be observed at the cylindrical axis.This is because of the radius direction change in cylindrical coordinate system and it does not appear in the actual wavefield.Fortunately,the cross-correlation imaging condition can directly eliminate the influence of the wavefield polarity reversal during wavefield multiplication.This is another reason why we prefer the cross-correlation imaging condition over other imaging conditions like excitation time or maximum amplitude imaging conditions.Moreover,illumination compensation and Poynting vector are also adopted in this study to improve the quality of CERTM result.

Fig.6.Flow chart of the proposed 3D CERTM in cylindrical coordinates to migrate seismic forward-prospecting data in tunnels.

Last but not least,multiple sources and survey lines are usually arranged on tunnel sidewalls,and the seismic data received via all survey lines will be imaged to present the single-shot PP-and SSimages.The final PP-and SS-images can be obtained via superposition of the PP-and SS-images from all sources,respectively.Integrated analysis will be performed on both final PP-and SSimages to interpret the geology in front of the tunnel face.

3.Numerical examples

Two sets of numerical experiments are conducted in this study.The first one compares the wavefield features simulated in Cartesian and cylindrical coordinates.After demonstrating the superior modeling results using cylindrical coordinates,its application effect in CERTM images is further verified on the basis of several typical geological models in tunnels.

3.1.Wavefield comparison in different coordinates

Currently,most researches in tunnels are based on Cartesian coordinates,where the 3D tunnel body is approximate to a cylinder body by cubes.As mentioned before,such approximation fails to accurately simulate the seismic wavefield around the tunnel body,leading to errors in seismic imaging and inversion.To clearly elaborate the advantage of using cylindrical coordinates to study the seismic forward-prospecting in tunnels,a comparison experiment is carried out based on the tunnel models and the seismic survey lines shown in Fig.7.

Fig.7.Diagram illustrating the tunnel model and the survey line setups in (a) Cartesian coordinates and (b) cylindrical coordinates.

The tunnel model size and data observation setup are similar for the experiment in both coordinate systems.Taking the model in Cartesian coordinates as an example,the tunnel body is 30 m long,12 m wide and the total model size is 60 m × 60 m × 140 m with grid spacing of 0.5 m.To record the seismic wavefield data,two survey lines are set on the boundary of the tunnel face and along the z axis 20 m ahead of the tunnel face.A 300 Hz source is placed on the top of the tunnel body with 5 m away from the tunnel face.The tunnel body is filled with air while the surrounding rock is set to be homogeneous with constant P-wave velocity of 4000 m/s and S-wave velocity of 2300 m/s.

Fig.8.Snapshots of wavefield in (a,b) Cartesian xoy plane and (c,d) cylindrical ro plane.(a) and (c) show the wavefield at time step of 10 ms and (b) and (d) show the wavefield at 15 ms.

By comparing the full elastic wavefield snapshots shown in Fig.8,one can clearly see that there is diffraction interference generated at the cubic corners and edges of the tunnel sidewall in Cartesian coordinates.In the tunnel seismic wave field,the energy of diffracted wave is similar to that of effective reflected wave.Therefore,in the process of RTM imaging,the reverse-time propagated wavefield will cross-correlates with not only the effective signal,but also the diffraction wave,leading to serious artifacts and misleading the interpretation of the migration imaging results.While in cylindrical coordinates,seismic waves propagate on the tunnel sidewall with no encounter of corners,thereby,no obvious diffraction around the tunnel face.Similar observation can be made in the seismic records shown in Fig.9.Compared with the Cartesian coordinates,cleaner seismic records can be obtained in cylindrical coordinates with no obvious diffraction.In a word,the clearer wavefield in cylindrical coordinates will naturally lead to a clearer migration results under the cross-correlation imaging condition and the cylindrical coordinate can be considered as a better choice to simulate seismic wavefield in tunnel environment than Cartesian coordinates.

3.2.Cylindrical ERTM imaging in tunnels

Fig.9.Seismic records of wavefield in(a,b)Cartesian xoy plane and(c,d)cylindrical ro plane.(a)and(c)show the seismic data recorded by the survey line 1 and(b)and(d)show the seismic records by the survey line 2.

To further validate the application effect of cylindrical coordinates in tunnel seismic imaging and demonstrate the imaging effect of CERTM,numerical examples are performed on the basis of two typical geological models in tunnels.We use the velocity of surrounding rock as the migration velocity and free-surface boundary condition is adopted to simulate the wavefield around tunnel body.

Fig.10.The designed tunnel models:(a)Model with a vertical interface and(b)Wedge fault model.

3.2.1.Tunnel models and data acquisition setups

In tunnel excavation,faults and fractured zones are common geological structures that could potentially lead to severe geohazards.Thus,tunnel models including one or multiple interfaces of large dip angels are often studied in the field of geotechnical engineering.In this experiment,a vertical interface model and a wedge fault model are designed based on a cylinder of 35 m in radius,as shown in Fig.10.In both models,the grid spacing is set to be Δr=Δz=0.5 m and the time step interval is 0.2 ms.The tunnel body in both models is 40 m in length and 12 m in diameter.Both the vertical interface and the wedge fault are placed 40 m in front of the tunnel face.The detailed physical parameters of different media in both tunnel models are summarized in Table 1.

Table 1 P-and S-wave velocities of different media in tunnel model.

The data observation system in 3D tunnel environment is shown in Fig.11.Considering the typical geological structures used in this experiment,only two hammer sources are excited on both tunnel sidewalls,respectively,and four survey lines are placed on the left,right,top and bottom of tunnel sidewall.Each survey line contains 12 receivers with spacing of 2 m.In order to record as much seismic data as possible,we set the left and right survey lines 18 m away from the tunnel face while the top and bottom survey lines are 13 m away.Both seismic sources are set 12 m away from the tunnel face and they are simulated by a 100 Hz Ricker wavelet.

3.2.2.CERTM results

For the vertical interface model,we show the single-shot CERTM imaging results and the final stacked CERTM imaging results in Fig.12.The single-shot result is obtained using the source excited on the right tunnel sidewall without illumination compensation,while the final results are stacked over all sources on both tunnel sidewalls with illumination compensation.From both the PP-and SS-imaging results,one can clearly see the energy arc corresponding to the vertical interface.In fact,given a correct surrounding rock velocity as migration velocity,the energy arc should be tangent to the actual interface right in front of the tunnel face.Moreover,the stack of single-shot images and the application of illumination compensation method can result in an improved image with more balanced amplitude and less artifacts,especially for S-waves.Generally,PP-images could demonstrate a superior imaging effect on simple geological interfaces,while the SS-images tend to have higher resolution on local structures.Thus,geological interpretation is recommended to mainly rely on the PPimages,supplemented by the SS-images.

Fig.11.Seismic observation setup in 3D tunnel environment.Four survey lines and two sources are placed on tunnel sidewall.

Fig.12.The CERTM results of the vertical interface model: (a) Single-shot PP-image,(b) Single-shot SS-image,(c) Stacked PP-image,and (d) Stacked SS-image.The white dashed lines represent the vertical interface.

Currently,most RTM in tunnels is conducted by directly loading one component of the recorded data into acoustic wave equation.Although such approximate treatment in Cartesian coordinates can result in an image successfully denoting the existence of geological structures,its application in cylindrical coordinates may suffer polarity reversal at the cylindrical axis (see Fig.13 for instance).Moreover,due to the lack of elastic information,the energy in the acoustic RTM (ARTM) result seems to be more focused on the arc tail,resulting in a limited ability to detect the geology right in front of the tunnel face.Thus,the proposed CERTM method can be regarded as a preferred choice for seismic imaging in tunnel geological forward-prospecting.

Fig.13.Single-shot ARTM imaging result of the vertical interface model.The white dashed line represents the vertical interface.

As for relatively complex geological structures,Fig.14 shows the final imaging results of the wedge fault model.From the PP-and SSimages on the vertical cross-section shown in Fig.14a and b,one can see that the proposed CERTM is capable of identifying the two fault interfaces.Particularly for the inclined interface,the energy arcs are focused on the lower part of the model and its tangent line corresponds to the actual interface.Thus,one can deduce the dip angle of the actual interface via these energy arcs.Moreover,the SSimages have higher resolution than the PP-images,which will provide important supplementary information about small-scale geological structures,such as the lower part of the wedge fault.From the images on the horizontal cross-section shown in Fig.14c and d,one may notice the positioning error regarding to the second interface,due to the inaccurate migration velocity.However,this type of error is acceptable since tunnel engineers often pay more attention to the location and shape of the first interface,while the imaging information from the later interfaces is used as an aid in making tunneling plans.If necessary,more detections should be arranged after tunneling through the first geological structure and reliable information about the subsequent geological conditions can be obtained.

Fig.14.The CERTM imaging results of the wedge fault model:(a)Final PP-image on the vertical cross-section(0°-180°),(b)Final SS-image on the vertical cross-section(0°-180°),(c)Final PP-image on the horizontal cross-section(0°-180°),and(d)Final SS-image on the horizontal cross-section(0°-180°).The white dashed lines represent the interface of the wedge model.

Fig.15.Observation setup used in practical seismic forward-prospecting (Liu et al.,2018).

4.Case study

To further validate the proposed CERTM method in practical seismic forward-prospecting situation,we conducted a field test based on a tunneling project in the Qinling Mountains,Shaanxi Province,China.This tunnel is part of the water diversion project that transports water from the Hanjiang River to the Weihe River.The overall geology of the tunneling area belongs to Caledonian granite with developed underground water system.It was excavated using the drilling-and-blasting method,and before this seismic detection,tunneling activities had been stopped due to several water-inrush accidents.Thus,to ascertain the geological structures between the mileage of K68+932-K68 +832,seismic forward-prospecting was conducted at the mileage of K68 + 932.Actually,seismic imaging result based on travel time has been reported in our previous work(Liu et al.,2020),and here we further use the proposed CERTM to image the same data.

Fig.15 shows the seismic data observation system used in the field case.Limited by the practical tunneling situation,12 hammer sources in two columns are triggered 5 m or 7 m away from the tunnel face,respectively,and totally ten seismic receivers are placed on both tunnel sidewalls,20 m away from the sources.From the received data (see Fig.16 for example),one can see multiple reflected waves in the separated P-wave record and only one strong seismic event can be observed in the S-wave record.Considering the P-and S-wave propagation features in water-bearing media,it was preliminarily predicted that there is a water-rich fissure zone ahead of the tunnel face,hence the strong S-wave reflection can be recorded.

Since it is impractical to obtain an accurate migration velocity,the P-and S-wave velocities of the surrounding rock are adopted to migrate the detection data.As shown in Fig.17,three main reflection areas can be observed in the final PP-image.The first reflection area is roughly within 20 m in front of the tunnel face.Considering the strong reflection in S-wave record,we assume it is quite possible to be a water-bearing area.Moreover,from the PP-image,we can observe two more reflection areas about 40 m and 80 m away from the tunnel face,respectively.Combined with the geologic information from preliminary surface exploration,developed joint fissures with high risk of water inrush may exist in these reflection areas and further measures are required to ascertain the water-bearing condition.

Fig.16.Practical seismic data corresponding to the receivers 1-5 and the source S11:(a)Raw data recorded,(b)Separated P-wave recorded,and(c)Separated S-wave recorded.The red dashed rectangles denote the reflected P-waves in groups.

Fig.17.Imaging result of the field case:(a)PP-image in cylindrical coordinates and(b)Its cross-section on roz plane.The red dashed lines mark the predicted geological structure in front of the tunnel face.

Fig.18 shows the geological sketch and corresponding photos taken after tunnel excavation exposure(Liu et al.,2020).There were water inrush accidents occurring on the tunnel face and sidewall when tunneling between the mileage of K68 + 932-K68 + 913,which can be regarded as the first predicted reflection area.When tunneling at K68 + 913,three boreholes were drilled and waterbearing structures were encountered near the mileage of K68 + 894.Finally,at the K68 + 851,tunneling encountered another fractured zone with water,which corresponds to the third reflection area shown in the PP detection image.In a word,all three major reflection areas predicted by the CERTM imaging result can be verified by practical exposure,which illustrates the effectiveness of the proposed imaging method in practical tunnel projects.

5.Discussion

In this work,to accurately image the adverse geology in front of tunnels,cylindrical coordinate system is applied in elastic wavefield simulation based on the decoupled non-conversion wave equation.Compared with the common Cartesian coordinates,cylindrical coordinates could provide more accurate modeling of the actual tunnel shape,leading to more realistic simulation of observation data and more reliable imaging results.Although the proposed method can be easily extended in many other cylindrical fields like seismic logging,its extensive application is limited by several numerical issues pending to be resolved.

Firstly,there is a mathematical singularity problem when calculating the wavefield at the cylindrical axis.In consideration of the separable nature of P-and S-wavefields in the decoupled nonconversion elastic wave equation,the singularity issue is solved by separating the original elastic wavefield into a scalar P-wavefield and a vectorial S-wavefield.Moreover,given the curl relation between the velocity and stress variables for S-waves,according to the Stokes’theorem,accurate value of the variables vsr,τsrθand τszrat r=0 can be obtained by properly choosing the adjacent variables into calculation.In fact,similar computation strategy has been reported in electromagnetic field (Chi et al.,2010) and the physical laws behind the singularity computation strategy are Ampere’s law and Faraday’s law,whose mathematical derivation is also based on Stokes’ theorem.However,we have not found a physically meaningful way to stably calculate the P-wave variable at the cylindrical singularity and the exact way to calculate the full elastic wavefield at the cylindrical axis is still under discussion.

Fig.18.Geological sketch and photos of the tunnel excavation exposure (Liu et al.,2020).

As for the computation cost of using cylindrical coordinates,current grid refinement strategy will lead to a large number of cells in the θ direction.Thus,the computation cost of a large 3D model would become more expensive in cylindrical coordinates than that in Cartesian coordinates.For example,on the workstation of Gold 6150 CPU with 512 GB memory,it takes about 24 h to obtain the 3D RTM imaging results.Another down side of cylindrical coordinates arises from the modeling of geological body.Although it is convenient to simulate a geological structure symmetric to the cylindrical axis,geological bodies placed off the axis may suffer constraints on shapes from the cylindrical cells.Thus,it is difficult to design a tunnel model with complex geological bodies.Recent development in mesh-free methods,such as the element-free Galerkin method(e.g.Katou et al.,2009) and the generalized FDTD method (e.g.Shragge,2014),could provide a new perspective in simulating irregular geometry in tunnel models.

6.Conclusions

Seismic imaging using decoupled non-conversion elastic wave has demonstrated superior application effect in tunnels.To better model the tunnel environment and simulate the observation data in seismic forward-prospecting,cylindrical coordinate system is adopted to numerically simulate the seismic wavefield propagation.To address the singularity issue inherited in the wavefield simulation at cylindrical axis,a stable singularity treatment with clear physical meanings is proposed based on P-and S-wavefield separation.With accurate elastic wavefield calculated,procedure of CERTM in tunnels is summarized with application effect further validated in numerical examples and field cases.Both PP-and SSimages can present the geological structure in front of the tunnel face without crosstalk artifacts and high-resolution information can be inferred from the SS-images.After revealing the features of the PP and SS migration results in tunnel environment,the field case study indicates the potential of applying the proposed method in practical seismic forward-prospecting.In the end,insights and potential improvement of seismic wave simulation in cylindrical coordinates are discussed.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The research work described herein was funded by the National Natural Science Foundation of China (Grant Nos.52021005 and 51739007) and the Key Research and Development Plan of Shandong Province(Grant No.2020ZLYS01).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2022.02.012.