Nonlinear Inversion for Complex Resistivity Method Based on QPSO-BP Algorithm

Weixin Zhang^1*, Jinsuo Liu², Le Yu², Biao Jin²

¹School of Resources and Geosciences, China University of Mining and Technology, Xuzhou, China.
²Anhui Construction Engineering Testing and Research Institute, Hefei, China.

DOI: 10.4236/ojg.2021.1110026                

Abstract

The significant advantage of the complex resistivity method is to reflect the abnormal body through multi-parameters, but its inversion parameters are more than the resistivity tomography method. Therefore, how to effectively invert these spectral parameters has become the focused area of the complex resistivity inversion. An optimized BP neural network (BPNN) approach based on Quantum Particle Swarm Optimization (QPSO) algorithm was presented, which was able to improve global search ability for complex resistivity multi-parameter nonlinear inversion. In the proposed method, the nonlinear weight adjustment strategy and mutation operator were used to enhance the optimization ability of QPSO algorithm. Implementation of proposed QPSO-BPNN was given, the network had 56 hidden neurons in two hidden layers (the first hidden layer has 46 neurons and the second hidden layer has 10 neurons) and it was trained on 48 datasets and tested on another 5 synthetic datasets. The training and test results show that BP neural network optimized by the QPSO algorithm performs better than the BP neural network without initial optimization on the inversion training and test models, and the mean square error distribution is better. At the same time, a double polarized anomalous bodies model was also used to verify the feasibility and effectiveness of the proposed method, the inversion results show that the QPSO-BP algorithm inversion clearly characterizes the anomalous boundaries and is closer to the values of the parameters.

Keywords

Complex Resistivity, Finite Element Method, Nonlinear Inversion, QPSO-BP Algorithm, 2.5D Numerical Simulation

1. Introduction

CR (Complex resistivity) method, also known as SIP (Spectrum Induced Polarization) method, is one of the geophysical prospecting methods that is based on the different induced polarization characteristics of rocks and ores in the frequency domain. Multi-parameter is the main feature of complex resistivity inversion [1]. The essence of inversion theory is to study the theory and method of mapping the observed data in Geophysics to the corresponding geophysical model [2]. For the complex resistivity inversion problem, its essence is a nonlinear problem. The nonlinear problem is easy to fall into the local extremum by using the traditional linear or quasi-linear method, so it needs to be solved by the nonlinear method.

In the inversion of electrical resistivity tomography, the neural network has become one of the most widely used complete nonlinear inversion methods because of its strong nonlinear mapping ability and easy construction.

Ahmad Neyamadpour et al. used the artificial neural network to study the inversion imaging of 2D DC resistivity data and created a Matlab application. The interpreted result shows that the trained network was able to invert 2D electrical resistivity imaging data obtained by a Wenner-Schlumberger configuration rapidly and accurately [3]. Dai Qian-Wei et al. Completed the Nonlinear inversion for electrical resistivity tomography based on chaotic DE-BP algorithm, and the results show that the proposed method has better performance in stability and accuracy and higher imaging quality than least-square inversion [4]. Trong Long Ho used a back-propagation neural network to research 3D inversion problem of borehole-to-surface electrical data. The result shows that the most successful learning algorithm in this network is the Resilient Propagation (RPROP) [5]. The above researchers have made a lot of contributions in resistivity nonlinear inversion, but they have not introduced the nonlinear method into the complex resistivity inversion.

However, for the complex resistivity method, which requires multi-parameter inversion, the application of a neural network is rare. The inversion methods are mostly linear or quasi-linear methods. Loke et al. introduced the regularization method and adopted the smooth constrained least square method to carry out the inversion of 2.5D complex resistivity method with ignoring the electromagnetic coupling effect, and obtained the parameters of Cole-Cole mode, but the inversion of the boundary is not as expected [6]. Son et al. proposed a new complex resistivity inversion algorithm based on DC resistivity, which was calculated in the complex domain. In order to directly invert the complex resistivity of multi-frequency points, the algorithm decomposes the Cole-Cole model parameters of multi-frequency complex resistivity [7]. Blaschek et al. introduced the minimum gradient regularization algorithm to study the inversion problem of complex resistivity [8]. Routh et al. studied the complex resistivity inversion of the smooth model under the condition of considering electromagnetic coupling [9]. Hönig et al. adopted a forward algorithm similar to DC resistivity, then introduced the zero-frequency smoothing constraint of complex resistivity, and studied the regularization inversion of complex resistivity [10]. Zhang Zhi-Yong et al. used MPI parallel algorithm to study 2D inversion problem of spectral induced polarization data. The tests result suggests that the proposed parallel algorithm is robust and efficient [11].

