Forward,prediction,for,tunnel,geology,and,classification,of,surrounding,rock,based,on,seismic,wave,velocity,layered,tomography

Bin Liu,Jiansen Wang,Senlin Yang,Xinji Xu,Yuxiao Ren

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

b School of Civil Engineering,Shandong University,17923 Jingshi Road,Jinan,250061,China

c Data Science Institution,Shandong University,27 Shandanan Road,Jinan,250100,China

d School of Qilu Transportation,Shandong University,17923 Jingshi Road,Jinan,250061,China

Keywords: Tunnel geological forward-prospecting Seismic wave velocity Layered inversion Surrounding rock classification Artificial neural network (ANN)

ABSTRACT Excavation under complex geological conditions requires effective and accurate geological forwardprospecting to detect the unfavorable geological structure and estimate the classification of surrounding rock in front of the tunnel face.In this work,a forward-prediction method for tunnel geology and classification of surrounding rock is developed based on seismic wave velocity layered tomography.In particular,for the problem of strong multi-solution of wave velocity inversion caused by few ray paths in the narrow space of the tunnel,a layered inversion based on regularization is proposed.By reducing the inversion area of each iteration step and applying straight-line interface assumption,the convergence and accuracy of wave velocity inversion are effectively improved.Furthermore,a surrounding rock classification network based on autoencoder is constructed.The mapping relationship between wave velocity and classification of surrounding rock is established with density,Poisson’s ratio and elastic modulus as links.Two numerical examples with geological conditions similar to that in the field tunnel and a field case study in an urban subway tunnel verify the potential of the proposed method for practical application.

The complex geological environment has posed great challenges to tunnel construction.Unfavorable geological structures,such as fault fracture zones,weak strata,and karst caves,are frequently encountered,which may cause the instability of surrounding rock,water and mud inrush,collapse,and other disasters (Liu et al.,2022).For the safety of tunnel construction,the forwardprospecting of tunnel geology is a necessary process,which can identify the unfavorable geological structure in front of the tunnel face more accurately (Li et al.,2017).The seismic method (Lüth et al.,2009) has the advantages of long detection distance,interface depiction of geological structure,inversion of rock mass parameters,i.e.wave velocity and rock density(Alimoradi et al.,2008;Li et al.,2019;Ren et al.,2022),and can further learn the classification of surrounding rock.Therefore,this method is widely used in the forward-prospecting of tunnel geology.

Seismic wave velocity is one of the most critical indicators in the forward-prospecting of tunnel geology.It can directly reflect the general situation of the rock mass in front of the tunnel face,and provide the basis for migration imaging and classification of surrounding rock.However,due to the narrow tunnel space and strong vibration from tunnel construction,it is challenging to obtain the wave velocity information of the surrounding rock in front of the tunnel face through the seismic forward-prospecting.At present,the wave velocity of the tunnel surrounding rock is acquired from the seismic exploration on ground surface.Two feasible methods have been developed: migration velocity analysis based on raytracing theory (Gong et al.,2010) and full waveform inversion(FWI) based on waveform fitting (Bharadwaj et al.,2017;Wang et al.,2019).The migration velocity analysis takes the maximum migration energy as the criterion for the accurate estimation of wave velocity,which has the advantage of high computational efficiency.However,it lacks the procedure of data fitting,which may result in low estimation accuracy of wave velocity,especially for complex geological structures.The FWI based on wave equation theory considers more waveform information,such as seismic amplitude and phase in seismic data,in addition to the travel time information (Yu et al.,2021),which has the advantage of high accuracy (Riedel et al.,2022).However,it has problems of high computational cost and has difficulty in waveform fitting due to the noise in seismic data.Therefore,a satisfactory FWI mainly depends on high-quality data,which are often difficult to be obtained in the field tunnel.

Tomography (Bording et al.,1987) combines the advantages of the aforementioned two methods.Its good performance has been demonstrated by field data (Zhang and Thurber,2005),which is expected to be applied to seismic velocity estimation in tunnels.However,the application of tomography to tunnels faces difficulties,such as small offset distance,small effective ray coverage area,and few reflection points.This means that the observed data may contain little effective information of velocity,making the problem of multi-solutions more serious when the method is applied to large-inversion areas.Therefore,it is necessary to improve the tomography according to the field tunnel environment.

