A fast approximation for 1-D inversion of transient electromagnetic data by using a back propagation neural network and improved particle swarm optimization

As one of the most active nonlinear inversion methods in transient electromagnetic (TEM) inversion, the back propagation (BP) neural network has high efficiency because the complicated forward model calculation is unnecessary in iteration. The global optimization ability of the particle swarm optimization (PSO) is adopted for amending the BP's sensitivity to its initial parameters, which avoids it falling into a local optimum. A chaotic-oscillation inertia weight PSO (COPSO) is proposed for accelerating convergence. The COPSO-BP algorithm performance is validated by two typical testing functions, two geoelectric models inversions and a field example. The results show that the COPSO-BP method is more accurate, stable and needs relatively less training time. The proposed algorithm has a higher fitting degree for the data inversion, and it is feasible to use it in geophysical inverse applications.

1 Introduction

The transient electromagnetic (TEM) method applies the secondary receiving voltage induced by the rapid switching of pulse current, and it then deduces the geoelectrical parameters consisting of the resistivities and thicknesses of the layers. The later is a typical TEM inversion issue with nonlinear features. The linear inversion method was simple and widely used through linearization processes, yet it is extremely dependent on the selection of initial parameters and results in poor inversion accuracy. Hence, the nonlinear inversion methods have attracted more geophysicists' attention in recent years.

The artificial neural network (ANN) is one of the most active nonlinear inversion methods, and it has very high computation efficiency because the complicated forward model calculation is unnecessary in this iteration. All the geoelectrical parameters and the forward model relations are implied in the weight and threshold parameters of ANN. And it is different from the nonlinear Monte Carlo method with a global space search solution (He et al., 2018; Jha et al., 2008; Pekşen et al., 2014; Sharma, 2012; Tran and Hiltunen, 2012). Srinivas et al. (2012) compared the inversion performance of back propagation (BP), the radial basis function (RBF) and the generalized regression neural network (GRNN) in vertical electrical sounding data, then established a 1-D inversion model with BP and finally realized the parameter inversion. Maiti et al. (2012) proposed a Bayesian neural-network training method in 1-D electrical sounding. Jiang et al. (2018) improved the training method for the kernel principal-component wavelet neural network and achieved the resistivity imaging. Jiang et al. (2016a) produced a learning algorithm based on information criterion (IC) and particle swarm optimization for RBF network which improves the global search ability. Johnson and Aizebeokhai (2017) utilized neural-network method to invert multi-layer georesistivity sounding. Jiang et al. (2016b) presented a pruning Bayesian neural-network (PBNN) method for resistivity imaging and solved the instability and local minimization problems. Raj et al. (2014) solved nonlinear apparent resistivity inversion problems with ANN. The ANN has been widely applied in electric-prospecting data interpretation for its powerful fitting ability. However, the neural-network method is sensitive to its initial parameter settings and falls easily into a local minimum. Many improved methods were proposed for balancing the convergence rate and inversion quality. Zhang and Liu (2011) proposed ant colony optimization for ANN, and they applied high-density resistivity and acquired smaller inversion errors and higher determinant coefficients. Dai et al. (2014) suggested a differential evolution (DE) for BP which enhanced the global search ability. Rosas-Carbajal et al. (2014) introduced the genetic algorithm for ANN.

Figure 1 Inertial weight curves comparison.

The particle swarm optimization (PSO) has a simple structure, fast convergence rate, high accuracy and global optimization ability. Fernández et al. (2010) successfully introduced the PSO in a 1-D resistivity inversion. Godio and Santilano (2018) applied it in a geophysical inversion and deduced a depth resistivity earth model. Due to the PSO's global searching performance, the BP's initial weights and thresholds can be trained by the PSO, and the BP's global optimization ability can be improved. Compared to the standard PSO (SPSO), a chaotic-oscillation inertia weight PSO (COPSO) can accelerate the convergence rate in the early stage, and it was proposed naturally (Shi et al., 2009).