Based on the above research experience, this paper introduces the optimized neural network to solve the multi-parameters inversion problem of the complex resistivity method and constructs a nonlinear mapping network between apparent complex resistivity and inversion parameters. We also propose a complex resistivity nonlinear inversion algorithm based on QPSO-BP algorithm to try to solve the multi-parameter nonlinear inversion problem.

2. Forward Modeling Theory of 2.5D Complex Resistivity Method

2.1. Forward Modeling Theory

Under the condition of ignoring the electromagnetic coupling effect, the complex resistivity forward calculation is similar to the DC method. The Cole-Cole model is introduced to characterize the induced polarization effect. After Fourier transform, 2D complex potential boundary value problem in wavenumber domain is [12]:

$\{\begin{array}{l}\nabla \cdot \left(\sigma \nabla U\right)-k^2\sigma U=-I\delta \left(A\right)\in \Omega \\ \frac{\partial U}{\partial n}=0\in {\Gamma }_s\\ \frac{\partial U}{\partial n}+k\frac{K_1\left(kr\right)}{K_0\left(kr\right)}cos\left(r,n\right)U=0\in {\Gamma }_{\infty }\end{array}$ (1)

In Equation (1), the area Ω is a three-dimensional study area, Boundary ${\Gamma }_s,{\Gamma }_{\infty }$ is the boundary of a two-dimensional region. Among them, ${\Gamma }_s$ is the surface boundary of area Ω, and ${\Gamma }_{\infty }$ is the underground boundary of area Ω, As shown in Figure 1, ${\sigma }_2$ is the complex conductivity of the polarized medium, r is the distance from the power source to the boundary, n is the outer normal direction of the infinity boundary, k is the wavenumber, $K_0$ is the zero-order Bessel function, and $K_1$ is the second kind of first-order Bessel function. The functional of Equation (1) is:

Figure 1. Mode schematic of two-dimensional structure.

$\begin{array}{c}F\left(U\right)={\displaystyle {\int }_{\Omega }\left[\frac{\sigma }{2}{\left(\nabla U\right)}^2+\frac{1}{2}{k}^2\sigma U^2-I\delta \left(A\right)U\right]\text{d}\Omega }\\ +\frac{1}{2}{\displaystyle {\int }_{{{\displaystyle \Gamma }}_{\infty }}\sigma k\frac{K_1\left(kr\right)}{K_0\left(kr\right)}\cos \left(r,n\right)U^2\text{d}\Gamma }\end{array}$ (2)

Discrete the calculation area and construct an interpolation function. The discrete form of Equation (2) is:

$F\left(U\right)=\frac{1}{2}{U}^{{\rm T}}KU-U^{{\rm T}}P$ (3)

where P is the column vector related to the excitation point power. Solving the variation of the Equation (3) and making it equal to zero, then we can get the system of linear equations:

$KU=P$ (4)

The wavenumber domain potential U of each point is obtained by solving the linear Equation (4), and the potential is obtained by inverse Fourier transform [13].

$u\left(r\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}U\left(r,k_j\right)g_j}$ (5)

where r is the position of the measuring point, $k_j$ is the wavenumber and $g_j$ is the weighting coefficient, the values of both are shown in Table 1 [14].

The complex resistivity produced by the induced polarization effect of rock and ore is expressed by Cole-Cole model [15]:

$\rho \left(i\omega \right)={\rho }_0\left\{1-m\left[1-\frac{1}{1+{\left(i\omega \tau \right)}^c}\right]\right\}$ (6)

where, ${\rho }_0$ is resistivity (0 Hz); m is charge ability; c is relaxation-constant; τ is a time constant. The complex resistivity forward calculation can be carried out by replacing the conductivity in Equation (1) and Equation (2) with the complex conductivity.

2.2. Accuracy Verification of Forwarding Modeling Algorithm

In order to verify the accuracy of the CR forward algorithm, the two-layer horizontal complex resistivity model is used to verify the calculation accuracy.