In addition,the wave velocity is related to the classification of surrounding rock (Kayabasi et al.,2003;Zhang,2017).The classification of surrounding rock is a main method for evaluating the quality of surrounding rock (Barton et al.,1974).In the process of tunnel construction,the classification of surrounding rock in front of the tunnel face can be continuously corrected based on the wave velocity results obtained by tunnel seismic forward-prospecting method,so as to facilitate tunnel excavation.Shi et al.(2014) and Bu et al.(2018)proposed an optimized classification of surrounding rock based on tunnel seismic prediction data.Various prediction algorithms have been proposed for the classification of surrounding rock at the construction stage.Zhou et al.(2015) used the gray evaluation model to assess the risk of surrounding rocks in front of the tunnel face,whose short-term prediction accuracy is high,and the long-term prediction tends to be unsatisfactory.Qiu et al.(2010) adopted the support vector machine algorithm to classify the surrounding rock of the tunnel,whose performance depends on the parameters of the learning machine,and the selection of parameters depends on personal experience.Artificial neural network(ANN)(Zhang and Phoon,2022;Zhu et al.,2022)has achieved great promotion in risk assessment in tunnel due to its powerful nonlinear fitting and classification capability (Hasegawa et al.,2019).Most of the aforementioned algorithms are the direct application of the existing networks,and it is difficult to achieve high-accuracy prediction of surrounding rock classification only by inputting the wave velocity information.Therefore,a proper network design should be carried out to realize a multiple perspectives comprehensive prediction of surrounding rock classification and improve the reliability for field tunnel.

To address the aforementioned issues (i.e.wave velocity estimation and classification prediction of surrounding rock),we presented a forward-prediction method for tunnel geology and classification of surrounding rock based on seismic wave velocity layered tomography.Firstly,in view of the problem that the wave velocity inversion has a strong multi-solution in the narrow tunnel space,a layered inversion for seismic reflection tomography and linear interface regularization were proposed,leading to a reduced number of variables for each inversion step,and subsequently an improved wave velocity inversion.Furthermore,an autoencoderbased surrounding rock classification network was developed,and the mapping relationship between wave velocity and classification of surrounding rock was established with density (ρ),Poisson’s ratio (μ),and elastic modulus (E) as links.Multi-parameter self-learning was carried out to further improve the accuracy of surrounding rock classification.Moreover,the feasibility of tunnel seismic layered tomography is verified by two numerical examples with geological conditions similar to that in the field tunnel.Finally,we conducted a detailed field case study to demonstrate the reliability of our approach in practical applications.

2.1.Layered tomography method for estimating seismic velocity in tunnel

In order to address the multi-solution of inversion and obtain accurate wave velocity of surrounding rock,a layered tomography method for estimation of tunnel seismic velocity was proposed,under the assumption of straight geological interface in front of the tunnel face.This is because the change of the geology in front of the tunnel along the axial direction has a great impact on the construction safety,and the engineers usually ignore the fluctuation of the geological interface and the radial change of the geology.Therefore,the results obtained by the assumption of a straight geological interface can satisfy the geological prediction in practical engineering.

2.1.1.Tunnel seismic reflection tomography

Ray-tracing is the basis of seismic tomography.According to the seismic ray-tracing theory,the relationship between the travel time of seismic wave and medium velocity is expressed in the following equation (Luxbacher et al.,2008):

whereTis the travel time of seismic wave,v(x,y) is the medium velocity,sis the slowness (inverse of velocityv),Sis the source position,Ris the receiver position,andlrepresents the sum of the distance from the source position to the reflection position and the distance from the reflection position to the receiver position in reflection tomography.Due to the advantages of high calculation speed and high accuracy(Moser,1991),the shortest path method is adopted to realize ray-tracing in this paper.

For the inversion problem of seismic reflection tomography,not only the slowness but also the reflection position(depth)needs to be updated.Thus,with a given initial velocity model,the relationship between the travel time residual δTand the slowness perturbations δscan be expressed according to Eq.(1)as(Korenaga et al.,2000):

wherezis the reflection position,and(∂T/∂z)δzmeans the change in travel time caused by the change of the reflection position.Discretization of Eq.(2) gives

where δd is the travel time residual vector;G=[GsGd] is the Fréchet derivative matrix (see Appendix A for details);δm=[δmsδmd]Tis the unknown model perturbation matrix;and subscripts s and d for the Fréchet submatrix and the model submatrix denote their slowness and depth components,respectively.

2.1.2.Layered tomography inversion

