Three-dimensional inversion analysis of transient electromagnetic response signals of water-bearing abnormal bodies in tunnels based on numerical characteristic parameters

1.   Introduction

The transient electromagnetic inversion method solves the physical parameters of the geoelectric section from the measured electromagnetic field on site and forwards modeling results, and it also transforms from a simple resistivity method to a modern and intelligent technology ^[1,2]. In recent years, tunnel defects have become increasingly prominent, and they are largely linked to underground water. Transient electromagnetic methods for tunnel defect detection are being proposed and continuously verified ^[3,4]. Due to the different tunnel environments, the problem of transient electromagnetic inversion is also different from that of other application fields. It is essential to carry out the inversion investigation in the three-dimensional full-space field.

In the field of transient electromagnetic inversion, a large number of geophysicists have carried out relevant investigations. Many methods are used for inversions, such as the Occam inversion, damped least squares method, fastest descent method, Newton method ^{[5,6,7,8,9,10]} and conjugate gradient method ^{[11,12,13,14,15]}. Christiansen et al. performed a nonlinear 1-dimensional inversion based on synthetic aviation data and field data ^[6]. Chang et al. researched the influence of the initial value on the one-dimensional inversion of experimental data ^[7]. Yogeshwar et al. reviewed the structural characteristics of two-dimensional conductivity cross-sections through the inversion of data on site ^[8]. Haber et al. proposed a modified Gauss-Newton strategy to recover a 3-dimensional distribution of electrical conductivity ^[9]. Zhao et al. adopted the nonlinear conjugate gradient inversion in magnetotelluric applications to predict the resistivity distribution based on original measurement data ^[10]. Zhou et al. studied the three-dimensional forward and inverse problems of transient electromagnetic data in full space ^[11]. Li et al. predicted the water-bearing cave in tunnels based on the inversion method of three physical parameter iteration ^[12]. Sun et al. investigated the 3-dimensional inversion data of geological anomalies based on the optimal transport of different parameters of the research target ^[13]. Wang et al. proposed a 3-D parallel inversion algorithm for mine transient electromagnetic method data based on the nonlinear conjugate gradient method and the moving footprint technique based on a large number of survey data ^[14]. Liu et al. adopted TEM data with distances of 50 and 200 m to carry out single inversion based on the Gauss-Newton method, and the inversion result realized the detection of geological bodies with complex terrain and geological structure ^[15]. In recent years, artificial intelligence has developed rapidly, and intelligent algorithms are also being used in transient electromagnetic inversion. Tan et al. used the database inversion expert system, neural network expert system and LSSVM expert system to predict the locations of the abnormal bodies in front of the tunnel face ^[16]. Ertan et al. used the particle swarm algorithm of a global search to quantitatively explain the potential transient electromagnetic data ^[17]. Wang et al. used the artificial bee colony algorithm to calculate the optimal global solution of the multi-layer geoelectric model inversion ^[18]. Liang et al. applied the PSO algorithm to the LSSVM model, and this method effectively improved the accuracy of tunnel detection inversion. Based on the above investigations, the inversion dimension of transient electromagnetic data has gradually developed from one-dimensional to three-dimensional, and the method has progressively changed from a single method to a joint inversion of multiple forms ^[19]. In terms of tunnel defect detection, favorable detection results have been achieved, and this will also provide an essential basis for transient electromagnetic three-dimensional inversion.

Based on the numerical characteristic parameters, a three-dimensional iterative inversion method is proposed in this paper. The correlation coefficients between characteristic parameters and inversion targets were fully investigated. Meanwhile, BP neural network inversion was also added. According to the inversion results of the two inversion methods, the advantages and disadvantages were analyzed. The results provide an essential idea for the subsequent interpretation of the three-dimensional spatial response signal of water-bearing abnormal bodies.

2.   Extraction and correlation analysis of numerical characteristic parameters

2.1.   Extraction of numerical characteristic parameters