For the two-layer horizontal formation, by solving the simplified Laplace equation of the cylindrical coordinate system and power series expansion, the calculation expression of the complex potential $\tilde{U}$ at each point on the surface can be obtained as:

Table 1. The inverse Fourier transform coefficient.

$\tilde{U}=\frac{I{\rho }_1\left(i\omega \right)}{2\pi }\left[\frac{1}{r}+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{K_{12}^n}{\sqrt{r^2+{\left(2n{h}_1\right)}^2}}}\right]$ (7)

In Equation (7), $K_{12}=\left[{\rho }_2\left(i\omega \right)-{\rho }_1\left(i\omega \right)\right]/\left[{\rho }_2\left(i\omega \right)+{\rho }_1\left(i\omega \right)\right]$, ${\rho }_1\left(i\omega \right)$, ${\rho }_2\left(i\omega \right)$ is the complex resistivity calculated by Equation (6) for the first and second layers, r is the distance between the measuring point and the origin, $h_1$ is the thickness of the first stratum. The calculation accuracy of Equation (7) depends on the size of the number of terms n. The larger the n is, the higher accuracy of the calculation is. This calculation takes n = 200. The parameters of the two-layer polarization model are shown in Table 2.

The numerical solution and analytical solution of 2.5D complex potential of the two-layer horizontal layer model are shown in Figure 2(a) and Figure 2(b). Both the amplitude and phase curves of complex potential are consistent with the analytical solution, which shows the accuracy of the finite element algorithm. The finite element algorithm provides a reliable basis for constructing BP neural network training database.

(a)(b)

Figure 2. Comparison of analytical solution of complex potential and numerical solution of finite element. (a) Excitation frequency is 1.25 Hz; (b) Excitation frequency is 125 Hz.

Table 2. Parameter table of two-layer polarization model.

3. Quantum Particle Swarm Algorithm for Training BP Neural Network

3.1. Quantum Particle Swarm Optimization Algorithm

Quantum particle swarm optimization (QPSO) absorbs the characteristics and advantages of PSO algorithm and quantum computing, combines the two methods, and presents the solution of the problem to be optimized in the form of qubit probability amplitude [16].

Initialization: The qubit phase θ in the range of $\left[0,2\pi \right]$ is randomly generated by using a random number function, and the qubit can be obtained, which can be expressed as a probability amplitude ${\left[\sin \theta ,\cos \theta \right]}^{{\rm T}}$. After that, the solution space is transformed by combining the upper and lower limits of each variable set by the algorithm, and the qubit is transformed into the solution space to calculate the appropriate value.

Update: In the QPSO algorithm, with the help of quantum rotation gates operation, the individual and group optimal positions (set as cosine positions) found by the particle $p_i$ in the current solution space can be expressed as:

$\begin{array}{l}{p}_{il}=\left[\cos \left({\theta }_{il1}\right),\cos \left({\theta }_{il2}\right),\cdots ,\cos \left({\theta }_{ilD}\right)\right]\\ p_g=\left[\cos \left({\theta }_{g1}\right),\cos \left({\theta }_{g2}\right),\cdots ,\cos \left({\theta }_{gD}\right)\right]\end{array}$ (8)

The updated rules of particle state are as follows:

The incremental update of qubit angle on particles is:

$\Delta {\theta }_{ij}\left(t+1\right)=\omega \Delta {\theta }_{ij}\left(t\right)+c_1{r}_1\left(\Delta {\theta }_l\right)+c_2{r}_2\left(\Delta {\theta }_g\right)$ (9)

where

