Tunnel Concealed Karst Cave Joint Detection by Tunnel Seismic and Transient Electromagnetic

The water-rich karst cave is the main hazard source during tunnel construction in the karst area; some tunnel accidents caused by huge water-rich caves have been reported [1]. A giant karst cave was encountered during the excavation of a water diversion tunnel in the karst region of southwestern China. Three pilot holes explore that the karst cave is more than 15 m in the direction of advance. The core samples of the pilot holes show that the filling material inside the cave is mainly underconsolidated flowing plastic yellow-colored soft clay. The direction of advance and the extended boundary on the left and right sides of the cave body is unknown. The giant cave and undiscovered cave will threaten the safety of tunnel construction. The tunnel seismic and transient electromagnetic (TEM) are the primary geophysical tools for tunnel advanced detection [2–4].

The tunnel seismic method can predict the fault fracture zone, weak interlayer, karst cave, and other unfavorable geological bodies in front of the tunnel face by analyzing the seismic signals received by the tunnel wall. Tunnel seismic methods mainly include horizontal seismic profile [5], true reflection tomographic [6], tunnel seismic tomography [7, 8], and tunnel seismic prediction (TSP) [9–11]. Compared with other prediction methods, the TSP method has been widely used in tunnel advanced detection due to its simple construction and higher detection precision. Due to the limitation of tunnel space, the migration profile obtained by the TSP method can only provide the boundary interface between the karst cave and the surrounding rock in the direction of advance but cannot determine the spatial distribution range of the cave body.

TEM method analyzed the primary and secondary field signals received by the coil, computes the apparent resistivity in front of the tunnel face and surrounding area using an inversion algorithm, and predicts unfavorable geological bodies such as filling water karst cave and water-rich fracture zones. Many scholars have made a great advance in terms of instruments, algorithms, and interpretations and achieved many successful cases. Xue and Li [12] proposed to draw an imaging profile based on the secondary conductance differential parameter data by applying the theory of “floating sheet.” Su et al. [13] summarized the interpretation methods, characteristics, and application effects of apparent longitudinal conductance. Sun et al. [14] proposed a tunnel TEM multipoint array detection method to analyze the 3D spatial distribution of water-bearing structures. S. Li et al. [15] proposed the parallel magnetic field response method for tunnel TEM advance prediction. Qi et al. [16] studied the whole-space whole-domain apparent resistivity interpretation method for the tunnel TEM arbitrary coplanar vertical magnetic field. K. Li et al. [17] investigated the effects of various device types on the data acquisition and imaging results. Zhang et al. [18] designed a homogeneous half-space model and a layered model to test the suitability of the TEM response in detecting the conversion of induced electromotive force into apparent resistivity. Xing and Jiang [19] designed a 3D multiturn small loop TEM device for tunnel advance prediction. Previous studies have revealed that the TEM method has a high-resolution to low-resistivity anomalous bodies such as water-rich structures and fissure water.

TSP and TEM have strengths and weaknesses in the detection imaging of karst cave. The TSP migration profile can provide information on the boundary between the cave and the surrounding rock in the direction of advance but cannot determine the extended boundary on the left and right side of the cave body. TEM method has a high resolution for low-resistivity bodies and can detect the spatial range of cave bodies from multiple angles by changing coil direction. The TEM inversion results rely heavily on the initial resistivity model. There are two types of TEM inversion methods: linear inversion and nonlinear inversion methods. Nonlinear inversion methods, such as genetic algorithm [20], simulated annealing algorithm [21], particle swarm optimization [22], and others, start with a random model. The nonlinear inversion method reduces the dependence of the initial model when compared to linear inversion methods such as least square, but the model update is stochastic and requires a lot of forwarding simulation and inversion calculation. The stochastic model may give rise to multisolution problems. An optimal initial model can reduce multisolution, save calculation time, and improve inversion accuracy. The pilot hole data and TSP results are taken as the prior information, the error of the initial model can be reduced, and the inversion accuracy of the karst cave body can be improved.

Previous studies mainly utilized a single method to detect the karst caves [10, 19], which is insufficient to describe the spatial distribution range of giant karst caves. To overcome the shortcoming of a single method, we developed tunnel seismic and TEM joint tunnel advanced prediction in this paper. To solve the problem of the initial resistivity layered model in TEM inversion, this paper proposes a generation method of an initial model: (1) the layer’s thicknesses are from the tunnel seismic migration data; (2) the resistivity values are from the pilot holes. The pilot hole data, TSP, and TEM data processing steps are described in detail, and a complete TSP data and TEM data processing flow are established. The inversion results of different initial models are compared to prove the rationality of the proposed initial model method. The TSP and TEM results will be presented and discussed to demonstrate the advantages of joint detection.