The paper structure is as follows: the principles of the PSO algorithm with different inertia weights schemes, the BP neural network and the proposed COPSO-BP algorithm are given in Sect. 2. Then, the COPSO-BP algorithm performance is validated by two typical testing functions in Sect. 3. And in a later section, inversion simulations of three-layer and five-layer geoelectric models are carried out; the hidden-layer neuron numbers determining method is put forward; and algorithm performance is compared.

2 Principle of COPSO-BP algorithms

2.1 Chaotic-oscillation PSO algorithm

For an n-dimensional optimization problem, it is supposed that the position (resistivity and thickness for layered-model parameter inversion) and velocity (update speed) of the ith particle (global search group number) at time t are x_i= (x_i1, x_i2, …, x_iN) and v_i= (v_i1, v_i2, …, v_iN), respectively. Then, at time t+1, they can be calculated by the iterations as

$$\begin{array}{}\text{(1)}& {\displaystyle }{v}_{\mathrm{id}}^{t+\mathrm{1}}=\mathit{\omega }\cdot v_{\mathrm{id}}^t+c_{\mathrm{1}}{r}_{\mathrm{1}}(p_{\mathrm{id}}^t-x_{\mathrm{id}}^t)+c_{\mathrm{2}}{r}_{\mathrm{2}}(p_{\mathrm{gd}}^t-x_{\mathrm{id}}^t),\text{(2)}& {\displaystyle }{x}_{\mathrm{id}}^{t+\mathrm{1}}=x_{\mathrm{id}}^t+v_{\mathrm{id}}^{t+\mathrm{1}},\end{array}$$

where r₁ and r₂ are random values evenly distributed in the interval (0,1), c₁ and c₂ are learning factors (usually equal to 2), and p_id and p_gd are the individual and global maximum values.

Figure 2 Three-layer BP neural-network structure.

Figure 3 Training error curves of the SPSO-BP and COPSO-BP algorithms.

The inertia weight parameter ω affects the algorithm performance seriously. A fixed weight always was used in the early time, and then various dynamic weights were proposed. Shi et al. (2010) have summarized several methods as

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{1}}\left(t\right)={\mathit{\omega }}_{\mathrm{s}}-({\mathit{\omega }}_{\mathrm{s}}-{\mathit{\omega }}_{\mathrm{e}})t/T_{max},\text{(4)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{2}}\left(t\right)={\mathit{\omega }}_{\mathrm{s}}-({\mathit{\omega }}_{\mathrm{s}}-{\mathit{\omega }}_{\mathrm{e}})(t/T_{max}{)}^{\mathrm{2}},\text{(5)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{3}}\left(t\right)={\mathit{\omega }}_{\mathrm{s}}-({\mathit{\omega }}_{\mathrm{s}}-{\mathit{\omega }}_{\mathrm{e}})\left[\mathrm{2}t/T_{max}-(t/T_{max}{)}^{\mathrm{2}}\right],\end{array}$$

where ω_s and ω_e are the start and end weight. The t and T_max are the current and maximum iterations. The above weights are smooth and monotonically decreasing. In this paper, we proposed a decreasing oscillation weight scheme which was based on the chaotic logistic equation. Its specific calculation formula is

$$\begin{array}{}\text{(6)}& {\displaystyle }{x}_{t+\mathrm{1}}=\mathit{\mu }{x}_t\left(\mathrm{1}-x_t\right)t=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },n,\text{(7)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{c}}\left(t\right)={\mathit{\omega }}_{\mathrm{e}}+\left({\mathit{\omega }}_{\mathrm{s}}-{\mathit{\omega }}_{\mathrm{e}}\right)\left({\mathrm{0.99}}^t\cdot x_t\right),\end{array}$$