$\Delta {\theta }_l=\{\begin{array}{l}2\pi +{\theta }_{ilj}-{\theta }_{ij}\left({\theta }_{ilj}-{\theta }_{ij}<-\pi \right)\\ {\theta }_{ilj}-{\theta }_{ij}\left(\pi \le {\theta }_{ilj}-{\theta }_{ij}\le \pi \right)\\ {\theta }_{ilj}-{\theta }_{ij}-2\pi \left({\theta }_{ilj}-{\theta }_{ij}>\pi \right)\end{array}$, $\Delta {\theta }_g=\{\begin{array}{l}2\pi +{\theta }_{gj}-{\theta }_{ij}\left({\theta }_{gj}-{\theta }_{ij}<-\pi \right)\\ {\theta }_{gj}-{\theta }_{ij}\left(\pi \le {\theta }_{gj}-{\theta }_{ij}\le \pi \right)\\ {\theta }_{gj}-{\theta }_{ij}-2\pi \left({\theta }_{gj}-{\theta }_{ij}>\pi \right)\end{array}$

And ω is a random number with inertia weight; $c_1,c_2$ is a self-factor and a global factor respectively; $r_1,r_2$ is a random number of (0, 1).

The probability amplitude of qubit on particles is updated as:

$\begin{array}{c}\left[\begin{array}{c}\cos \left({\theta }_{ij}\left(t+1\right)\right)\\ \sin \left({\theta }_{ij}\left(t+1\right)\right)\end{array}\right]=\left[\begin{array}{cc}\cos \left(\Delta {\theta }_{ij}\left(t+1\right)\right)& -\sin \left(\Delta {\theta }_{ij}\left(t+1\right)\right)\\ \sin \left(\Delta {\theta }_{ij}\left(t+1\right)\right)& \cos \left(\Delta {\theta }_{ij}\left(t+1\right)\right)\end{array}\right]\left[\begin{array}{c}\cos \left({\theta }_{ij}\left(t\right)\right)\\ \sin \left({\theta }_{ij}\left(t\right)\right)\end{array}\right]\\ =\left[\begin{array}{c}\cos \left({\theta }_{ij}\left(t\right)+\Delta {\theta }_{ij}\left(t+1\right)\right)\\ \sin \left({\theta }_{ij}\left(t\right)+\Delta {\theta }_{ij}\left(t+1\right)\right)\end{array}\right]\end{array}$ (10)

Among them $i=1,2,\cdots ,n;j=1,2,\cdots ,d$, n and d are the number of population and the dimension of unknown variables, respectively.

Mutation: In the QPSO algorithm, quantum non-gate is used to mutate particles. First, the mutation probability $K_m$ is set, and each particle is given a random number $ran{d}_i$ between (0, 1). If $ran{d}_i<K_m$, then [n/2] qubits are randomly selected, and the mutation operation is carried out by using Equation (11), and the memorized optimal position and rotation angle vector remains unchanged.

$\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}\cos \left({\theta }_{ij}\right)\\ \sin \left({\theta }_{ij}\right)\end{array}\right]=\left[\begin{array}{c}\sin \left({\theta }_{ij}\right)\\ \cos \left({\theta }_{ij}\right)\end{array}\right]=\left[\begin{array}{c}\cos \left({\theta }_{ij}+\frac{\pi }{2}\right)\\ \sin \left({\theta }_{ij}+\frac{\pi }{2}\right)\end{array}\right]$ (11)

where $i=1,2,\cdots ,n;j=1,2,\cdots ,d$.

Weight adjustment: In particle swarm optimization, inertia weight ω is an important parameter. If the value of ω is increased, the global search ability will be enhanced. If the value of ω is decreased, the local search ability will be enhanced. Equation (12) is used to realize the nonlinear dynamic adjustment of inertia weight:

$\varpi =\{\begin{array}{l}{\varpi }_{\min }-\frac{\left({\varpi }_{\max }-{\varpi }_{\min }\right)\ast \left(f-f_{\min }\right)}{f_{avg}-f_{\min }},f\le f_{\text{avg}}\\ {\varpi }_{\max },f>f_{avg}\end{array}$ (12)

where f represents the real-time objective function value of particles; $f_{\min }$ and $f_{\text{avg}}$ are the minimum and average moderate values of all particles.

3.2. Back Propagation Neural Network Construction

BP neural network modeling method: take the horizontal position, vertical position, amplitude and phase information of measuring apparent complex resistivity as the input information of neural network, and take the parameter information of model as the output information of the neural network.

In order to adapt to the multi-parameter characteristics of complex resistivity inversion, the BP neural network model obtained through numerical test is shown in Figure 3. The overall neural network has 12 input nodes in the input layer, 56 hidden neurons in two hidden layers (the first hidden layer has 46 neurons and the second hidden layer has 10 neurons), 4 output nodes in the output layer, and the overall topology is 12-46-10-4.

3.3. Training Process of Back Propagation Neural Networks

QPSO-BP algorithm mainly uses the global optimization ability of QPSO to optimize the initial weight and threshold of BP neural network. In order to obtain the best overall effect, BP neural network uses the optimized initial weight and threshold to continue training.

The inversion algorithm flow is as follows:

Figure 3. Schematic diagram of neural network structure.

Step 1: Determine the topology of BP neural network and the number of parameters to be optimized, and set the particle population dimension. Each individual qubit phase θ is a D-dimensional parameter vector corresponding to the weights and thresholds of BP neural network. Initialize the particle population and set the basic parameters of the algorithm, such as $\varpi$ adjustment range $\left[{\omega }_{\min },{\omega }_{\max }\right]$, self-factor $c_1$, global factor $c_2$, mutation probability $K_m$.

Step 2: According to the upper and lower limits of each particle, the solution space is transformed. Then, according to Equation (13), calculate the particle fitness value, select and record the individual optimal position and group optimal position found in the solution space.

$\begin{array}{l}{F}_k=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{l}{\sum }}{\left(e_{kj}-y_{kj}\right)}^2}\\ f_i=\frac{1}{cn}{\displaystyle \underset{k=1}{\overset{cn}{\sum }}{\left(F_k\right)}^2}\end{array}$ (13)

where, cn is the total number of inversion parameters, k is the number of inversion parameters, $k=1,2,\cdots ,cn$, n is the total number of training samples, j is the number of training samples, $j=1,2,\cdots ,n$. $e_{kj}$ and $y_{kj}$ are the expected output and actual neural network output of the j-th training sample of the k-th parameter, respectively. $F_k$ is the mean square error of the k-th parameter. $f_i$ is the fitness value of the i-th particle.

Step 3: particle state update, nonlinear dynamic adjustment of inertia weight, and qubit mutation calculation operations.

Step 4: Compare the appropriate value of each particle in order. Then the optimal particle position is recorded and updated as an individual optimal position. At the same time, the optimal position of each generation is sorted forward in the current evolution generation, and the optimal position of the individual is recorded and updated as the global optimal position.

Step 5: Identify whether the iteration meets the end condition, if it is satisfied, exit the algorithm and output the result, otherwise return to step 2 to continue the iterative calculation.

Step 6: The global optimal solution is used as the initial weight and threshold of BP neural network. Then we continue to train the neural network by using the BP elastic algorithm (RPROP). When the number of algorithm iterations reaches the maximum or the output error satisfies the requirements, output and save the network. At the last, we input the data volume which needs to be inverted into the network to obtain the QPSO-BP algorithm inversion result.

3.4. Performance Analysis of Inversion Algorithm

In this paper, we construct an apparent complex resistivity database with 48 single-polarization anomaly models. The forward simulation of the database model was calculated out by the finite element method described in Section 2. All the data are used for QPSO-BP neural network training. There are two sizes of abnormal models in the database; each of them has six different spatial locations. We also set a variety of contrast characteristics for these abnormal bodies, such as low resistance with high polarization, high resistance with low polarization. Five independent models are used to test the generalization ability of the algorithm. The test data is not involved in neural network training. Part of the sample model used for training is shown in Figure 4.

At the same time, QPSO algorithm and basic PSO algorithm are used to optimize the initial value and threshold of BP neural network. Set the particle population size of QPSO algorithm to 80, set the evolution generation to 100, ${\omega }_{\max }=0.6$, ${\omega }_{\min }=0.4$, $c_1=c_2=1.2$, $K_m=0.35$. Set the particle population

Figure 4. Schematic diagram of partial training sample model.

size of basic PSO algorithm to 80, set the evolution generation to 100, $\omega =0.5$, $c_1=1.5$, $c_2=2.5$. The fitness decline curves of the two algorithms are shown in Figure 5(a).

QPSO algorithm is easier to jump out of the local extreme value in the early calculation and converge to the global extreme value in the later calculation than the basic PSO algorithm with fixed weight because it has more particle number (sine solution and cosine solution) and nonlinear dynamic adjustment weight.

Using the global optimal solution as the initial weight and threshold of the BP neural network, and using the BP elastic algorithm (RPROP) to continue training the neural network, finally, the convergence curve of algorithm training error is shown in Figure 5(b). The optimized BP neural network can converge to the global optimal value more quickly, and its final training MSE is lower than that of the non-optimized network. At the same time, Figure 6 shows the difference of MSE distribution histogram between training and test samples under the two algorithms.

