Research on Inversion Resolution for ERT Data and Applications for Mineral Exploration

This paper presents the effects of two factors on electrical resistivity tomography (ERT) resolution using numerical simulations. The results obtained from the first group of synthetic resistivity model numerical simulations with different electrode configurations illustrate that the inversion resolution can be improved if the current electrodes and potential electrodes are installed at a subsurface depth. The results obtained from the second group of synthetic resistivity models with different mesh densities illustrate that the inversion resolution can be improved if the mesh grid density is increased. The increased number of degrees of freedom necessary to fit the data to a prescribed target misfit leads to the inverted imagery deviating from the true resistivity model. Throughout the literature the inversion procedure is implemented in the conjugate gradient relaxation method and the model calculation speed is enhanced using the parallel processing mode. ERT surveys were employed for mineral exploration at Chifeng in Inner Mongolia, China. The inverted image clearly describes the vein behavior down to 200 m in depth in the prospect area and is consistent with the drilling data results.


INTRODUCTION
Electrical resistivity tomography (ERT) is a geophysical technique in which direct current (DC) is injected into the ground between a pair of electrodes with the voltage measured between a second pair of electrodes.ERT surveys have been used in resource exploration (Tang et al. 2007;Yi et al. 2011;Fikos et al. 2012;Hermans et al. 2012b), hydrogeology surveys (Lee et al. 2008;Hermans et al. 2012a;Maiti et al. 2012;Muchingami et al. 2012), engineering geology surveys (Chang et al. 2012;Lehmann et al. 2013;Springman et al. 2013;Perrone et al. 2014) and environmental geology surveys (Martínez-Pagán et al. 2010;Metwaly et al. 2013;Sirhan and Hamidi 2013;Sonkamble 2014).Numerous approaches have been developed by geophysicists including the finite-difference (FD) method (Dey and Morrison 1979a;Mizunaga and Ushijima 1991;Mizunaga et al. 2003;Zhang et al. 2011;Demirci et al. 2012), the finite-element (FE) method (Bing and Greenhalgh 2001;Zhou and Greenhalgh 2002;Zhou et al. 2009;Vachiratienchai et al. 2010;Vachiratienchai and Siripunvaraporn 2013) and the boundary-integral (BI) method (Schulz 1983(Schulz , 1985;;Hvoždara 2012).The electric potential (difference) in twodimensional (2-D) and three-dimensional (3-D) conductivity structures can be obtained using the linear system of equations solution arising from partial differential equation discretization using these numerical approaches.In general, FD methods are much more flexible than BI methods in treating conductivity structures.The FD method calculation requirements are less stringent than those for FE methods.Throughout the literature, FD methods are used for numerical inversion.
The incomplete Cholesky conjugate gradient (ICCG) method (Meijerink and van der Vorst 1977;Poole and Ortega 1987) were employed for the DC resistivity forward problem to solve a large sparse symmetrical positive definite linear system of equations.The conjugate gradient relaxation method (Mackie and Madden 1993;Spitzer 1995;Wang and Mezzatesta 2001) is applied to the resistivity inversion problem to solve a maximum likelihood system of equations rather than solving the system directly using matrix inversion routines.
Based on the principles outlined above a rapid inversion procedure for ERT surveys is developed in this work using the C# language.Using the parallel computation method for the resistivity inversion program, the computational efficiency is significantly improved.ERT method uses multiple current sources at multiple separations to generate multiple sets of ERT data.The parallel methods are perfectly adequate for solving such problems.
This paper focuses on a study of the effects of several factors associated with the inversion resolution for ERT data using synthetic resistivity models with different electrode configurations and different mesh densities.ERT surveys are carried out for mineral exploration in Chifeng, Inner Mongolia, China.Although many other methods have been used to investigate metallic minerals, such as geochemical exploration analysis, surface geology and magnetotelluric (MT) surveying, geochemical exploration analysis cannot provide much detailed information for a large area.Geological methods can provide land surface information in cases in which it is difficult to precisely describe the subsurface resistivity distribution.MT surveying also has a worse anti-jamming ability than DC resistivity methods (Ushijima et al. 1990).By contrast, ERT surveys can provide highresolution images of the subsurface resistivity distribution.ERT is a rapid, effective and economical method for mineral exploration.Moreover, ERT using buried transmitters or receivers using boreholes have become increasingly important exploration tools for deep mineral deposit exploration.

INVERSION STRATEGY
Equations governing the DC response of a point current source in a steady state problem are given by (Telford et al. 1990;Bing and Greenhalgh 2001;Long et al. 2006) where z is the electrical potential (V), I is the source current (A), r r s d -^h is a delta function, r and r s are the observation point and point source current locations (m), respectively, and v is the electrical conductivity (S m -1 ).
There is a singular term defined by Eq. ( 1).The electrical potential is divided into two terms to delete the singular term.The first term ( 0 z ) is the well-known formula for the potential in a uniform material with conductivity 0 v , and the second term ( S z ) is a correction made due to the non-uniformity of a given material (Mizunaga and Ushijima 1991).The potential can therefore be defined as follows: By substituting Eqs.(1) into ( 2) and (3), we obtain the following equation: The resistivity problem discretization discussed here is based on the theory outlined by Dey and Morrison (1979a, b).Using the FD method potentials are calculated for a discrete number of transform variables at the nodes of a quadrilateral element mesh.The FD method is employed for forward modeling, in which the modeling routine accounts for 2-D source current electrodes.The discretization equation is obtained as follows: Boundary conditions are applied at the ground surface, with z = 0; specifically, the Neumann boundary condition represented by Eq. ( 6) is applied.The other model space boundary for inversion is implemented by applying the mixedboundary condition (Dey and Morrison 1979b) indicated by Eq. ( 7): where i is the angle between radial distance r and outward normal n.
Solving for a set of equations normally consumes most of the computing time in the forward problem.Compared with direct methods the ICCG algorithm is a moreefficient iterative solver for large sparse linear equations.In the ICCG algorithm only the nonzero elements are stored and the algorithm converges very quickly (Fu et al. 2002).Therefore, the ICCG algorithm is employed to solve the forward problem.
The smoothness-constrained least-squares optimization method (Ellis and Oldenburg 1994;Rodi and Mackie 2001;Marescot and Loke 2004) is frequently used for the 2-D resistivity problem.The iteratively reweighted smoothness-constrained least-squares optimization method is used for inversion of ERT data.The equation used is where g i is the data misfit vector containing the difference between the observed values and the calculated potential or the apparent resistivity obtained using forward modeling.
i t D is the change in the model parameters for the ith iteration, and i 1 tis the model parameters vector for the previous iteration.J is the Jacobian matrix of partial derivatives and W is the roughness filter.R d and R m are weighting matrices that are introduced so that different elements of the data misfit and model roughness vectors are given equal weights in the inversion process.i m is the damping factor after the ith iteration, as defined by Eq. ( 9), and the initial damping factor is set to a large value ( 0 m ).
The inverse problem calculation efficiency depends upon the forward calculation counts per inversion iteration.There are two methods to solve J for the DC resistivity inversion including directly solving method (Yorkey et al. 1987) and employing a similar procedure for solving J by multiplying an arbitrary vector x (Rodi 1976;Mackie and Madden 1993).We prefer employing Rodi's method over Yorkey's method.Because directly solving method require solving one forward problem for each model parameter for each inversion iteration on one current source, the amount of calculation is great when the grid mesh density is high, so the mesh grid density influence on the efficiency inverse is enormous.However, the approach for computing Jx and J T y requires quasi forward calculation once for inversion iteration on one current source.Rodi's method is more efficient for the resistivity inversion problem than the traditional method for resistivity inverse problem with multiple current sources.Parallel techniques have been employed to solve J, which greatly increases the calculation efficiency.Model roughness must be minimized to suppress any model structure not required by the data.For a 2-D structure with x lying along the strike axis direction a model roughness may be given by (De Groot-Hedlin and Constable 1990) where m is the model parameters vector, ~x 2 is a roughening matrix that differentiates the model parameters of laterally adjacent prisms, and ~z 2 is a roughening matrix that differentiates the model parameters of vertically adjacent prisms.
For solving a set of equations developed using Eq. ( 8) the conjugate gradient relaxation method is employed with implementation in parallel processing mode.The parallel methods are applicable in this ERT method with multiple current sources problem.The forward modeling of different sources can be thrown into the different threads for calculations using a personal micro-computer or workstation.Solving the Jacobi matrix for different sources can be implemented using parallel processing.
The Root mean square (RMS) relative error is defined as shown in Eq. ( 11), where, d obs is measured potential in field study; d forward is calculated by the inverted modeling; M is the number of observations.

SYNTHETIC TEST RESISTIVITY MODELS AND RESOLUTION
This paper mainly studies two factors affecting the inversion resolution including electrode arrangements and mesh grid density for numerical inversion.Numerical simulation for inversion of ERT data from the synthetic resistivity modeling can guide the field study.We then can achieve better inverted imageries at minimal cost for reasonable electrode configurations and other resources.Several synthetic resistivity models are carried out to study the resolution, as described in the following sections.

Double Anomalous Bodies with Comprehensive Configurations Data Sets
The electrode configuration is an important factor affecting inversion resolution.The resolution properties of different data configuration sets are explored.One conductive 10 Ω•m and one resistive 500 Ω•m body are inserted into the 100 Ω•m half-space (Fig. 1).Synthetic resistivity models with four kinds of electrode configurations are performed as shown in Figs.1a, c, e, and g.The first electrode configuration (Fig. 1a) is performed with both current electrodes and potential electrodes installed at the surface.The second electrode configuration (Fig. 1c) is carried out with current electrodes installed in the subsurface while potential electrodes are installed at the surface.Compared to the second electrode configuration installed for the second synthetic resistivity model, the third electrode configuration (Fig. 1e) adds many more potential electrodes in two boreholes.In the fourth electrode configuration (Fig. 1g) the added potential electrodes in two boreholes for the third electrode configuration are replaced with current electrodes.The study region volume is 30 × 20 m.The region is meshed into 30 × 20 equal grids, and the number of nodes for the grid mesh is 31 × 21.The unit grid block spacing in this study model is 1 m.
Images (Figs. 1b,d,f,h,i,and j) were obtained from different electrode data configurations using the same inverse iteration of 20.The four inverted images illustrate that the first image (Fig. 1b) obtained from the first electrode configuration cannot describe the subsurface resistivity distribution well from the true model and it can only resolve the shallow conductive body and not the body in depth.Inverse modeling with the second electrode configuration is illustrated in Fig. 1d.Compared with the Figs.1b and d images the second electrode configuration can achieve a better resolution and discover the resistive body in depth.The images obtained from third and fourth electrode configurations describe the subsurface resistivity distribution well, which is close to the true resistivity model.The interface boundaries are described better than in the first two images obtained from the first two electrode data configurations set inversion.Two vertical cross sections at x = 7.5 m and x = 22.5 m respectively are shown in Figs.1i and j to highlight the inverse resolution capabilities.The two cross sections illustrate that inverted image with the third and fourth electrode configurations can perform better than the first two configurations.The two anomalous bodies (conductive 10 Ω•m body and resistive 500 Ω•m body), present inversion curves for the two latter data configuration sets that are closer to the true resistivity model than the curves from the first two configurations.The more electrodes (current and potential electrodes) installed at subsurface, the better resolution will be achieved.Figure 1k shows the convergence plots for inversion while employing the first electrode configuration, second electrode configuration, third configuration and the fourth electrode configuration, respectively.From the RMS curves with the four kinds of electrode configurations, the inversion reaches convergence with about 1.0% RMS data misfit, presenting satisfying results after about six iterations.

An Extreme Synthetic Resistivity Model with Limited Data Sets
The above described synthetic resistivity model tests show that the electrode configuration is one of the most important factors affecting the inversion resolution.We conclude that current electrodes or potential electrodes installed at subsurface will achieve better subsurface resistivity distribution description.However, a new problem is proposed: which factor is more influential, the current electrodes or the potential electrodes?To address this problem an extreme synthetic resistivity model is employed (see Figs. 2a, c, e, and g).One conductive 10 Ω•m body is inserted into the 100 Ω•m half-space (Fig. 2).The study region volume and the number of grid blocks is the same as described in section 3.1.1.
Four kinds of electrode configurations were performed as shown in Figs.2a, c, e, and g.The first electrode configuration (Fig. 2a) employed multiple pairs of current electrodes and a single pair of potential electrodes.The second electrode configuration (Fig. 2c) is carried out using single pair of potential electrodes installed at the subsurface with potential electrodes installed at the surface.The third (Fig. 2e) and fourth (Fig. 2g) electrode configurations employed multiple pairs of potential electrodes with a single pair of current electrodes.
Images obtained from four extreme electrode data configurations are illustrated in Figs.2b, d, f, and h.All of these configurations used the same inverse iteration of 20.Comparing Figs.2b with d the potential electrode installed at the subsurface can achieve better resolution at the interface boundary, particularly in Fig. 2j.Comparing Figs.2f with h the current electrode installed at the subsurface can achieve better resolution at the interface boundary, particularly in Fig. 2j.These two conclusions are the same as those obtained from section 3.1.1.Comparing Figs.2b with f the potential electrode installed at the subsurface can obtain better resolution, because of better interface boundaries description, as shown in Fig. 2j. Figure 2k is the convergence plot for inversion employing the first electrode configuration, second electrode configuration, third configuration and fourth electrode configuration, respectively.From the RMS curves with the four kinds of electrode configuration the inversion reaches convergence with about 4.0% RMS data misfit and results in satisfactory results after about twelve iterations.

Mesh Grid Density Effects
ERT surveys are normally carried out in the field using certain electrode spacing.This spacing can be employed as the unit grid block spacing for field ERT data inversion.However, mesh refinement is a better choice for obtaining better resolution imaging.Theoretically, the greater the mesh grid density is, the better the inversion resolution will become.There will be many more resistivity blocks used to describe electrical characteristics underground, allowing for better interface boundary description than that obtained using low mesh grid density.Assuming the same data density is used, the higher mesh grid density means that the number of degrees of freedom for the inversion will be greater, which could lead to much more redundant physical structures.Based on the foregoing discussion a group of synthetic resistivity models were employed with different mesh densities for modeling and inversion.
Figure 3a shows the synthetic rectangular low-resistivity models with one conductive 10 Ω•m body inserted into the 100 Ω•m half-space.The region is meshed into 30 × 20 and 15 × 10 for the high and low mesh density models, respectively.The unit grid block spacing's in this study are 1 and 2 m, respectively.Current electrodes are installed at the subsurface while potential electrodes are installed at the surface and subsurface for measuring the electrical potential.
Inverse modeling is illustrated in Figs.3b, c, d, e.Compared with images obtained from two different of mesh density data set inversions, Fig. 3b shows that the interface boundaries can be described better than Fig. 3c, particularly at the center of the rectangular conductive body and the top boundary of the anomalous body.Compared to the imagery (Fig. 3b) obtained from the high density data set inversion, the center of the rectangular conductive body drifts upward in Fig. 3c.There are fewer redundant physical structures appearing along the anomalous body interface.Figure 3f is the convergence plot for inversion when high and low mesh grid density modes are employed, respectively.From the RMS curves with different mesh density models, the inversion reaches convergence with about 0.5% RMS data misfit and produces satisfying results after five iterations.
Based on the aforementioned comparative experiments, we can infer that the resolution will be better and the number of degrees of freedom for the inversion will be greater, leading to many redundant physical structures if the mesh grid density is greater.

COMPARISON OF THE PARALLEL AND NONPARALLEL COMPUTATION
Parallelism has been employed in high-performance computing for many years.For the numerical simulation of ERT with the high data density and high mesh density, the traditional method for solving the inverse problem has difficulty meeting the requirements of large data sets, particularly for inverting real-time images in some projects, such as real-time groundwater flow monitoring and rainfall infiltration process imaging and so on.
Using the Microsoft Visual Studio development environment, the .NET Framework 4 enhances support for parallel programming by providing a new runtime, new class library types, and new diagnostic tools.These features simplify parallel development so that we can scale parallel code in a natural idiom without having to work directly with threads or the thread pool.The code modules that use parallel computing for the inversion procedure include the forward module and Jacobi matrix calculating module.
A synthetic resistivity model (Fig. 1e) is used to invert ERT data for parallel and nonparallel method comparison.The numbers of current and potential electrodes are 900 and 906, respectively.The performance was rated on a Z820 WorkStation with 12 Interl (R) Xeon (R) CPU E5-2620 v2 @2.1 GHz processors under a 64-bit operating system.Using a progressive reduction in the number of CPUs to run the inverse procedure, we obtained different calculation times in 1, 2, 4, 8, 12, 16, 20, 24 threads, respectively (Fig. 4).
Figure 4 illustrates that the model calculation speed is enhanced by parallel processing.From Table 1 when using the parallel method the computation time is significantly reduced with the number of threads increased.However, without using the parallel method, the computation time is reduce by only about 3 hours with the number of threads increased from 1 -24.

ERT INSTRUMENT
The ERT instrument employed for field study includes the transmitter (Fig. 5a) and receiver (Fig. 5b).The DC is injected using a borehole electrode (Fig. 5c) and the voltage is measured using the recording geometry (Fig. 5d).The maximum emissive power is 10 KW and the output current can be set to 10 A (the maximum supply current).The voltage range for the input and output are 20 Vpp and 20 Vpp, respectively.The offset voltage is less than 2 mV.The channel-to-channel cross-talk is better than 80 dB.The determination accuracy is ±0.5% of the measurement range.

CASE STUDY
ERT surveys were carried out in Chifeng, Inner Mongolia, China (see Fig. 6a).The ERT survey line is perpendicular to the metalliferous vein strike according to the geological and well logging data.The ERT survey line is installed in the north-south direction for measuring the potential to inverse subsurface resistivity distribution.Using the limited information obtained from the geological and well logging data is difficult to describe the vein behavior with high resolution and estimate the mineral reserves with high accuracy.However, the ERT method can achieve the subsurface resistivity distribution with high resolution.Moreover, using the inverted imagery the vein scale can be estimated with high accuracy.
The electrode array is shown in Figs.6b, c, d.Current electrodes are located at various depths in the borehole and potential electrodes are located at the surface.Because the study area is in low dip and subdued topography there is no need to correct the topography effect in this case study.There were 98 electric potential difference observations obtained in this field exploration.As previously mentioned, ERT imaging using ERT data inversion can provide highresolution images of the subsurface electrical resistivity distribution, which could be used to describe the vein distribution and estimate mineral reserves.
Imagery obtained from ERT data inversion is shown in Fig. 7a.The study region volume is 700 × 320 m.The region is meshed into 28 × 12 grids, and the unit grid block spacing in this study model is 25 m.The imagery illustrates that the metalliferous vein tendency is southward.The inverted image is in accord with the drilling data determined by the borehole (ZKW7-2), and except for a slight upward bias, the first three metal mine discovery depths agree well with the borehole drilling data (Fig. 7b).Using the pseudo section maps from the observations (Fig. 7d) and calculated apparent resistivity (Fig. 7e), the inverted image is reliable for investigating mineral reserves.

DISCUSSIONS AND CONCLUSIONS
Two influencing factors were investigated to resolve the ERT data inversion by conducting a numerical simulation of a 2-D resistivity inversion.The current and potential electrode effects were examined.The analysis shows that current and potential electrodes located in the subsurface   would achieve better inversion resolution under the same data density.Moreover, more installed electrodes will obtain more information for inversion.The mesh density effect was also investigated.The results obtained from the numerical simulation illustrate that the resolution will be better if the mesh grid density is higher.Additionally, the number of degrees of freedom for the inversion will be greater, which will lead to many more redundant physical structures in the inverted imagery.
Overall, the current and potential electrodes should be located as deep below the surface as possible to achieve better resolution for deep prospecting.The traditional resistivity tomography methods with the surface array are sensitive to anomalous bodies near the surface but not to deep mineral reserves.With the development of new computer technologies current and potential electrodes can be installed at the subsurface using boreholes and the developed inverse procedure can be applied to solving such problems.If topographic conditions permit and the survey instruments can be applied in the required arrangements, the current and potential electrodes should be located as deep below the surface as possible.
Throughout literature the inverse procedure is implemented in parallel processing mode because there is no procedure to solving such problems with low mesh density and low data density.However, working with large data sets, particularly for inverting real-time images in some projects such as real-time large scale groundwater flow monitoring, rainfall infiltration process imaging and so on is difficult.Comparing parallel and nonparallel inversion methods, the speed of model calculation is improved with parallel processing mode implementation.
ERT surveys were carried out for mineral exploration in Chifeng, Inner Mongolia, China.The imagery obtained from ERT data inversion supports drilling data and achieved high-resolution subsurface resistivity image distribution.ERT surveying offers its own advantages for mineral exploration that cannot be achieved using drilling data and geological methods.
According to the above mentioned analysis an ERT survey can allow for high-resolution inverse modeling.Therefore, ERT surveying can be applied in other fields, allowing for static and dynamic monitoring of variations in subsurface resistivity distribution characteristics.Examples include monitoring the behavior of groundwater and chemical pollution or detecting landslides.

Fig. 3 .
Fig. 3. True resistivity models (a).Red dots and blue triangles represent current electrodes and potential electrodes, respectively.Tomograms resulting from different mesh density data sets (b), (c).(d) and (e) show vertical and horizontal cross sections at x = 15 m and z = 5.5 m, respectively.(f) RMS of inversions with different electrode configurations.

Fig. 4 .
Fig. 4. Comparison of the inversion calculation time using parallel and nonparallel processing.

Fig. 6 .
Fig. 6.Location and topographical map of the research area.(a) Location of the study site.(b) Google Earth images of survey stations; the yellow label represents the location of exploratory borehole (ZKW7-2), whereas the red labels represent locations of potential electrodes.(c) Topography of the research area; the ERT survey line is perpendicular to the strike of the metalliferous vein.(d) Elevation of the electrodes; the current electrode C 0 and potent electrode P 0 are 1.5 km away from the well, and the current through C 0 , P 0 forms the return circuit.

Fig. 7 .
Fig. 7. Inverted imagery from inversion of ERT data in Inner Mongolia (a), drilling data (b) and iterations of inversion (c), pseudo section maps from the observations (d), and calculated apparent resistivity (e).

Table 1 .
Computing time with different methods.