where μ is the control parameter. A complete chaos state is established for x∈(0,1) and μ=4, and an inertia weight is then obtained from Eq. (7). Numerical experiments were carried out correspondingly and showed that the initial value of x₀ has little effect on the inertia weight ω. The inertia weight comparison was shown in Fig. 1, where x₀=0.234 and μ=4 for chaotic oscillation.

2.2 BP neural network

The BP neural network has a multi-layer feed-forward structure, and a typical three-layer network is shown in Fig. 2 (Li et al., 2009).

For Fig. 2, x₁,x₂, …, x_n are the input values, y₁, y₂, …, y_m are the predicted outputs, and w_ij and w_jk are the network weights. The threshold parameter α is defined in the hidden layer with its output,

$$\begin{array}{}\text{(8)}& H_j=f\left(\sum _{i=\mathrm{1}}^n{w}_{ij}{x}_i-{\mathit{\alpha }}_j\right)j=\mathrm{1},\mathrm{2},\mathrm{\dots },l,\end{array}$$

where l is the hidden-layer node numbers, f is the activation function with different expressions, and the most widely used is a sigmoid-type function. The predicted output for the kth unit is calculated by

$$\begin{array}{}\text{(9)}& O_k=\sum _{j=\mathrm{1}}^l{H}_j{w}_{jk}-b_k,\end{array}$$

and parameter b is the output threshold. Then the prediction error can be determined based on the predicted output O_k and the expected output T_k, which is $e_k=(T_k-O_k)O_k(\mathrm{1}-O_k)$. The updated formula for weights and thresholds is the following:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{c}{w}_{ij}=w_{ij}+\mathit{\eta }{H}_j(\mathrm{1}-H_j)x_i{\sum }_{k=\mathrm{1}}^m{w}_{jk}{e}_k\\ w_{jk}=w_{jk}+\mathit{\eta }{H}_j{e}_k\\ {\mathit{\alpha }}_j={\mathit{\alpha }}_j+\mathit{\eta }{H}_j(\mathrm{1}-H_j){\sum }_{k=\mathrm{1}}^m{w}_{jk}{e}_k\\ b_k=b_k+e_k\end{array},\right.\end{array}$$

where i is 1,2, …, n, j is 1,2, …, l, k is 1,2, …, m and η is the learning rate.

2.3 BP neural network with the COPSO algorithm

The initial parameters are chosen randomly, which affects the convergence rate, learning efficiency and perhaps falling into a local minimum. The chaotic-oscillation PSO (COPSO) has a much better global optimization capability; therefore, the COPSO algorithm is proposed to optimize the initial weight and threshold of the BP. The COPSO-BP pseudo-codes are briefly described in Algorithm 1.

The formula for calculating the ith particle fitness is defined as

$$\begin{array}{}\text{(11)}& f_i={\displaystyle \frac{\mathrm{1}}{S}}\sum _{s=\mathrm{1}}^S\sum _{j=\mathrm{1}}^m{\left(Y_{sj}-{\widehat{Y}}_{sj}\right)}^{\mathrm{2}},\end{array}$$

where S is the number of training set samples, m is the output neurons number, Y_sj is the jth true output of the sth sample, and ${\widehat{Y}}_{sj}$ is the corresponding predicted output.

3 Algorithm testing

In order to investigate the COPSO-BP performance and reliability, Rosenbrock and Bohachevsky testing functions were adopted, which are typical non-convex functions and mainly evaluate the performance of unconstrained algorithms. However, due to the random nature of the function, it is not easy to solve and has a global minimum function value of zero.

3.1 Rosenbrock function

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{f}_{\mathrm{1}}\left(x\right)& =\mathrm{100}\times {\left(x_{\mathrm{1}}^{\mathrm{2}}-x_{\mathrm{2}}\right)}^{\mathrm{2}}+{\left(\mathrm{1}-x_{\mathrm{1}}\right)}^{\mathrm{2}},x_i\\ & \in \left[-\mathrm{10},\mathrm{10}\right],i=\mathrm{1},\mathrm{2}\end{array}\end{array}$$