In the narrow tunnel,the numbers of sources and receivers are limited,which results in a small number of rays in reflection tomography.Moreover,the small offset may lead to a group of ray paths with high similarity,which may increase the ill-conditioned degree of the inversion equation.In order to address this problem,a layered inversion is proposed in this work.The number of reflective layers is firstly determined based on the picked-up reflections,and then the inversion objective function is established as follows:

whereKis the number of reflective layers;the superscriptkdenotes the inversion of thekth layer;is the observed travel time matrix of thekth layer,which is usually obtained by automatic pickup of travel time;is the theoretical synthetic travel time of thekth layer through ray-tracing forward modeling;Gkis the Fréchet derivative matrix of thekth layer;and mkis the model matrix of thekth layer.

The slowness and interface positions of each layer can be inverted via the following equation:

This inversion process is conducted layer by layer,starting from the area closest to the tunnel face and advancing in the direction of the construction excavation till the furthest layer of the inversion is reached.In the layered inversion process,each layer is inverted only using the ray path and travel time residuals of that layer.The inversion between layers does not interfere with each other.This is because after the updated layer is fitted,unknown ray paths may cause the deteriorated inversion results of the updated layer if the inversion region contains the updated layer.

The initial model for the tomography inversion consists of a slowness model and a depth model.The slowness model is built by picking up the velocity of the first arrival wave from seismic data.For the depth model,the observed seismic data need to be first processed by migration velocity analysis and diffraction-stack migration (Bellino et al.,2013;Zhan et al.,2014).Then,the position of reflections is picked up on the diffraction-stack migration as the initial depth model.

For further improvement of the inversion convergence,the linear system needs to be regularized to stabilize the inversion of tunnel seismic reflection layered tomography for tunnel.After regularization (details in Appendix B),inversion of Eq.(5) can be written as

whereIand subscriptirepresent the total number of iterations and theith inversion,respectively;Lsis the smoothing matrix for slowness perturbations,which conforms to the law of geological slow change;Ldis a post-process operator for the interface depth,which fits the interface coordinates obtained from theith inversion as a straight line;λsand λdare the regularization coefficients,which determine the relative importance of the regularization with respect to the resolution and also play the role of normalization;ω is the depth kernel weighting parameter,which adjusts the relative weighting of depth sensitivity in the Fréchet matrix;andis the depth model after regularization.

The inversion is a sparse system that can be efficiently solved by the sparse matrix solver least-squares QR-factorization (LSQR)(Paige and Saunders,1982) to obtain the model updated amount δm.

Although solving a single sparse system is computationally inexpensive,obtaining the optimal combination of regularization coefficients and depth kernel weighting parameter requires multiple trials per iteration.Therefore,the recommended values of regularization coefficients and depth kernel weighting parameters that meet the conventional geological conditions are given in Appendix B.For complex geological structures,we recommend first testing the optimal combination of regularization coefficients and depth kernel weight parameter through a preliminary one-step inversion,and then fixing them in subsequent iterations.At the later stages of the iteration,the model perturbations tend to become small and the model may become rough.To address this,we add post-inversion smoothing to the slowness model in the iteration:

2.1.3.Tunnel layered tomography process

For the raw data obtained from seismic forward-prospecting for tunnel,they contain numerous noises and the P-and S-waves are coupled,which makes it difficult to pick up the travel time of reflections.Thus,the raw seismic data need to be pretreated via data denoising and wavefield separation (F-K or τ-pfiltering)(Greenhalgh et al.,1990;Duncan and Beresford,1994).The first arrival wave is further picked up from the raw data for building the initial model.The reflections are obtained from the pretreated data for calculating the travel time residuals during the inversion process.Then,the proposed layered tomography is adopted to obtain the wave velocity information of the surrounding rock.Fig.1 demonstrates its specific process,which can be summarized as follows:

Fig.1.Flow chart of the inversion of tunnel seismic reflection layered tomography for tunnel.

(1) Pretreat the raw seismic data;

(2) Pick up first arrival wave and reflections;

(3) Build the initial model;

(4) Perform iteration of each layer;

(5) If travel time residuals or number of iterations meet the requirements,output the results of each layer;

(6) If all layers are traversed,output the final inversion results.

Note that the number of picked-up reflections,the wave velocity of the first arrival,and the inversion results of the previous layer can all be considered as constraints in inversion to ensure the stability of the iteration.

2.2.Estimated classification of surrounding rock