The BP neural network optimized by QPSO algorithm performs better than the BP neural network without initial optimization on the inversion training and test models, and the mean square error distribution is better. This shows that the BP neural network optimized by QPSO has a better training effect.

(a)(b)

Figure 5. Decline curve of moderate value and training error. (a) Fitness decline curve of different algorithms; (b) Training error convergence curve.

(a)(b)

Figure 6. Mean square error distribution Histogram of different algorithms. (a) The QPSO-BP algorithm; (b) The BP elastic algorithm (RPROP).

4. Research on Inversion Imaging of Typical Model

In order to verify the feasibility and accuracy of the QPSO-BP algorithm inversion, we established a model with two complex resistivity anomalies, and input the finite element simulation calculation data of the model as the test data into the algorithm for inversion; Figure 7 shows the model structure. Abnormal body No.1: size is 2 m × 3 m (Length × depth), ${\rho }_0=50\Omega \cdot \text{m}$, $m=0.6$, $c=0.2$, $\tau =10\text{s}$ ; Abnormal body No.2: size is 2 m × 3 m (Length × depth), ${\rho }_0=500\Omega \cdot \text{m}$, $m=0.4$, $c=0.3$, $\tau =20\text{s}$, the horizontal distance between the two polarization anomalies is 10 m. A total of 37 electrodes are set on the surface, the electrode spacing is 1m, and the observation frequency points are 10^-2 Hz, 10^-1 Hz, 1 Hz, 10 Hz, 100 Hz, take the apparent complex resistivity amplitude and phase corresponding to the observed frequency points as the input of the QPSO-BP network and start the inversion calculation.

From the inversion results, as shown in Figure 8, the QPSO-BP inversion algorithm can more accurately describe the spatial shape of the anomaly and the four parameters to be inverted by the complex resistivity method. At the same time, the QPSO-BP inversion algorithm has a clearer description of the boundary of the abnormal body, and the inversion results are more consistent with the model, which reflects the strong nonlinear inversion ability of the algorithm.

In order to analyze the inversion effect of the algorithm more intuitively, we extracted the result data of each parameter in the inversion result data at $z=-3.5\text{m}$ and $8\text{m}\ge x\ge 29\text{m}$, and compared them with the model data. As shown in Figure 9, the QPSO-BP inversion result corresponds well to the model value, and the boundary of the inversion result is clearly described, except for the slight deviation of the inversion result near the anomaly.

Figure 7. Schematic diagram of double polarized anomalous bodies model.

Figure 8. Inversion results of single double polarization anomalous body by using QPSO-BP algorithm. (a) The Resistivity; (b) The Chargeability; (c) The Relaxation-constant; (d) The Cole-Cole Time constant.

(a) (b) (c) (d)

Figure 9. Comparison of inversion results of QPSO-BP algorithm at z = −3.5 m. (a) The Resistivity; (b) The Chargeability; (c) The Relaxation-constant; (d) The Cole-Cole Time constant.

5. Conclusions

1) The BP neural network optimized by the QPSO algorithm performs better than the BP neural network without initial optimization on the inversion training and test models, and the mean square error distribution is better. This shows that the BP neural network optimized by QPSO has a better training effect.

2) The QPSO-BP algorithm can effectively invert the four parameters of complex resistivity, which corresponds well to the model value. At the same time, the inversion results clearly describe the abnormal boundary.

3) BP neural network has high recognition and association ability for trained samples or samples with similar characteristics, but poor recognition ability for models with different feature types. For the actual complex data inversion, the method proposed in this paper has limitations. The classical algorithm should be used as the basis for the construction of a neural network training sample database in the early inversion stage.

4) Despite this contribution, there are many remaining challenges for future work. First of all, due to the complexity of multi-parameter inversion of complex resistivity, the selected training samples should be close to the type of model to be inverted as much as possible, and it is essential to choose the right training models with appropriate characteristics. Secondly, the inversion method proposed in this paper has space limitations. The inversion of anomalous bodies outside the observed apparent resistivity profile needs to be further studied. Finally, using the algorithm in this paper to process the inversion of the actual complex resistivity data would be challenging work.

Acknowledgements