The transient electromagnetic response curves of the numerical models were calculated by COMSOL Multiphysics (Version 5.4, AC/DC module). The physics parameters of the surrounding rock medium and the coils can be referenced in previously published literature ^[3,12,21]. Therefore, the transient electromagnetic response signals can be obtained by the numerical simulation model, and the results are demonstrated in Figure 1.

As illustrated in Figure 1, the response signals of the full-space surrounding rock in the case of double logarithmic coordinates are generally linear, and the variation law is consistent with the calculation results in many application fields of TEM ^[12,20]. The detailed physics parameters of the water-bearing abnormal bodies are as follows: The distance between the water body and the coils equals 4 m, the radius of the water body equals 4 m, the thickness of the water body equals 3 m, and the resistivity of the water body equals 5.0 Ω·m. The induced voltage response curve of the water-bearing abnormal bodies has a sudden change in the middle stage, as shown in Figure 1. In addition, this curve has distinct characteristic points, including the uplift point, the peak point and the regression point.

Based on the above comparison, the response processes of the two models are different. To study the difference quantitatively, Li et al. extracted characteristic parameters using an experience equation, but the method has a certain degree of randomness ^[12]. The response curve of the full-space surrounding rock model in Figure 1 shows a linear variation. To better investigate the response of the water-bearing abnormal bodies, a mathematical approach is applied, and the mathematical method of the water-bearing abnormal bodies is obtained as Eq (1):

where ${t}$ means the transient electromagnetic response time, $NC\left(t\right)$ represents the ratio of $M\left(t\right)$ to $N\left(\mathrm{t}\right)$, $M\left(t\right)$ denotes the induced voltage of the water-bearing abnormal bodies at a certain moment, and $N\left(t\right)$ is the induced voltage of the full-space surrounding rock at a certain moment.

The linear transformation result is exhibited in Figure 2. Through the linear transformation of Eq (1), seven parameters can be extracted, and the extraction process is obtained as Eq (2):

where $MEM$ represents the extraction equation of the seven typical parameters, $\mathrm{U}\mathrm{T}$ means the uplift time, $\mathrm{M}\mathrm{V}$ represents the maximum peak value, $\mathrm{M}\mathrm{T}$ denotes the time corresponding to the maximum peak value, $\mathrm{R}\mathrm{T}$ represents the regression time, $\mathrm{L}\mathrm{T}\mathrm{S}$ represents the time span on the left, $\mathrm{R}\mathrm{T}\mathrm{S}$ means the time span on the right, and $\mathrm{T}\mathrm{T}\mathrm{S}$ denotes the time span overall.

Equation (3) shows the mathematical extraction equation of the characteristic parameters. Seven numerical characteristic parameters are used, and four key points can reflect the entire response curve well. The detailed explanation of the numerical characteristic parameters can be found in the work by Qian et al. ^[21].

where $ϴ$ represents the condition for estimating $\mathrm{U}\mathrm{T}$, ${\mathit{\Psi}}$ means the condition for judging $\mathrm{R}\mathrm{T}$ of the characteristic parameter curve, ${t}_{1}$ denotes the time of $\mathrm{U}\mathrm{T}$, ${t}_{2}$ means the time to reach $\mathrm{M}\mathrm{V}$, ${t}_{3}$ represents the time of $\mathrm{R}\mathrm{T}$, and $\mathrm{l}\mathrm{o}\mathrm{g}$ means the logarithm of the time with 10 as the base.

2.2.   Introduction of Pearson's correlation coefficient

Pearson's correlation coefficient reflects the similarity degree between two variables, and it is also introduced to judge the correlations between numerical characteristic parameters and physical parameters of the water-bearing abnormal bodies in the transient electromagnetic detection ^[13]. The calculation of the correlation coefficient is obtained as Eq (4).