The water diversion tunnel is located in Chongqing, southwest China. The age of the surrounding rock stratum is the Maokou Formation of the Lower Permian (P₁m). The lithology is dominated by limestone, with medium-thick stratiform. The attitude of strata is 290°-300°<32°-38°. The tunnel buried depth is from 120 m to 130 m. The surrounding rock has fewer fissures. According to the classification of surrounding rock integrity, the integrity of the surrounding rock belongs to integrity and is the highest. According to the classification of surrounding rock type, the rock belongs to III type. The estimated resistivity value of the surrounding rock is 680 Ω•m. The tunnel face is 2.2 m wide and 2.62 m high. The results of the geological description of the three pilot holes are shown in Table 1. Studies conducted in the same environment demonstrated that the soft plastic filling material in the gigantic karst cave acts as a highly effective attenuation medium for seismic and electromagnetic waves [23, 24]. The modest current TEM equipment, as well as the low-energy seismic source, are having difficulty penetrating the massive karst cave. Tunnel seismic equipment necessitates the use of a powerful energy source and a geophone with great sensitivity [25, 26]. Coils with high sensitivity and a considerable current are required for TEM devices [27].

According to the geological description of the pilot holes, the spatial distribution of the giant karst cave can be inferred as shown in Figure 1(b). The theory of 1D TEM inversion is based on the layered stratum, so we need to regard the surrounding rock and karst cave as a virtual layered stratum. Figure 1(c) shows the relationship between the concept of layer and the surrounding rock, cave, and the direction of advance. The direction of the layer is perpendicular to the direction of tunnel excavation. The symbols T1 to T4 are the spatial interfaces of the karst cave. The surrounding rock and karst cave are not layered stratum; too many layers are unreasonable. According to the three borehole data, the medium in front of the tunnel face can be divided into three layers. The surrounding rock between the tunnel face and the interface L1 is the first layer, named layer 1. The cave is the second layer, named layer 2. The surrounding rock in the front of the cave body is the third layer, named layer 3. The concept of layer thickness is also shown in Figure 1(c).

This paper proposes TSP and TEM joint detection methods and establishes a complete data processing flow. An optimal initial model can save iteration time, reduce multisolution, and obtain high-precision inversion results. The initial apparent resistivity model is established based on TSP and pilot hole data. The detailed procedures are shown in Figure 2. The main data processing steps are as follows:

According to tunnel site conditions and the principle of the tunnel seismic method, an optimum tunnel seismic geometry is designed. Twenty-four sources were placed on one side of the tunnel wall. The depth of the source holes was 2.0 m. The shot interval was 1.5 m. The depth of the geophones borehole was 1.5 m. The last shot is 5 m away from the tunnel face. Two geophones are used in the conventional TSP geometry. Increasing the number of geophones can enhance the illumination of the target bodies. Four three-component acceleration geophones are used in this paper (Figure 3, R-1 to R-4). The bandwidth is between 0.1 Hz and 5000 Hz [25]. The distance between R-2 and R-4 and between R-1 and R-3 is 20 m. The installation of the geophone is shown in Figure 4(a). The TETSP-2 tunnel seismic advanced detection system was used to collect seismic data (Figure 4(b)). There are obvious differences between TSP and TETSP-2 devices in geophone installation. The TETSP-2 device is used a seismic-automatic coupling geophone that can be installed in a few minutes [25]. The source is a 30 g explosive source. The sampling interval is 41.7e^-6 s. The number of samples per trace is 8,192.

The raw seismic data has a high signal-to-noise ratio (SNR). Figures 5 and 6 show the three-component seismic data processed by AGC at R-1 and R-3, respectively. Figures 7 and 8 are locally enlarged waveforms from 0 ms to 30 ms. Wavefield identification is mainly based on slope and amplitude. The direct wave is a P-wave with the highest velocity. The surface wave and converted S-waves are slower than P-wave, their time is greater than P-wave, and their energy is stronger. The velocity of the acoustic wave is 340 m/s. The numerical simulation and measured data by several scholars [26, 28] also prove that the tunnel seismic data contain surface waves. According to these rules, it is considered that symbol A denotes the direct P-wave, B the surface wave, and C the acoustic wave.