3.2 Bohachevsky function

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{f}_{\mathrm{2}}\left(x\right)& =x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{3}}-x_{\mathrm{1}}{x}_{\mathrm{2}}{x}_{\mathrm{3}}+x_{\mathrm{3}}-\sin \left(x_{\mathrm{2}}^{\mathrm{2}}\right)\\ & -\cos \left(x_{\mathrm{1}}{x}_{\mathrm{3}}^{\mathrm{2}}\right),x_i\in \left[-\mathrm{2}\mathit{\pi },\mathrm{2}\mathit{\pi }\right],i=\mathrm{1},\mathrm{2},\mathrm{3}\end{array}\end{array}$$

Using the standard PSO-BP (SPSO-BP) with linear decreasing inertia weight in Eq. (3), the COPSO-BP tests were carried out. The three-layer BP with an n-s-1 structure is constructed with different hidden nodes. The PSO parameters are the population size M=60, learning factors $c_{\mathrm{1}}=c_{\mathrm{2}}=\mathrm{2.0}$, maximum iteration T_max=30, inertia weight ω_s=0.9, ω_e=0.4, x₀=0.234 and μ=4 for chaotic parameters, and the search dimension $D=n\times s+s\times \mathrm{1}+s+\mathrm{1}$, which includes all the neuron weights and thresholds. For BP network, 150 training samples and 50 testing samples were randomly produced within the variable range. The training error is defined as

$$\begin{array}{}\text{(14)}& E={\displaystyle \frac{\mathrm{1}}{S}}\sum _s^S{\left(T_s-O_s\right)}^{\mathrm{2}},\end{array}$$

where S is the training samples number and T_s and O_s are the expected and predicted outputs for training sample s, respectively. The network structures with minimum training errors for the Rosenbrock and Bohachevsky functions are 2-7-1 and 3-6-1, respectively. The simulation is performed 20 times for each testing function with the SPSO-BP and COPSO-BP algorithms. The numerical result was shown in Table 1. One of the evolutionary training error curves (randomly selected once in 20 times) was shown in Fig. 3, and the fitting curves of the COPSO-BP algorithm were shown in Fig. 4.

Table 1 Comparison of the SPSO-BP and COPSO-BP algorithms for testing functions.

Figure 4 Fitting curves of the COPSO-BP algorithm.

Figure 5 Influence of hidden-layer nodes on R² for different geoelectric models.

Figure 6 Distribution of resistivity ρ₁ and thickness h₁ in training samples.

It can be seen in Table 1 that although both the SPSO-BP and COPSO-BP algorithms can acquire relatively high accuracy for testing functions, the COPSO-BP algorithm is slightly better than the SPSO-BP algorithm. However, the COPSO-BP algorithm has a better convergence rate and optimization efficiency in the early stage in Fig. 3. Therefore, the SPSO-BP and COPSO-BP algorithms have a strong learning ability, good stability and generalization ability, which will be suitable for TEM inversion.

4 Layered-model and parameter analysis

4.1 Forward model

According to the derivation of Kaufman and Keller (1983), the frequency response of the central loop source for the layered model takes the following Hankel transform:

$$\begin{array}{}\text{(15)}& H_z(\mathit{\rho },\mathit{\omega })=Ia\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{m^{\mathrm{2}}}{m+m_{\mathrm{1}}/R_{\mathrm{1}}^{\ast }}}{J}_{\mathrm{1}}\left(m\mathit{\rho }\right)\mathrm{d}m,\end{array}$$