Although the P-and S-wave velocities,which can be obtained by tunnel seismic layered tomography in Section 2.1,are considered to have a definite relationship with the classification of surrounding rock(Kayabasi et al.,2003),it is difficult to satisfy the classification of surrounding rock under practical conditions with these data alone.We introduced correlated rock property parameters (i.e.density,Poisson’s ratio,and elastic modulus)to improve estimation accuracy.Through the training process of numerous data,ANN can extract features from the data and summarize the nonlinear mapping between input and output(Guo et al.,2022).After the training of ANN,direct data input can quickly predict the results in 1 s.Therefore,the property parameters of rock mass obtained from field boreholes are used as labels.The network training is conducted in a supervised manner.The trained network can be used to achieve a rapid classification of surrounding rock at the construction stage of the tunnel.

In this work,an autoencoder-based ANN is designed to predict the classification of surrounding rock.In order to incorporate the information of density,Poisson’s ratio,and elastic modulus into the classification prediction process of ANN,two encoders and training manner with self-learning are used to extract the features of highdimensional data,which would be subsequently beneficial for an accurate classification of surrounding rock.As shown in Fig.2,the designed ANN includes two encoders and two decoders.One of the decoders is for self-learning and the other will output the probability of different surrounding rock classifications.

Specifically,the first encoder (Enc1) is used to predict the density of surrounding rock,which is extremely difficult to be obtained by direct inversion or inaccurate obtained by empirical petrophysical formula (Gardner et al.,1974;Gaviglio,1989).Therefore,the P-(VP) and S-wave (VS) velocities are regarded as input to estimate the density(ρ)of the rock mass to offer more information for further classification of surrounding rock.For the other encoder(Enc2),five relevant parameters are regarded as input to extract latent features(fea),including ρ obtained by Enc1,Poisson’s ratio μ;elastic modulusE,VPandVS.Among them,onceVP,VSand ρ are determined,the Poisson’s ratio μ and elastic modulusEcan be obtained by the following equation:

Fig.2.Designed ANN structure.

In this way,the number of rock property parameters is increased from two(VPandVS)to five with the assistance of Enc1,and Enc2is used to understand the inherent characteristics among the parameters.

There are also two decoders.One (Dec1) is used to return the input data for extracting the latent featurefea,while the other(Dec2)is used to predict the classification of surrounding rock usingfea.By employing the decoder Dec1to return the input data directly,different from the standard ANN(only Enc2and Dec2),it is anticipated that the network would extract the indispensable features of the input to obtain the probability of each classification.With a limited number of training sets,it may assist the network to enhance the learning impact to obtain an accurate prediction of rock classification.

The whole workflow can be described as follows:

where θi(i=1;2;3;4) represents the updateable network parameters;ρ′′is the density predicted by Enc1;,,ρ′,μ′andE′are the predicted results to return the input data;andPis the expected probability of each classification of surrounding rock.

Moreover,to achieve the nonlinearity of ANN,the rectified linear unit(ReLU)is used as the activation function.Table 1 shows neural network parameters in detail.A loss function is used to guide the learning of neural networks and achieve the desired outcomes.For our proposed ANN,a hybrid loss function is used to provide constraints for the ANN,including cross-entropy and mean square error(MSE),as shown in Eq.(12):

Table 1 The number of nodes in the ANN.

Table 2 MSE values of tomographic results in the fault model.

whereRCis a vector with five unit of length,whose value equals 1 at the true category and 0 elsewhere;α is the weighting factor,usually determined empirically.Loss1 (cross entropy),the most commonly used loss function in classification,determines the category by predicting the probability distribution,which also provides us with the probability of surrounding rock classification.For the autoencoder and self-learning mentioned above,the Loss2 function is able to be interpreted as a constraint imposed on ANN training by recovering the input data itself,so as to obtain the critical featuresfeaand make the predictionPaccurate.In addition,it should be noted that Enc1should be pre-trained with the sameL2-norm loss function,whose parameters would be employed directly in the following network:

where Lossρ is the loss function for the density prediction.

The relative geological borehole information near the field tunnel is collected to build the training datasets.The size of the training dataset would affect the final prediction of the surrounding rock classification.We recommend that the borehole data near the field tunnel should exceed 100 groups to ensure the stability of the training network.If numerical simulation is used to construct training dataset,it is necessary to consider transfer learning and generalization problems,and more subjective factors from manual labels are carried.