where $X$ and $Y$ are defined as two random variables, $\rho (X, Y)$ is the correlation coefficient between the random variables, $\mathit{cov}\left(X, Y\right)$ is the covariance between the random variables, and ${\sigma }_{X}$ and ${\sigma }_{Y}$ are the standard deviations of random variables $X$ and $Y$, respectively.

Based on the statistical analysis module in MATLAB software (Version 2021b), three parameters are selected as the variable objects. Data samples are shown in Table 1, and calculation results are shown in Table 2.

2.3.   Correlation analysis of multifactor parameters

The greater the absolute value of the correlation coefficient is, the stronger the relationship between the two variables. Among the characteristic parameters, $\mathrm{U}\mathrm{T}$ has a very strong correlation with $L$, a weak correlation with $R$ and a very weak correlation with $H$.$\mathrm{R}\mathrm{T}$ has a weak correlation with $L$, a strong correlation with $R$ and a strong correlation with $H$.$\mathrm{M}\mathrm{T}$ has a very weak correlation with $L$, a strong correlation with $R$ and a strong correlation with $H$.$\mathrm{M}\mathrm{V}$ has a strong correlation with $L$, a medium correlation with $R$ and a weak correlation with $H$.$\mathrm{L}\mathrm{T}\mathrm{S}$ has a strong correlation with $L$, a weak correlation with $R$ and a medium correlation with $H$.$\mathrm{R}\mathrm{T}\mathrm{S}$ has a strong correlation with $L$, a very weak correlation with $R$ and a weak correlation with $H$.$\mathrm{T}\mathrm{T}\mathrm{S}$ has a strong correlation with $L$, a weak correlation with $R$ and a medium correlation with $H$.

3.   Three-dimensional inversion based on the iterative principle

3.1.   Proposal of the iterative direction algorithm

Based on the theory of iterative inversion ^[12], the iterative inversion of transient electromagnetic data obtains the physical parameters of water-bearing abnormal bodies based on the response curve of the water-bearing abnormal bodies. First, an initial transient electromagnetic physical model is determined, and the response curve of the model is calculated. Then, the difference between the calculated value and the observed value is used to modify the initial model, and the modified model is run again. Furthermore, the model is modified based on the comparison result, and the iterative process is repeated until the difference (or mean square error) reaches the preset accuracy. Finally, the inversion result of the water-bearing abnormal bodies is obtained when the difference meets the requirement of accuracy.

The critical issue of iterative inversion is the choice of model modification method, which will directly affect the speed of iterative convergence and the correctness of the results. In this paper, numerical characteristic parameters of the water-bearing abnormal bodies are taken as the comparison object. The differences between the numerical characteristic parameters of different models are used to determine the update directions and correction values of the following iteration parameters and to make the next iteration closer to the actual calculated values.

Among the seven numerical characteristic parameters, the first four parameters independently describe the rule of the response curve of water-bearing abnormal bodies, and the last three parameters can be converted from the first three parameters. Simultaneously, the order of correlations between the numerical characteristic parameter $\text{UT}$ and physical parameters of the water-bearing abnormal bodies is as follows: ${L}$ >> ${R}$ > ${H}$. The order of correlations between the numerical characteristic parameter $\mathrm{M}\mathrm{V}$ and physical parameters of the water-bearing abnormal bodies is as follows: ${L}$ > ${R}$ > ${H}$. The order of correlations between the numerical characteristic parameter $\mathrm{R}\mathrm{T}$ and physical parameters of the water-bearing abnormal bodies is as follows: ${H}$ = ${R}$ > ${L}$. The order of correlations between the numerical characteristic parameter $\mathrm{M}\mathrm{T}$ and physical parameters of the water-bearing abnormal bodies is as follows: ${R}$ = ${H}$ >> ${L}$ (>> , > and = , respectively, represent that the correlation distribution interval between the two variables is very large, generally large and none). Therefore, the numerical characteristic parameter $\mathrm{U}\mathrm{T}$ can be used to invert parameter $L$, the numerical characteristic parameter $\mathrm{M}\mathrm{V}$ can be used to invert parameter $R$, and the numerical characteristic parameters $\mathrm{R}\mathrm{T}$ and $\mathrm{M}\mathrm{T}$ can be used to invert parameter $H$.