where a is the radius of the transmitting coil, I is the excitation current, ρ is the center distance between the transmitting coil and the receiving coil, J₁(mρ) is the 1st-order Bessel function, m is the integral variable, m₁ is equal to (m²-$k_{\mathrm{1}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$, k₁ is the conduction current, σ₁ is the conductivity, k₁ is equal to −iωμσ₁, and $R_{\mathrm{1}}^{\ast }$ is the first-layer apparent resistivity conversion function, which can be obtained by the following recurrence formula:

$$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}{R}_n^{\ast }=\mathrm{1}\\ R_j^{\ast }=\frac{m_j{R}_{j+\mathrm{1}}^{\ast }+m_{j+\mathrm{1}}\mathrm{th}\left(m_j{h}_j\right)}{m_{j+\mathrm{1}}+m_j{R}_{j+\mathrm{1}}^{\ast }\mathrm{th}\left(m_j{h}_j\right)}\end{array}.\right.\end{array}$$

There is no analytical solution for the time-domain response for the layered model; it can only be solved by numerical calculation. The Hankel transform in formula (15) is calculated by an improved digital filtering algorithm with 47 points, the J₁ filter coefficient, and then the time response can be obtained using the Gaver–Stehfest transform as follows:

$$\begin{array}{}\text{(17)}& H_z(\mathit{\rho },t)={\displaystyle \frac{\ln \mathrm{2}}{t}}\sum _{n=\mathrm{1}}^{\mathrm{N}}{K}_n{H}_z(\mathit{\rho },s_n),\end{array}$$

where s_n is equal to $(\ln \mathrm{2}/t)\times n$, K_n is the coefficient and N is determined by the computer bits, generally N =12.

The ramp excitation current of TEM is

$$\begin{array}{}\text{(18)}& I\left(t\right)=\left\{\begin{array}{c}\mathrm{0},t<\mathrm{0}\\ t/T_{\mathrm{1}},\mathrm{0}\le t<T_{\mathrm{1}}\\ \mathrm{1},T_{\mathrm{1}}<t\end{array},\right.\end{array}$$

where T₁ is the turn-off time, and the Laplace transform is

$$\begin{array}{}\text{(19)}& I\left(s\right)={\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}{s}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}{s}^{\mathrm{2}}}}{e}^{-T_{\mathrm{1}}s}={\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}{s}^{\mathrm{2}}}}\left(\mathrm{1}-e^{-T_{\mathrm{1}}s}\right).\end{array}$$

Therefore, for a specific layered model, the apparent resistivity conversion function $R_{\mathrm{1}}^{\ast }$ is firstly calculated by the recurrence in Eq. (16) based on geoelectric structure parameters. And then the frequency response at fixed point H_z(ω) is calculated by a Hankel transform as in Eq. (15). For ramp excitation, the Laplace transform of H_z(s) should multiplied by I(s). Finally, the time response H_z(t) is obtained by a Gaver–Stehfest transform in Eq. (17). So the H_z(t) is obtained by a Gaver–Stehfest transform, a Hankel transform and a recurrence calculation, and it is somewhat computationally consuming.

However, the vertical magnetic field H_z(t) is the actual observed signal in the transient electromagnetic method in engineering applications. It is the inversion input, and the outputs are geoelectric structure parameters. A method which can avoid the complicated forward model calculation is of great importance in algorithm efficiency.

4.2 BP network design and the COPSO algorithm

For BP structure, the output nodes are determined by the number of inversion geoelectrical parameters, the input nodes are determined by the samples number of H_z(t), and the hidden nodes vary according to approximation performance. As a three-layer or five-layer geoelectric model, its geoelectrical parameters are five (three resistivity and two thickness parameters) or nine (five resistivity and four thickness parameters); the output nodes are five or nine, correspondingly. The characteristic samplings of H_z(t) are chosen as 10 or 20, which are determined by the model's complexity, with more layers meaning more sampling points needed. The 10 samplings were selected in this paper, hence with 10 input nodes. The hidden-layer neuron number is directly related to the weight and threshold parameter amount and greatly affects the BP performance. An appropriate hidden-node number is necessary, and a determination coefficient R² is defined for evaluation as