To verify the feasibility of the proposed tunnel seismic layered tomography,two types of numerical examples were designed:fault model and multi-layer interface model.For each model,we compared the conventional tomography and the layered tomography.In order to better demonstrate the superiority of the layered inversion strategy tomography,the inversion convergence curve and MSE of the tomographic results are also presented.For a quick estimation of classification of surrounding rock in tunnel,supervised training of the network is performed using field borehole information as labels,which is difficult to simulate by numerical methods.Therefore,the practical application of surrounding rock classification network will be verified with field cases in Section 4.

3.1.Model parameter

In tunnel construction,faults,fractured zones,weak intercalated layers,and lithologic interfaces at different angles have significant impacts on excavations.Thus,we designed a fault model (Fig.3a)and a multi-layer interface model (Fig.3b) to investigate the applicability of the layered inversion strategy tomography.Both models contain 45 ×140 grid cells with a uniform grid spacing of 1 m in both directions.The length and width of the tunnel are 40 m and 6 m,respectively,and the wave velocity of the tunnel is 340 m/s.For the fault model,the wave velocity of the tunnel surrounding rock (area A in Fig.3a) is 3500 m/s,the wave velocity within the fault zone (area B in Fig.3a) is 2000 m/s,and the dip angle of the fault is about 75°.For the multi-layer interface model,the wave velocity of the tunnel surrounding rock (the first layer,area A in Fig.3b)is 3500 m/s,the wave velocity of the second layer(area B in Fig.3b)is 3000 m/s,and the wave velocity of the third layer(area C in Fig.3b) is 2500 m/s.

In addition,two types of observation layout in the tunnel are shown in Fig.4.In the drilling-blasting tunnel,the tunnel face can be used to fix the seismic source.Thus,sources and receivers are generally placed at the tunnel face and side wall,respectively (i.e.observation 1).Three sources are arranged with an interval of 2 m on the tunnel face,six receivers are arranged with an interval of 2 m on each side wall,and the receiver array is 13 m away from the tunnel face.However,in the tunnel construction using tunnel boring machine (TBM),the tunnel face is covered by the cutter head.Hence,sources and receivers are usually mounted at the tunnel side wall (i.e.observation 2).Two sources and six receivers are arranged at intervals of 1 m and 2 m,respectively,on each side wall,and the distance between the receiver array and the source array is 2 m.Therefore,we carried out tomography numerically for the two construction conditions.

The termination condition of the inversion iteration is that the travel time residual is less than 10-5or the number of interactions is no more than 10.

3.2.Analysis of tomographic results

Fig.5 displays the tomographic results in the fault model under different observation layouts and different inversion strategies.For conventional tomography,the interface information is vague,the inversion results of wave velocity inside the fault are not uniform,and the wave velocity in some areas is inaccurate.For layered tomography,the results are satisfactory under both observations 1 and 2,and the wave velocity and structure are generally consistent with that in the original model.

The convergence curves in the fault model obtained by theL2-norm objective function (Eq.(4)) are displayed in Fig.6,which includes the inversion of conventional tomography and that of tomography with layered strategy at each layer.It can be seen that the convergence curve of conventional tomography fluctuates greatly,and the final travel time residual is around 0.02.The convergence curves of the first and second layers under observations 1 and 2 are relatively stable,and their travel time residuals can all converge to less than 0.001.Note that since there is no reflection interface at the third layer,the wave velocity of the third layer cannot be inverted in reflection layered tomography.

Fig.5.Tomographic results in the fault model:Conventional tomographies under(a)observation 1 and(b)observation 2,and layered tomographies under(c)observation 1 and(d)observation 2.A and B represent the rocks with different velocities.The unit of color bar is m/s.

Fig.6.Convergence curves of tomography in fault model under (a) observation 1 and (b) observation 2.

In order to further evaluate the tomographic results in the fault model,the MSE of different tomographic results is calculated by the following formula (Liu et al.,2021):

Fig.7.Tomographic results in multi-layer interface model: conventional tomographies under (a) observation 1 and (b) observation 2;(c) layered tomographies under (c) observation 1 and (d) observation 2.A,B and C represent the rocks with different wave velocities.The unit of color bar is m/s.

where mcaland mtruare the calculated tomographic results and the true values,respectively;andnx∙nzrepresents the size of the model.The MSE values of different tomographic results in Fig.5 are shown in Table 2.The MSE of layered tomography is one order of magnitude lower than that of conventional tomography,and the MSE value of layered tomography under observation 1 is the smallest(3.8029).From the perspective of MSE,it also can be seen that the tomographic results under observation 1 are slightly better than that under observation 2.