Based on this, an approximate linear recurrence equation is used to solve ${L}$, ${R}$ and ${H}$. The iterative algorithms are obtained as Eqs (5)–(7).

where ${L}_{i}$, ${R}_{i}$ and ${H}_{i}$ are the distance, radius and thickness of the water-bearing abnormal bodies of the i-th iteration, respectively. ${\mathrm{U}\mathrm{T}}_{input}$, ${\mathrm{M}\mathrm{T}}_{input}$, ${\mathrm{R}\mathrm{T}}_{input}$ and ${\mathrm{M}\mathrm{V}}_{input}$ are input parameters of the iterative inversion. $j$ and $k$ are the step length control parameters of distance $L$ and radius ${R}$, respectively. ${m}_{1}$ and ${m}_{2}$ are the step length control parameters of thickness $H$. $l$ is the lag factor.

The step length control parameters affect the computational efficiency of the whole iterative process. When the step length control parameter is too large, the optimal solution is easily missed. When the step control parameter is too small, the iteration speed is very slow, and the number of iteration steps is too large. The lag factor $l$ mainly controls the iterative process of radius ${R}$ and thickness $H$ to ensure that the distance $L$ can meet the iterative accuracy first, and its value depends on the difference between ${\mathrm{U}\mathrm{T}}_{input}$ and ${\mathrm{U}\mathrm{T}}_{\mathrm{i}-1}$. The larger the difference is, the closer the lag factor is to 0. The smaller the difference is, the closer the lag factor is to 1. After trial calculations of a large number of sample data, the step length control parameter of distance $L$ is determined to be 2, the step length control parameter of radius ${R}$ is to be 2 under the inversion model and parameter conditions, and the two step length control parameters of thickness $H$ are 0.2. The value range of lag factor $l$ is from 0 to 1, which is similar to the activation function of the $sigmoid$ function or $tanh$ function. The value of the lag factor can be selected by referring to the relevant literature ^[12]. The relative error of the inverted numerical characteristic parameters also affects the iterative process. Through multiple trial calculations of sample data, it is determined that the allowable relative error of the numerical characteristic parameters of the transient electromagnetic inversion is 5%. To facilitate the iteration of different inversion objects, the lag factor $l$ is assumed to be 0 before the distance $L$ meets the accuracy. After the distance $L$ meets the accuracy, the value is considered to be 1.

3.2.   Inversion process and code implementation

In the process of transient electromagnetic iteration, partial iteration is used for inversion iteration. According to the correlation relationship, distance $L$ and numerical characteristic parameter $\text{UT}$ have the greatest correlation, so the distance $L$ is first iterated according to Eq (2). Combined with the lag factor $l$, the radius ${R}$ and thickness $H$ are simultaneously iterated according to Eqs (3) and (4) after the relative error of the iteration of the distance $L$ reaches the accuracy requirement. After the iteration parameter meets the relative error, the iteration process of the corresponding parameter can be stopped, and the iteration process of the following iteration parameter can be carried out by analogy. The flow chart of the transient electromagnetic iteration process is shown in Figure 3. The main logic flow used in the three-dimensional iterative inversion program is represented in Table 3.

4.   Inversion analysis of transient electromagnetic response signals

4.1.   Selection of virtual examples based on numerical simulation

Based on the distribution ranges of $\mathrm{L}$, $\mathrm{R}$ and $\mathrm{H}$, five groups of virtual data were regularly selected. The five groups of virtual example samples are shown in Table 4. At the same time, the transient electromagnetic response signals under different physical parameters can be calculated ^[21,22,23]. The numerical characteristic parameters of the water-bearing abnormal bodies are extracted according to Eqs (2) and (3), and the corresponding results are described in Table 4.