$$\begin{array}{}\text{(20)}& R^{\mathrm{2}}={\displaystyle \frac{{\left(n\sum _{i=\mathrm{1}}^n{Y}_i{\widehat{Y}}_i-\sum _{i=\mathrm{1}}^n{Y}_i\sum _{i=\mathrm{1}}^n{\widehat{Y}}_i\right)}^{\mathrm{2}}}{\left(n\sum _{i=\mathrm{1}}^n{\widehat{Y}}_i^{\mathrm{2}}-{\left(\sum _{i=\mathrm{1}}^n{\widehat{Y}}_i\right)}^{\mathrm{2}}\right)\left(n\sum _{i=\mathrm{1}}^n{Y}_i^{\mathrm{2}}-{\left(\sum _{i=\mathrm{1}}^n{Y}_i\right)}^{\mathrm{2}}\right)}},\end{array}$$

where Y_i is the true value, ${\widehat{Y}}_i$ is the predicted value for the ith training data and n is the training data number. A larger determination coefficient means better approximation performance. The simulations on hidden-nodes effects were carried out for three-layer and five-layer geoelectric models. The BP structure is 10-s-5 and 10-s-9, and its transfer, training and learning functions are the “log sigmoidal”, “Levenberg–Marquardt” and “gradient descent momentum”, respectively. The average, minimum and maximum values of R² were obtained after running 20 times for each simulation. The R² curves were shown in Fig. 5.

It can be seen that the optimal neural-network structures were 10-2-5 and 10-5-9 for three- and five-layer models based on the maximum R² values. Then, the PSO-BP algorithms with different inertia weight were implemented and compared for the three-layer model. The BP structure was chosen as 10-2-5, and four types of inertia weight as in Eqs. (3)–(7) in the PSO were compared in Table 2.

Table 2 Comparison of different inertia weights in PSO algorithms (ω_s=0.9 and ω_e=0.4).

The simulation was implemented on a Core (TM) i5-7500 processor with 8 GB of memory. It is obviously found in Table 2 that the COPSO algorithm has a much faster convergence rate and a lower iteration number and is time consuming.

4.3 Layered-model inversion

A three-layer and five-layer geoelectric models were investigated, and the PSO parameter values are the same as those of the “Algorithm testing” section in this paper. In order to simulate actual TEM applications, the ramp turn-off is taken into account. Considering the probability distribution characteristic of the above algorithms, the average of 20 simulation results was chosen. The BP, SPSO-BP and COPSO-BP algorithms and a nonlinear programming genetic algorithm (NPGA) (Li et al., 2017) were compared.

4.3.1 Three-layer H-type model

The central loop TEM parameters were set as follows: the transmitting coil radius was set to a=100 m; the ramp emission current was set to 100 A; and the turn-off time was set to 1 µs. In the geoelectric model, the resistivity is ρ₁=100, ρ₂=10 and ρ₃=100 Ω m and the thickness is h₁=100 and h₂=200 m.

The BP training samples, which are a series of H_z(t) for different geoelectrical parameters, were generated by the TEM forward model. The resistivity ranges were ρ₁∈(50, 150), ρ₂∈(5, 15) and ρ₃∈(50, 150) and the thickness ranges were h₁∈(50, 150) and h₂∈(100, 300); 1000 random groups were chosen. The resistivity and thickness distributions of ρ₁ and h₁ were shown in Fig. 6. The relative error is defined as

$$\begin{array}{}\text{(21)}& {\mathrm{Err}}_{\mathit{\_}\mathrm{rel}}=\left|\left.{\displaystyle \frac{T_{\mathit{\_}\mathrm{cal}}^{\ast }-O_{\mathit{\_}\mathrm{ref}}^{\ast }}{O_{\mathit{\_}\mathrm{ref}}^{\ast }}}\right|,\right.\end{array}$$