The main factors affecting the SNR of seismic data at R-1 and R-3 are geophone sensitivity, geophone coupling, and offset. The sensitivity of all three-component geophones is the same. The coupling degree of two geophones may be different, but the coupling difference is so small that it can be discarded. In general, the lower the offset, the less reflection wave attenuation and the higher the SNR. The surface wave with strong energy can suppress random noise and improve SNR. The larger the offset, the stronger the surface wave and the higher the SNR. The lower the offset, the stronger the acoustic wave and the lower the SNR. The distribution of surface waves varies with different components [29]. The surface wave is the most obvious in $Y$-component (Figures 7(b) and 8(b), B).

There are surface waves from 0 ms to 30 ms and acoustic waves from 30 ms to 200 ms in Figures 5 and 6. To analyze the influence of surface wave and acoustic wave on SNR, the SNR of seismic data in different time ranges is calculated, and the results are shown in Table 2. The offset of R-1 is less than that of R-3; the SNR of $X$- and $Z$-component data at R-1 is greater than that of R-3. For the $Y$-component surface wave from 0 ms to 30 ms, R-3 (Figure 8(b), B) is more obvious than R-1 (Figure 7(b), B), especially in region B1. Table 3 shows the statistical results of SNR from 0 m to 9 m and 10.5 m to 36 m in 0 ms to 30 ms. The SNR of data at R-1 is lower than that of R-3, and the slight difference is between 10.5 m and 36 m, while the significant difference is between 0 m and 9 m (Table 3). There are reflection waves and noise from 30 ms to 40 ms. The offset of R-1 is lower than that of R-3, and the SNR of data at R-1 is greater than that of R-3. There are reflection waves, acoustic waves, and noise from 30 ms to 200 ms, and the SNR of $Y$-component data at R-1 is greater than that of R-3. The direction of acoustic propagation in the tunnel is consistent with the direction of the $Y$-component, so the $Y$-component has the strongest acoustic interference, while the $X$- and $Z$-component have less influence. Compared with the range from 30 ms to 40 ms, acoustic waves existing from 30 ms to 200 ms will reduce the SNR of the three-component data. By comparing their SNR values, it is found that the effect of the acoustic wave on the SNR of $X$- and $Z$-components is small, but that of the $Y$-component is significant. Based on the above analysis, the SNR of the $Y$-component at R-1 is lower than that of R-3 from 0 ms to 200 ms; the main reason is the influence of the surface wave, especially region B1.

The spectrum analysis results within 0 ms to 60 ms show that the main frequency bandwidth of data at R-1 is from 50 Hz to 1500 Hz (Figure 9(a)), while that of data at R-3 is from 50 Hz to 1000 Hz (Figure 9(b)). The filtering frequency bandwidth is used from 50 Hz to 800 Hz in the data processing. With the propagation distance increased, the amplitude energy of the seismic wave decreases continuously, and the amplitude value at R-3 is lower than that at R-1 (Figure 9).

The raw seismic data at R-1, R-2, R-3, and R-4 were processed by the RMS, AGC, band-pass filtering, F-K transformation, predictive deconvolution, polarization filtering, etc. [26]. The high SNR P-wave data from the direction of the advance can be obtained. The continuous reflected events with high amplitude in the wavefield are considered as the reflective wave from the direction of advance (Figure 10, the red-dotted lines).

The first arrival travel time of $Y$-component data at four geophones was picked successively. The direct wave velocity was obtained by least square fitting, and the average P-wave velocity of the surrounding rock is 6288 m/s. The average P-wave velocity is used as the migration velocity model; equal-time plane migration imaging (Figure 11) is achieved. Stack the migration data at R-1 and R-2 to obtain the migration profile in Figure 11(a). Stack the migration data at R-3 and R-4 to obtain the migration profile in Figure 11(b). The migration profiles (Figure 11) show that the position of the strong reflection interfaces is highly similar. The stack processing can suppress interference and highlight the interface of a real anomalous body. Figure 11(c) is the stack processing result of Figures 11(a) and 11(b). The real reflection interfaces (A, B, and C) are enhanced, while the pseudoanomalous reflection interfaces (A1 to A4) are suppressed. There are some strong reflection interfaces at mileage K0+594 (A), K0+629 (B), and K0+669 (C). Combined with the three pilot hole data (Table 1), it is preliminarily predicted that the scale of the giant cave is from K0+594 to K0+629 m, and K0+669 is the interface of the second cave. The interfaces T1, T2, and T5 can be obtained from the migration profiles.