This paper is supported by the National Natural Science Foundation of China (Nos. 41974149, 41804106).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zhang, Z.-Y. and Zhou, F. (2014) 2.5D Modeling of Complex Resistivity Using Finite Element Method. Progress in Geophysics, 29, 2326-2331. (In Chinese)
[2]	Wang, J.Y. (2007) Introduction to Geophysical Inverse Problems. Chinese Journal of Engineering Geophysics, 4, 1-3.
[3]	Neyamadpour, A., Taib, S. and Wan Abdullah, W.A.T. (2009) Using Artificial Neural Networks to Invert 2D DC Resistivity Imaging Data for High Resistivity Contrast Regions: A MATLAB Application. Computers & Geosciences, 35, 2268-2274. https://doi.org/10.1016/j.cageo.2009.04.004
[4]	Dai, Q.-W., Jiang, F.-B. and Dong, L. (2014) Nonlinear Inversion for Electrical Resistivity Tomography Based on Chaotic DE-BP Algorithm. Journal of Central South University, 21, 2018-2025. https://doi.org/10.1007/s11771-014-2151-9
[5]	Ho, T.L. (2009) 3-D Inversion of Borehole-to-Surface Electrical Data Using a Back-Propagation Neural Network. Journal of Applied Geophysics, 68, 489-499. https://doi.org/10.1016/j.jappgeo.2008.06.002
[6]	Loke, M.H., Chambers, J.E. and Ogilvy, R.D. (2006) Inversion of 2D Spectral Induced Polarization Imaging Data. Geophysical Prospecting, 54, 287-301. https://doi.org/10.1111/j.1365-2478.2006.00537.x
[7]	Son, J.S., Kim, J.H. and Yi, M.J. (2007) A New Algorithm for SIP Parameter Estimation from Multi-Frequency IP Data: Preliminary Results. Exploration Geophysics, 38, 60-68. https://doi.org/10.1071/EG07009
[8]	Blaschek, R., Hordt, A. and Kemna, A. (2008) A New Sensitivity-Controlled Focusing Regularization Scheme for the Inversion of Induced Polarization Data Based on the Minimum Gradient Support. Geophysics, 73, F45-F54. https://doi.org/10.1190/1.2824820
[9]	Routh, P.S. and Oldenburg, D.W. (2001) Electromagnetic Coupling in Frequencydomain Induced Polarization Data: A Method for Removal. Geophysical Journal International, 145, 59-76. https://doi.org/10.1111/j.1365-246X.2001.00384.x
[10]	Honig, M. and Tezkan, B. (2007) 1D and 2D Cole-Cole-Inversion of Time-Domain Induced-Polarization Data. Geophysical Prospecting, 55, 117-133. https://doi.org/10.1111/j.1365-2478.2006.00570.x
[11]	Zhang, Z.Y., Tan, H.D., Wang, K.P., Lin, C.H. and Xie, M.B. (2016) Two-Dimensional Inversion of Spectral Induced Polarization Data Using MPI Parallel Algorithm in Data Space. Applied Geophysics, 13, 13-24. https://doi.org/10.1007/s11770-016-0530-8
[12]	Coggon, J.H. (1971) Electromagnetic and Electrical Modeling by the Finite Element Method. Geophysics, 36, 132-155. https://doi.org/10.1190/1.1440151
[13]	Zonge, K.L. and Wynn, J.C. (1975) Recent Advances and Applications in Complex Resistivity Measurements. Geophysics, 40, 851-864. https://doi.org/10.1190/1.1440572
[14]	Xu, S.Z. (1988) Selection of Wave Number K in Fourier Inverse Transformation for Point Source and 2-D Electric Field Problems. Computing Techniques for Geophysical and Geochemical Exploration, No. 3, 235-239.
[15]	Cai, J.T., Ruan, B.Y., Zhao, G.Z. and Zhu, G.P. (2007) Two-Dimensional Modeling of Complex Resistivity Using Finite Element Method. Chinese Journal of Geophysics, 50, 1615-1624. https://doi.org/10.1002/cjg2.1182
[16]	Han, K.H. and Kim, J.H. (2002) Quantum-Inspired Evolutionary Algorithm for a Class of Combinatorial Optimization. IEEE Transactions on Evolutionary Computation, 6, 580-593. https://doi.org/10.1109/TEVC.2002.804320