For the multi-layer interface model,its tomographic results under different observations and different inversion strategies are shown in Fig.7.A low-velocity region can be seen in conventional tomographic results under the observations 1 and 2.Nevertheless,the interfaces in both the second and the third layers are difficult to identify,and the wave velocity value of the third layer is inaccurate.Moreover,the interface angle inversion of the third layer in the observation 2 conventional tomography is wrong.In the tomographic results with layered inversion strategy under observations 1 and 2,the wave velocities and structures of the three layers are more accurate.Additionally,the inversion result of the wave velocity at the third layer under observation 1 seems to be closer to the original model than that under observation 2.

Convergence curves and MSE values of different tomographic results are also calculated,as shown in Fig.8 and Table 3,respectively.Both conventional tomography and layered tomography can converge stably.The travel time residuals of conventional tomography can converge to less than 0.01.For layered tomography,thanks to the known direct wave information of the surrounding rock,tomography at the first layer converges to 10-5in only four iterations.The travel time residuals of tomography at the second and third layers converge to 0.001.By observing the MSE of the tomographic results,although there is no order of magnitude difference between the two methods,the MSE of the layered inversion strategy tomography is smaller.The minimum MSE value is still from the layered tomographic results under observation 1(1.1978).

Table 3 MSE values of tomographic results in the multi-layer interface model.

To further verify the proposed forward-prediction method for tunnel geology and classification of surrounding rock based on seismic wave velocity layered tomography,a case study was conducted in a field tunnel in Qingdao,Shandong Province,China.This tunnel is located between Chuangzhigu Station and Shishanlu Station of Qingdao Metro Line 6,which is constructed by doubleshield TBM.Influenced by the tectonic fracture,the surrounding rock is mostly fragmented.Joints and fissures have been developed in this area.In addition,numerous clay minerals are often filled along the joint surfaces,resulting in low bond strength between rock masses.The geological profile of the ZCK22+300-ZCK22+700 m segment in the preliminary survey is shown in Fig.9.The surrounding rock of the tunnel in this segment is mainly slightly weathered,moderately weathered and strongly weathered granite.In particular,in the range of ZCKZ22+550 m to ZCK22+450 m,passing through the Caijiazhuang Reservoir,the surrounding rock is fragmented and uneven in hardness and softness,which may cause abnormal damage to the TBM cutters,and even cause the TBM to jam in severe cases.When the tunnel was excavated near ZCK22+530 m,the seismic forward-prospecting for the tunnel was carried out in time to ensure the safe construction of TBM.

The seismic forward-prospecting was conducted at ZCK22+530 m,and the detection distance was 100 m from ZCK22+430 m.Fig.10 illustrates the observation layout in the field tunnel constructed by the double-shield TBM.Two aerodynamic sources were mounted between the TBM telescopic shield and support shield.Twelve receivers were distributed on both sides of the TBM in two survey lines.Restricted by the tunnel segment,the receiver can only be placed at the grouting hole of the segment,and the receivers all use plug geophones.The source was 5 m away from the tunnel face,the interval between receivers was 3 m,and the distance between the source and the receiver was 20 m.The sampling frequency was 5 kHz and the sampling time was 0.512 s.

The observed seismic data were processed by data denoising,separation of P-and S-waves,etc.,and then the tomography was performed on the P-and S-wave seismic data,respectively.The tomographic results are shown in Fig.11,which shows an area of 30 m×100 m in front of the tunnel face.Around 20-70 m in front of the tunnel face,there is a low-velocity area in both the P-and Swave tomographies,which may reflect the moderately weathered granite in the preliminary survey data.The previous borehole data show that the wave velocity of the slightly weathered granite in this area is 4198 m/s,and the wave velocity of the moderately weathered granite is 3536 m/s,which is roughly consistent with the layered tomographic results of wave velocity.The borehole acoustic test was conducted at ZCK22+465 m (location A) and ZCK22+510 m (location B) before excavation in the studied area,and the boreholes went across the tunnel.The comparison between the borehole acoustic test results and the tomographic results is shown in Fig.12.The blue curve is the borehole velocity.The red curve is the tomographic velocity at the corresponding position,which roughly matches the variation trend of the borehole velocity,especially in the excavation area of the tunnel.Moreover,the velocity curve of the borehole fluctuates greatly,while the tomographic velocity curve is relatively smooth.The average wave velocity of the borehole is generally consistent with that of the tomographic results,which is sufficient to meet the construction requirements.