According to the proposed data processing method of tunnel seismic data and TEM, the layer’s thicknesses are extracted from the tunnel seismic migration profile (Figure 11(c)). The initial resistivity values are determined based on both borehole data and prior geological data. The background resistivity value of the tunnel rocks (limestone) is given as 680 Ω·m. The resistivity value of the karst cave is given as 520 Ω·m. In this paper, the TEM inversion is 1D single-point inversion. The initial 1D resistivity layered model is shown in Figure 12. Figure 12 reveals that the relatively low-resistivity areas are at the mileage from K0+599 to K0+629 (B), and the relatively high resistivity areas are at other mileage. The average apparent resistivity model (650 Ω·m) and the constraint apparent resistivity model were used as the initial model for TEM inversion.

According to tunnel site conditions and the principle of the TEM detection method, the TEM geometry was designed. A horizontal survey line (Figure 13(a), L1) and a fan-shaped scanning survey line (Figure 13(a), L2) were designed on the tunnel face. To investigate the left and right extension of the giant karst cave, fan-shaped scanning survey lines were designed on both sides of the tunnel wall 15 m away from the tunnel face (Figure 13(b), L3 and L4). The FCTEM60-1 device was used to collect TEM data. The main technologies of the FCTEM60-1 device include “constant voltage clamping” fast linear turn-off technology and lossless mutual inductance elimination technology. The device has a variety of coils and has been used in domestic engineering fields in China. The coil structure is discussed in detail in the literature [27]. The distance between the two observation points is 0.5 m, the power supply voltage is 21 V, and the transmitting current is 60 A. Turn-off time is 70 μs. The shielding coil diameter is 0.9 m, and the transmitting magnetic moment of the coil is 1500 Am². The coil was placed on the tunnel face and the tunnel left wall as shown in Figure 14.

To verify the effect of the different initial model on TEM inversion results, we take L1 line as an example. The average apparent resistivity model and constraint resistivity (Figure 12) are used as the initial model for inversion calculation, respectively. Figure 15 shows that the inversion results have significant differences, especially the mileage that goes from K0+639 to K0+689. The average apparent resistivity model was used as the initial model; the inversion profile as shown in Figure 15(a) was obtained, with three low-resistivity anomalous areas (A, B, and C). When the constrained apparent resistivity model was used as the initial model, the inversion profile as shown in Figure 15(b) was obtained, with only one low-resistivity anomalous area (B). The low-resistivity areas are from K0+594 to K0+635, while the relatively high-resistivity areas are from K0+649 to K0+ 689 (Figure 15(b)). The pilot hole data show that the mileage of the surrounding rock is from K0+ 589 to K0+594, and the giant karst cave is from K0+594 to K0+ 601. The TEM method has a volume effect; the resolution will decrease with detection distance. Figure 15(b) is more consistent with this theoretical basis than Figure 15(a). Combined with the boreholes data and the theoretical basis of the TEM method, Figure 15(b) is closer to actual geological characteristics. The comparison results in Figure 15 demonstrate the rationality of the proposed initial model.

The inversion calculations are similar to Figure 15; the inversion profiles of the L2 line are obtained when using the average resistivity model and constrain resistivity model as the initial model. For the L3 and L4 lines, we cannot obtain constrained resistivity models due to a lack of TSP and borehole data. Therefore, the average apparent resistivity model is used as the initial model for TEM inversion. The inversion profiles of the L1 to L4 lines are pieced together; the final profile is shown in Figure 16. When the average resistivity model was the initial model, Figure 16(a) shows three low-resistivity anomalies (A1, A2, and A3). When the constrain resistivity model was the initial model, Figure 16(b) shows a low-resistivity area (A6 and A7). The inversion profiles of the L1 and L2 lines are consistent in the coincidence area (Figure 16(b), red arrow). Combined with inversion profiles in different directions, Figure 16(a) shows multiple low-resistivity areas (A1 to A5). The inversion profiles of the L2, L3, and L4 lines are inconsistent in the adjacent areas (Figure 16(a), red circles). The low-resistivity anomaly area of the giant karst cave can be visually displayed, and the boundary of the cave (dark circle) can be determined in Figure 16(b). Compared with the boreholes data and TSP migration profile, Figure 16(b) is closer to the actual karst cave.