4.2.   Inversion result of response signals

Based on the virtual example sample in Table 4, the final inversion result is obtained by iterative inversion of numerical characteristic parameters, as shown in Table 5. According to the iterative inversion results, different types of relative error can be obtained by comparing the theoretical value and the inversion value. For the inversion of distance, the minimum relative error can reach 1.35%, the maximum relative error can reach 4.34%, and the average relative error is 2.55%. For the inversion of the radius, the minimum relative error is 3.31%, the maximum relative error is 13.04%, and the average relative error is 9.08%. For the thickness inversion, the minimum relative error can reach 1.23%, the maximum relative error can reach 9.76%, and the average relative error is 6.38%. In general, the relative error of the iterative inversion results is relatively small, and the accuracy of the prediction results is relatively high, so the accuracy of the iterative inversion results is very high.

Taking virtual sample 4 as an example, the iterative process is analyzed, and the results are shown in Table 6. The variation process of the response curves of the iterative parameters in different iteration steps is shown in Figure 4.

As illustrated in Table 6 and Figure 4, each iteration gradually approaches the target value (the known transient electromagnetic response curve). When the iteration process reaches Step 9, the relative error between the numerical characteristic parameters of this step and the numerical characteristic parameters of the target value is lower than the accuracy tolerance value, so the iteration parameters of Step 9 are taken as the final inversion result.

4.3.   Inversion method of neural network

A neural network is a modeling technology for abstracting a model, from data, that is suitable for solving problems involving intelligent learning ^[24,25]. The neural network can study the internal laws of problems through learning and realize an approximation of the issues through repeated iteration. These characteristics give it unique advantages for solving nonlinear problems and complex influencing factors ^[21,25]. This paper will combine neural network theory to carry out inversion analysis of transient electromagnetic detection data.

A BP neural network is a kind of mapping relationship that can learn and store many input-output patterns without establishing a mathematical model of the system ^[24]. Its learning rule is to adopt the fastest descent method. The weight and threshold of the network are constantly adjusted by comparing the predicted value with the target value after repeated back propagation, so as to minimize the sum of squares of error of the network. The corresponding topological structure of this model includes an input layer, a hidden layer and an output layer. The characteristic mathematical parameters of the water-bearing abnormal bodies are selected as input values, and the physical parameters of the water-bearing abnormal bodies are chosen as output values. The transient electromagnetic BP neural network structure is in the form of one input layer, one hidden layer and one output layer. A schematic diagram of the structure is shown in Figure 5. The number of neurons in the input layer is seven, the number of neurons in the output layer is four and the number of neurons in the hidden layer is fifteen.

Using the self-learning mathematical method of the BP neural network, 125 groups of data results are used as samples, the output value is repeatedly compared with the ideal value until the accuracy reaches a specific range, and then a trained data system is obtained. The transient electromagnetic response curve under the characteristic condition of unknown water-bearing abnormal bodies can be obtained through on-site, actual indoor measurement or virtual data. The numerical characteristic parameters of the response curve can also be obtained through normalization processing, and they are used as input values and put into the successfully trained data system. Furthermore, the physical parameters of the water-bearing abnormal bodies can be obtained, including the distance, radius and thickness.

Similarly, based on the virtual example samples in Table 4, the trained data system is used, seven numerical characteristic parameters are input and then the corresponding three physical parameters can be obtained. Meanwhile, the correlation error between the theoretical value and the predicted inversion value is calculated, and the results are represented in Table 7.

According to the BP neural network inversion results of the five calculation examples in Table 7, the minimum relative error can reach 0.32%, the maximum relative error can reach 9.09%, and the average relative error is 5.48% for the inversion of distance. For the inversion of the radius, the minimum relative error is 2.47%, the maximum relative error is 12.96%, and the average relative error is 9.53%. For the thickness inversion, the minimum relative error can reach 4.39%, the maximum relative error can reach 10.98%, and the average relative error is 6.94%. In general, the relative error of the BP neural network inversion result is minor, and the accuracy of the prediction result is relatively high, so the accuracy of the three-dimensional inversion result based on this method is better.