where $T_{\mathit{\_}\mathrm{cal}}^{\ast }$ and $O_{\mathit{\_}\mathrm{ref}}^{\ast }$ are the calculated and reference values for the geoelectric models.

The inversion results were shown in Table 3 and Figs. 7–8. The BP type algorithms were superior to the NPGA inversion in Table 3. Moreover, the inversion accuracy, convergence rate and optimization ability of the COPSO-BP algorithm were better than others.

Table 3 Inversion comparison of three-layer H-type geoelectric models.

Figure 7 Fitness curves of SPSO-BP and COPSO-BP.

Figure 8 Comparison of mean square error curves.

Additional results showed that the solution range of ρ₁ and h₁ in the 20 simulations for the above algorithms were ρ₁∈(97.980, 103.102) and h₁∈(96.962, 102.480) for BP, ρ₁∈(98.954, 101.137) and h₁∈(96.955, 101.829) for SPSO-BP, ρ₁∈(99.382, 100.989) and h₁∈(97.877, 101.044) for COPSO-BP, respectively. Therefore, the COPSO-BP can acquire a higher accuracy and is more stable.

4.3.2 Five-layer KHK-type model

A five-layer KHK-type geoelectric model was adopted, and its resistivities were ρ₁=100, ρ₂=300, ρ₃=50, ρ₄=200 and ρ₅=30 Ω m, and its thickness were h₁=100, h₂=200, h₃=300 and h₄=500 m.

The training samples with parameter ranges were ρ₁∈(50, 150), ρ₂∈(150,450), ρ₃∈(25, 75), ρ₄∈(100, 300) and ρ₅∈(15, 45) for resistivity and h₁∈(50, 150), h₂∈(100, 300), h₃∈(150, 450) and h₄∈(250, 750) for thickness. The 1000 group training samples were generated within the above ranges. The inversion results were shown in Table 4 and Figs. 9–10. It can be seen that the COPSO-BP algorithm has better global optimization performance.

Table 4 Inversion comparison of five-layer KHK-type geoelectric models.

Figure 9 Fitness curves of SPSO-BP and COPSO-BP.

Figure 10 Comparison of mean square error curves.

4.3.3 Inversion comparison

Three kinds of BP methods, the traditional BP, SPSO-BP and COPSO-BP algorithms, were compared in Table 5. Hence, the training times of COPSO-BP were obviously less than those of SPSO-BP and were almost equal to BP; it can obtain better precision especially for its global optimization performance.

Table 5 Simulation comparison of different algorithms.

Figure 11 Inversion comparison of different geoelectric models.

Figure 12 Forward data of Hz and data with 5 % noise.

The inversion of COPSO-BP and NPGA were compared in Fig. 11. The fitting ability of COPSO-BP was much better than NPGA.

Table 6 Comparison of inversion results of three-layer H-type (with noise) models.

Figure 13 1-D inversion of forward results for (a) BP and (b) COPSO-BP.

Figure 14 Inversion results for BP (a) and COPSO-BP (b).

4.3.4 Robust performance analysis

In order to verify the algorithm robustness, 5 % (26 dB) and 10 % (20 dB) Gaussian random noise was added in TEM data for the three-layer geoelectric model. Three kinds of inversions were implemented respectively. The results and a comparison were shown in Table 6. The H_z(t) and data with 5 % noise were shown in Fig. 12.

As can be seen in Table 2, after applying 5 % and 10 % Gaussian noise, the COPSO-BP inversion has a higher robust ability. The accuracy was obviously improved based on the total relative-error data.

4.4 Field example