According to the advanced drilling data, the TSP migration profiles (Figure 11(c)), and the TEM inversion results (Figure 16(b)), the spatial distribution of the giant karst cave is delineated (Figure 17). The mileage of the cave in the direction of advance is from K0+594 to K0+ 630. The distance in the direction of advance is about 36 m. The karst cave extends 28 m to the left side and 43 m to the right side. The TSP migration profile shows a strong impedance interface at the mileage K0+669, so it is speculated that the mileage from K0+669 to K0+ 684 may be a karst body (Figure 11(c)). The interfaces T1 to T4 are shown in Figure 16(b).

To enhance the energy of the anomaly body in the migration profile, two geophones were added to the traditional TSP geometry. There are four geophones in total in this paper (Figure 3). The weighted stacking method was adopted to stack the migration data of all geophones, suppress the interference wave, and improve the energy of the actual anomalous interface (Figure 11(c)). A 2D linear TSP geometry was designed to accommodate the limitation of tunnel space. The acquired three-component seismic data come from the tunnel 3D space. We assumed that after the seismic data processing (Figure 2), reflected waves (Figure 10) come from the observation plane. Under this circumstance, only the position of the strong reflection interface can be extracted from the TSP migration profile, and the boundary interface between the karst cave and surrounding rock (T1, T2, and T5) can be delineated, but the spatial distribution of the karst cave cannot be described.

In this paper, a particle swarm nonlinear inversion algorithm was used to achieve TEM inversion. The average resistivity model was used as the initial model; a lot of forwarding simulation and inversion calculation are needed to search for the best fitting parameters. Parallel computing methods such as CPU and GPU can only improve the computational efficiency but cannot solve the multisolution problem. The same fitting curve can be obtained by forwarding modeling with different layers’ thicknesses and resistivity values. However, the inversion profiles are inconsistent in the adjacent areas (Figure 16(a), red circles). Therefore, it is difficult to determine actual anomalous bodies when the average model is the initial model.

TEM inversion results are very dependent on the initial apparent resistivity model. The ground geophysical data cannot provide a reliable initial apparent resistivity model for the tunnel TEM inversion [30, 31]. This paper proposes a generation method of an initial model: (1) the layer’s thicknesses are from the tunnel seismic migration data; (2) the resistivity values are from the pilot hole data. The optimized initial model can reduce the search area and quickly search the optimal layer’s thickness and apparent resistivity value, which is beneficial to save calculation time and improve the inversion accuracy. With the average model as the initial model, a point takes 1 minute to calculate. While the constraint model is the initial model, the calculation time of one point is only 5 s.

The comparison of TEM inversion profiles shown in Figure 15 proves the rationality of the proposed initial model. The splicing results of the inversion profiles in different directions correspond well, which can visually display the spatial distribution range of the karst cave (Figure 16(b), T1 to T4). We cannot collect sufficient TSP and TEM data due to tunnel space limitations. A small amount of data from a single detection method cannot delineate the range of cave bodies. The results of multiple methods verify each other, and the multidata constraint each other can improve the prediction accuracy. If the results of the pilot borehole, TSP, TEM, and actual geological data are consistent, we can consider the delineated cave to be credible.

To investigate the spatial distribution of the giant karst cave, tunnel seismic and TEM methods were used for tunnel advanced joint detection in this paper. A complete tunnel seismic and TEM data processing flow is established. To solve the problem of the initial resistivity model in TEM inversion, this paper proposes a generation method of an initial model. Based on the imaging results of the giant cave, the following points are summarized:

• (1)

  Four geophones were used to enhance the energy of the real wave impedance interface. The layers’ thicknesses were determined by the impedance interface from the TSP migration data. The resistivity values of different layers are determined from pilot hole data

• (2)

  The constraint apparent resistivity model can reduce the range of low-resistivity anomalous bodies, save calculation time, and improve the precision of inversion imaging

• (3)

  The spatial distribution of the giant cave can be detected by changing the coil angle. The inversion results from all directions can be verified each other, which can improve the precision of the inversion result

It is necessary to accurately delineate the spatial distribution of giant caves. To overcome the shortcoming of a single method, multimethod detection and data processing were carried out. The tunnel seismic results provide prior information for TEM and help to establish an optimal initial apparent resistivity model. The multimethod imaging results can be verified and referenced to improve the precision of the target body.

The seismic and TEM data used to support the findings of this study were supplied by Xinglin Lu under license and so cannot be made freely available. Requests for access to these data should be made to Zhihong Fu (fuzhihong@cqu.edu.cn).

The authors declare no conflict of interest.

This work was supported by the National Key Research and Development Program (2018YFC0406904) and the Science and Technology on Near-Surface Detection Laboratory (6142414200613).