Fig.8.Convergence curves of tomography in multi-layer interface model under (a) observation 1 and (b) observation 2.

In this case study,400 borehole datasets are selected near the research location and the corresponding classifications of surrounding rock are tagged as labels,to build datasets for networks.They are divided into the training set and validation set in the ratio of 7:1.Under the workstation with an Intel Xeon(R)gold 6150 CPU and two Tesla P100 GPUs,network training of 200 epochs takes about 10 min (including density network and surrounding rock classification network).After training,it usually takes less than 1 s to estimate the classification of surrounding rock.As shown in Fig.13,all three loss functions converge to a low level,indicating that the trained network has the ability to accurately predict the classification of surrounding rocks.

Fig.9.Geological profile and geographic location of the tunnel.

Fig.10.Observation layout in the field tunnel constructed by the double-shield TBM.

Fig.11.Tomographic results in field case:(a)P-wave and(b)S-wave.The unit of color bar is m/s.

By inputting the P-and S-wave velocities in front of the tunnel face obtained by the layered tomography to the trained network,the density,Poisson’s ratio and elastic modulus can be obtained,as shown in Fig.14.The property parameters of rock mass vary greatly in the range of ZCK22+455 m to ZCK22+510 m,indicating that the type of surrounding rock may change.The predicted results obtained by further classification of the surrounding rock are shown in Fig.15.It shows the predicted probability of surrounding rocks of classes I-V in the detection area.As shown in Fig.16b,the surrounding rock classification corresponding to the maximum probability at each mileage is selected as the predicted classification of surrounding rock for that mileage.Actually,at the design stage,the classifications of surrounding rock(Fig.16a)contain Class III (ZCK22+400-ZCK22+450 m) and Class V (ZCK22+450-ZCK22+560 m).After the correction by the rapid estimation of the surrounding rock,in the range of ZCK22+510 m to ZCK22+530 m,the Class V surrounding rock at the design stage is corrected to the Class III surrounding rock.Meanwhile,there is a transition of the Class IV surrounding rock in the process of changing from the Class V to Class III surrounding rock.Accordingly,it is speculated that the surrounding rock may become more fragmented in the range of ZCK22 +460 m to ZCK22+510 m.

In the construction process of double-shield TBM,the segments need to be assembled after TBM excavation,and the classification of the surrounding rock is unavailable due to the inability of direct observation on the surrounding rock.Therefore,we show the discharged rock slag of this section (Fig.17) during TBM excavation.From the slag,it can be seen that the surrounding rock changes from slightly weathered granite to moderately weathered granite to strongly weathered granite and then back to slightly weathered granite,which is roughly consistent with the variation of wave velocity and the estimated classification of surrounding rock from our seismic forward-prediction.In addition,in the vicinity of ZCK22+500 m,the speed of the TBM cutter head is increased from 4 r/min to 6 r/min,the torque is reduced from 1000 kN m to 800 kN m,the total propulsion force is reduced from 6000 kN to 4500 kN,and the penetration is increased from 18 mm/r to 22 mm/r.This also verifies the prediction of the deterioration of the surrounding rock in this section from another aspect.

This case study illustrates the reliability of the proposed method for the forward prediction of the tunnel geology and the estimated classification of surrounding rock.It also demonstrates the potential of this method in keeping the tunnel construction safe and optimizing construction measures.

A seismic layered tomography method was adopted to obtain the wave velocity of the surrounding rocks in front of the tunnel face.Furthermore,the estimated classification of surrounding rock was achieved by ANN.This method can dynamically correct the initial measurement results at the design stage in time at the construction stage.

In seismic reflection tomography,the accuracy of picking up travel time is crucial,which greatly affects the final inversion results.The existence of data noise will increase the difficulty in picking up travel time.After the data are denoised,if the data are relatively clean,the automatic pick-up method (i.e.STA/LTA algorithm (Li et al.,2016),kurtosis algorithm (Galiana-Merino et al.,2008) and ANN (Yuan et al.,2020)) can be used;if the reflection is difficult to identify,only manual pick-up can be used.Errors are unavoidable in the picking up travel time process of seismic data in field tunnel.Therefore,in the future,Monte Carlo and other methods will be used to carry out uncertainty analysis on pick-up errors of travel time.