4.4.   Comparison of the two methods

Based on the calculation processes and results of the above two inversion methods, the advantages and disadvantages of different inversion methods are comprehensively evaluated with algorithm complexity, calculation speed, inversion accuracy, updated sample data and hardware requirements as the analysis objects. The inversion comparison results are shown in Table 8.

As shown in Table 8, the two inversion methods have advantages and disadvantages in different performances. In terms of algorithm complexity, the iterative inversion algorithm is more complicated. In contrast, the BP inversion is relatively uncomplicated. In terms of calculation speed, iterative inversion requires multistep iteration, so the calculation speed of iterative inversion is slow. However, the BP inversion is based on prediction or query parameter inversion, so the calculation speed of BP inversion is faster. In terms of inversion accuracy, the iterative inversion accuracy is relatively higher by comparing the average relative errors of the two methods, and the BP inversion accuracy is relatively lower under the condition of a simple neural network structure. In terms of updated sample data, the iterative inversion and the BP inversion cannot update database samples automatically. For hardware requirements, the iterative inversion needs to use the COMSOL simulation software for the forward calculation in the subsequent iterations and the MATLAB software to obtain the numerical characteristic parameters, so the iterative inversion has high hardware requirements. BP inversion only uses MATLAB software to realize parameter inversion based on intelligent parameter learning, so the BP inversion has relatively low hardware requirements.

In summary, both inversion methods can realize the inversion of physical parameters of water-bearing abnormal bodies, and each method has advantages and disadvantages. Therefore, it is necessary to fully utilize the benefits of different inversion methods and use multiple inversion methods to realize the three-dimensional inversion of water-bearing abnormal bodies.

5.   Conclusions and discussion

This paper mainly investigates the three-dimensional inversion of water-bearing abnormal bodies by two different inversion methods. The main work carried out in this paper is as follows:

1) Among the seven numerical characteristic parameters, the first four parameters describe the rule of the response curve of the water-bearing abnormal bodies independently. The last three parameters can be converted from the first three parameters, so only the first four parameters are independent. Therefore, these four parameters are mainly selected as inversion calculation values in the iterative inversion method.

2) The inversion results of the physical parameters of the water-bearing abnormal bodies are obtained using two inversion methods, and the advantages and disadvantages are analyzed. Both inversion methods can obtain the inversion results of the physical parameters of the water-bearing abnormal bodies. However, different inversion methods have various errors in the inversion results of various physical parameters.

3) A single inversion method cannot realize the three-dimensional inversion of the physical parameters of water-bearing abnormal bodies quickly, effectively and intelligently. Therefore, it is necessary to evaluate the three-dimensional inversion of water-bearing abnormal bodies jointly by integrating multiple inversion methods.

In this paper, two different three-dimensional inversion methods were mainly carried out based on virtual data, and the two methods have certain inversion abilities. In future work, the judgment criterion of the iterative inversion and the structural form of the BP inversion need to be further investigated. Meanwhile, the response signal data from the laboratory test and the actual field test will be applied to the joint inversion, which can further improve the accuracy and applicability of the three-dimensional inversion of transient electromagnetic response signals of water-bearing abnormal bodies of tunnel defects.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 52208395), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. KYCX22_3350), the Systematic Project of Guangxi Key Laboratory of Disaster Prevention and Engineering Safety (Grant No. 2021ZDK015), the Nantong Science and Technology Plan Project (Grant Nos. JC2021169, JC2020121, and MS12020080) and the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 18KJB580015). The authors appreciatively acknowledge the financial support of the abovementioned agencies.

Conflict of interest

The authors declare that there are no conflicts of interest.