In order to test the effectiveness of the method, a transient electromagnetic vertical magnetic field (Hz) with 10 measuring points at 380 to 1280 m of the no. 1 line from a mining area in Anhui Province were selected. After the data processing, the inversion was performed using the three-layer neural-network model in the previous section, and the results of the BP and COPSO-BP inversion were compared. Figure 13 shows the comparison between the surveyed data and the inversion data at 380 m of the no. 2 line in the mining area. Figure 14 displays the pseudo-sections of the 10 sets of inversion data combined with the geological data interpolation smoothing. It can be seen from Fig. 14 that the first layer has low resistivity (100–200 Ω m), which is inferred to be the second-layer (T2g22) gray dolomite of the Middle Triassic old Malague section, with a thickness of about 200 m. The second layer has the second-highest resistivity (300–400 Ω m), which is surmised to be the first-layer (T2g21) dolomite of the Middle Triassic old Malague section, with a thickness of about 400 m. The third layer has high resistivity (600–800 Ω m), which is speculated to be the sixth-layer (T2g16) limestone dolomite of the Middle Triassic old group. The results are basically consistent with the geological conditions of the mining area, indicating the feasibility and effectiveness of the neural-network method. And the results of COPSO-BP inversion are better than those of the BP, for which the inversion position is more accurate, the shape and spacing are clearer, and the resistivity of each layer is more consistent with the those of the actual geological model.

5 Discussion

The inversion is performed for three-layer (H-type) and five-layer (KHK-type) geoelectric models in this paper. The results show that the BP neural network is better than the NPGA algorithm because the BP method does not need to use the forward algorithm repeatedly, and its calculation time is short, which is different from the nonlinear heuristic method based on a global space search solution.

The main advantage of the BP is that it can interpret the transient electromagnetic sounding results quickly after training the network. Furthermore, the BP algorithm could automatically obtain the “reasonable rules” between input and output data by learning, and it can adaptively store the learning content in the network weight, for which the BP neural network has a high self-learning and self-adaptation ability. In addition, the superior simulation results of the test function indicate that the BP algorithm can approximate any nonlinear continuous function with arbitrary precision, which means it has a strong nonlinear mapping ability; the inversion results of the layered geoelectric model with uncorrelated noise data prove that the BP algorithm has strong robustness, which means it has the ability to apply learning results to new knowledge. However, the BP neural-network weight is gradually adjusted by the direction of local improvement, which causes the algorithm to fall into a local extremum, and the weight converges to a local minimum that leads to the network training failure. Moreover, the BP is very sensitive to the initial network weight, and the initialization network with different weight values tends to converge at different local minimums, so it obtains different results each time. In addition, the BP algorithm is essentially a gradient descent method, which leads to a slow convergence rate.

From the results of the layered-model and parametric analysis part, it can be seen that the single BP algorithm has a higher error value than SPSO-BP because the BP method is sensitive to initial weight and easy to fall into local minimum values; thus a heuristic global search particle swarm optimization algorithm with a simple structure, rapid convergence and high precision is applied to optimize the weight and threshold of the BP neural network, which improves the global optimization performance of the algorithm. Furthermore, the PSO algorithm adjusts the inertia weight adaptively based on the chaotic-oscillation curve that is similar to the annealing process in the simulated annealing algorithm (SA), which jumps out the local extremum faster in the early stage and accelerates the convergence and reduces the training times. Therefore, compared with the SPSO-BP and BP algorithms, the inversion results of COPSO-BP are closer to the theoretical data with smaller error fluctuations, stronger anti-noise controls, a better generalization of performance and higher stability, which it is effective in solving geophysical inverse problems.

From the simulation experiment, it is not clear how the weight organization affects the BP neural-network weight learning process. It is necessary to conduct a more systematic study on this problem to improve our understanding of how the BP neural network handles training data.

6 Conclusions

The nonlinear COPSO-BP method was proposed for TEM inversion. The BP's initial weight and threshold parameters were trained by the COPSO algorithm, which makes it easy to not fall into a local optimum. The chaotic-oscillation inertia weight for PSO was proposed so as to improve the PSO's global optimization ability and fast convergence in the early stage. The layered geoelectric model inversion showed that the COPSO-BP method is more accurate and stable and requires relatively less training time.