The proposed tomography method based on straight-line assumption can depict the structural interface with good resolution in the lateral direction (i.e.tunnel axial direction) and poor resolution in the longitudinal direction.In practical engineering,large-scale geological structures (i.e.faults and fracture zones) are usually unbeneficial to tunnel construction.Geological changes along the axial direction in front of the tunnel are the main concern of engineers.In this sense,the assumption of a straight geological interface is reasonable,and the wave velocity results obtained based on this can meet the needs of geological forward-prospecting in field engineering.In addition,due to the small number of sources and receivers in tunnel seismic forward-prediction,there are few seismic rays covering the area ahead of the tunnel,resulting in a small amount of effective information available for tomography.With the straight-line interface assumption,the tomography inversion problem will benefit from a reduced solution space,alleviating the inversion artifacts associated with multi-solution problems,and thus giving better inversion results of wave velocity.While inverting a curved geological interface,it would also fit its boundary as a straight line,which is a forced simplification and approximation.These errors are acceptable in the field tunnel,and the inversion velocity results can also provide guidance for tunnel construction.However,it is still difficult to invert for karst caves and discontinuous interfaces,and it may need to be combined with other geophysical methods for joint inversion.Possible solutions to this issue include proposing a better assumption of interface morphology,improving the wave propagation mode,and designing a proper inversion objective function.

For the layered inversion strategy,while it can successfully improve inversion accuracy,the computational cost rises with the number of layers.We will consider parallel acceleration in the future to ensure the efficiency of application in engineering.Moreover,since the last layer in the inversion region has no reflective interface,the ray paths within its layer cannot be obtained.The layered tomography method is difficult to accurately invert the wave velocity for the last layer.Therefore,this may lead to low reliability of inversion results at long distances in practical applications.Combining with FWI may ameliorate this deficiency.In addition,existing borehole data and other preliminary survey data might be utilized to improve the inversion accuracy and convergence of the field data.

Fig.12.P-wave velocity comparison between tomographic results and borehole results at (a) location A and(b) location B.

Fig.13.Convergence curves for (a) Loss1,(b) Loss2 and (c) Lossρ functions.

Fig.14.Property curves of surrounding rock along the axial direction of the tunnel.

Fig.15.Prediction results of classification of surrounding rock.The indicator bar represents the probability of Classes I-V surrounding rocks.

In the ANN structure,a good prediction of the surrounding rock classification was achieved by extending the input data to multiple parameters (i.e.P-and S-wave velocities,density,Poisson’s ratio,and elastic modulus).These parameters are selected based on the personal experience in engineering,in which Poisson’s ratio and elastic modulus depend on conventional empirical formulae,and whether more key rock properties or mechanical parameters(rock mass integrity factor,shear modulus,etc.) need to be added is worth further exploration.In addition,the network is trained in a supervised manner and relies on the quantity and quality of field borehole data.This network is poorly applicable to places where no borehole data are available.The inclusion of physical laws is an effective means to improve the generalizability of the network(Song et al.,2021).In this regard,how to incorporate the multiparameter physical laws of petrophysics to improve the prediction accuracy of the network under different geological conditions is a hot issue for future research.

A forward-prediction method for tunnel geology and classification of surrounding rock based on seismic wave velocity layered tomography was developed,which can identify the unfavorable geological structure in front of the tunnel face and correct the classification of surrounding rock during the construction,thereby facilitating safe tunnel excavation.To address the strong multisolution of tomography caused by few ray paths in the narrow tunnel space,we presented a regularization-based layered inversion strategy.By reducing the wave velocity inversion area of each iteration step and applying straight-line interface assumption,the convergence and accuracy of inversion are effectively improved.Furthermore,we proposed the surrounding rock classification network based on the autoencoder,and established the mapping relationship between the wave velocity results and the classification of surrounding rock using more rock property parameters as links,so as to improve the prediction accuracy and credibility of the surrounding rock classification.Numerical examples suggest that layered tomography can yield inversion results with clear structure and accurate wave velocity.A detailed field case study demonstrates the potential application of the proposed estimation method for tunnel geology and classification of surrounding rock in the future.

Fig.16.(a) Classification of surrounding rock at the design stage,and (b) Corrected classification of surrounding rock at the construction stage.

Fig.17.Rock slag during TBM excavation.

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 No.51922067),The Key Research and Development Plan of Shandong Province of China(Grant No.2020ZLYS01),and Taishan Scholars Program of Shandong Province of China (Grant No.tsqn201909003).

Appendices A and B.Supplementary data

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

推荐访问:Classification surrounding ROCK