Electrical Resistivity Tomography: A Subsurface-Imaging Technique

Abstract

Electrical resistivity tomography (ERT) is a popular geophysical subsurface-imaging technique and widely applied to mineral prospecting, hydrological exploration, environmental investigation and civil engineering, as well as archaeological mapping. This chapter offers an overall review of technical aspects of ERT, which includes the fundamental theory of direct-current (DC) resistivity exploration, electrode arrays for data acquisition, numerical modelling methods and tomographic inversion algorithms. The section of fundamental theory shows basic formulae and principle of DC resistivity exploration. The section of electrode arrays summarises the previous study on all traditional-electrode arrays and recommends 4 electrode arrays for data acquisition of surface ERT and 3 electrode arrays for cross-hole ERT. The section of numerical modelling demonstrates an advanced version of finite-element method, called Gaussian quadrature grid approach, which is advantageous to a numerical simulation of ERT for complex geological models. The section of tomographic inversion presents the generalised standard conjugate gradient algorithms for both the l1- and l2-normed inversions. After that, some synthetic and real imaging examples are given to show the near-surface imaging capabilities of ERT.

1. Introduction

Direct-current (DC) resistivity exploration is a traditional geophysical method. It employs two electrodes to inject electric current into the ground and other two electrodes to measure the electric potential difference. The measurements are often carried out along a line or in an area on the earth surface, and then the observed potential differences are converted into sounding curves or pseudo-sections of apparent resistivities, which indicate the resistivity changes of subsurface rocks. Analyses of these data enable us to find the underground resistivity anomalies or outline the subsurface geological structure. With development of computer technology and numerical computational techniques, accurate numerical simulations of subsurface electrical field and acquiring a large amount of data in fields become possible [1, 2, 3], so that the traditional DC resistivity exploration was developed to a computerised geotomography technique, called electrical resistivity tomography (ERT), which employs a multielectrode equipment or system to automatically acquire a large number of data [4, 5] and applies a computer software to the reconstruction of subsurface resistivity structure with the observed data [6, 7, 8, 9, 10]. Due to its conceptual simplicity, low equipment cost and ease of use, ERT is now widely applied in mineral exploration, civil engineering, hydrological prospecting and environmental investigations, as well as archaeological mapping [11]. This chapter provides an overall review of ERT techniques, which consists of four sections: (1) fundamental theory, (2) electrode arrays, (3) numerical modelling and (4) tomographic inversion. In each section, diagrams and formulations are used to illustrate basic concepts and principles of ERT techniques. Some synthetic experiments and practical imaging applications are also given to show the imaging capability of ERT.

2. Fundamental theory

According to the continuity of electrical current, the following integral equation is satisfied at any point in a conductive medium:

$${∯}_{\mathrm{\Gamma }}\mathit{J}\cdot d\mathrm{\Gamma }=-I,$$E1

where Γ is a full or half spherical surface that encloses an electrode that injects electric current J with a magnitude of I (Figure 1). According to Ohm’s law:

$$\mathit{J}=\mathrm{\sigma }\mathit{E}=\mathit{-}\sigma \nabla U,$$E2

and the property of the delta-function δ(x-x_s), Eq. (1) gives

$$-{∯}_{\mathrm{\Gamma }}\sigma \nabla U\cdot d\mathrm{\Gamma }=-I{\iiint }_{\mathrm{\Omega }}\mathrm{\delta }\left(\mathit{x}\mathbf{-}{\mathit{x}}_s\right)d\mathrm{\Omega }$$E3

Figure 1.

Electric current density J around an electrode in (a) a full space and (b) a half-space.

Here σ is the conductivity of medium, U is the electric potential, Ω represents the volume of the spherical surface Γ, and x_s = (x_s,y_s,z_s) is the location of the current electrode. Applying the divergence theorem to Eq. (3) leads to

$${\iiint }_{\mathrm{\Omega }}\left[\nabla \cdot \left(\sigma ,\nabla ,U\right)+I\mathrm{\delta }\left(\mathit{x}\mathbf{-}{\mathit{x}}_{\mathit{s}}\right)\right]=0.$$E4

Note that Eq. (4) satisfies everywhere in a medium, so that the following governing equation of electric field is obtained:

$$\nabla .\left(\sigma ,\nabla ,U\right)=-I\mathrm{\delta }\left(\mathit{x}-{\mathit{x}}_s\right),.\mathit{x}\in \Omega .$$E5

In general, the conductivity σ(x) = σ(x,y,z) changes with three coordinates so that it defines a 3D geological model, and the electric potential U=U(x,y,z) becomes a 3D function of a given current injection I$\mathrm{\delta }\left(\mathit{x}-{\mathit{x}}_s\right)$ and conductivity σ(x) (Figure 2a). If the conductivity is a constant with the y-coordinate, i.e. σ(x) = σ(x,z), it defines a 2D geological model (Figure 2b), but the electric potential U(x,y,z) is still a 3D function of the coordinates due to the point current injection I$\mathrm{\delta }\left(\mathit{x}-{\mathit{x}}_s\right)$. Applying cosine Fourier transform F_c{⋅} to the y-coordinate of Eq. (5), the governing equation is changed into

$$\nabla .\left(\sigma \nabla \stackrel{\sim }{\mathrm{U}}\right)+k_y^2\stackrel{\sim }{\mathrm{U}}=-I\mathrm{\delta }\left(\mathit{x}-{\mathit{x}}_s\right)/2,\mathit{x}\in \Omega ,$$E6

Figure 2.

Sketch of (a) a 2D and (b) a 3D geological model defined by conductivity σ(x,z) and σ(x,y,z), respectively.

which is often named 2.5D governing equation. Here k_y is the wavenumber, and $\stackrel{\sim }{\mathrm{U}}\left(x,{\mathrm{k}}_y,z\right)=F_c\left\{\mathrm{U}\left(x,y,z\right)\right\}$ becomes the spectrum of electric potential in the wavenumber domain. The coordinate vector x = (x,z) and gradient ∇ = (∂_x, ∂_z) are the 2D versions. If a geological model σ(x,y,z) or σ(x,z) and current injection I$\mathrm{\delta }\left(\mathit{x}-{\mathit{x}}_s\right)$ are given, solving the governing equation (5) or (6), one obtains electrical potential U(x,y,z) for a 3D conductivity model or its spectrum $\stackrel{\sim }{\mathrm{U}}\left(x,k_y,z\right)$ for a 2D conductivity model and then performs the inverse cosine Fourier transform U(x, y, z) = $F_c^{-1}\left\{\stackrel{\sim }{\mathrm{U}}\left(x,{\mathrm{k}}_y,z\right)\right\}$ to obtain electric potential U(x,y,z). These computations are called forward modelling. In theoretical computation or numerical forward modelling, Green’s function G or $\stackrel{\sim }{\mathrm{G}}$—the electric potential response to a unit current injection (I = 1 Am)—is often applied, so as to remove the magnitude of electric current I in Eqs. (5) and (6). One may focus on the computation of Green’s function, and then the electric potential U or $\stackrel{\sim }{\mathrm{U}}$ is computed by the multiplication of Green’s function G with the magnitude of practical electric current I, i.e. U = I⋅G.

The simplest geological model is a homogenous half-space. Applying Eq. (3) (equivalent to Eq. (5)) to a constant medium, the surface integral is calculated by

$${∯}_{\mathrm{\Gamma }}\sigma \nabla U\cdot d\mathrm{\Gamma }=\sigma \nabla U\left(2\mathrm{\pi }{r}^2\right)=-I,$$E7

which gives the electric potential at distance r from a current source:

$$\mathrm{U}=\frac{I}{2\mathrm{\pi \sigma }r}=\frac{\mathrm{\rho }I}{2\mathrm{\pi }r},$$E8

where ρ = 1/σ is resistivity. In practice, to inject electric current into the ground, a pair of current electrodes C1 and C2 must be employed; one is positive (+I) and another is negative (-I). Thus, the electric potential at a point P with a pair of current electrodes is calculated by.

$$\mathrm{U}=\frac{\mathrm{\rho }I}{2\mathrm{\pi }}\left(\frac{1}{r_{C1P}}-\frac{1}{r_{C2P}}\right).$$E9

Here, r_C1P and r_C2P are distances of the observed point P to two current electrodes C1 and C2. To measure the potential on the earth surface, a pair of potential electrodes P1 and P2 is also required (Figure 1). According to Eq. (9), one has the following expression for the electric potential difference between two potential electrodes:

$$\mathrm{\Delta }U=\frac{\mathrm{\rho }I}{2\mathrm{\pi }}\left[\left(\frac{1}{r_{C1P1}}-\frac{1}{r_{C2P1}}\right)-\left(\frac{1}{r_{C1P2}}-\frac{1}{r_{C2P2}}\right)\right],$$E10

from which apparent resistivity ρ_a is defined as follows:

$${\mathrm{\rho }}_a=K\frac{\mathrm{\Delta U}}{I}=K\mathrm{\Delta G},$$E11

where ΔG = ΔU/I and K is called geometry factor of electrode array given by

$$K=2\mathrm{\pi }{\left[\left(\frac{1}{r_{C1P1}}-\frac{1}{r_{C2P1}}\right)-\left(\frac{1}{r_{C1P2}}-\frac{1}{r_{C2P2}}\right)\right]}^{-1},$$E12

which depends on the positions of four electrodes. Different layouts of four electrodes have variable geometry factors and are often called electrode arrays. In the traditional electrode arrays, C2 and P2 may be set up very far away from C1 and P1, so that C2 and P2 are treated as remote electrodes in theory, and distances $r_{C1P2},r_{C2P2}$, $r_{C2P1}$, and $r_{C2P2}$ are supposed to be infinite (∞) in Eq. (12). In these cases, one can find that geometry factors are still applicable for these electrode arrays that involve one or two remote electrodes, which are named pole-pole, pole-dipole and dipole-pole arrays.

If subsurface resistivity is homogenous (ρ₀), Eq. (11) shows that no matter which electrode array is used, apparent resistivity is constant (ρ_a = ρ₀). Otherwise, ρ_a indicates resistivity variation of the underground. For a certain range of apparent resistivity ρ_a, Eq. (11) also reveals that the geometry factor is inversely proportional to the potential difference ΔG = ΔU/I. It implies that a big value of geometry factor K will cause a small reading of the electric potential difference ΔG in fields. Such a small reading is easily obscured or contaminated by background noise. For a good data quality, one should avoid very big values of geometry factors in data acquisition for ERT. Therefore, Eqs. (11) and (12) are the fundament formulae of the traditional DC electrical resistivity exploration.

3. Electrode arrays

In order to obtain apparent resistivity in fields, many electrode arrays were developed in the traditional DC resistivity exploration. In principle, ERT requires a high data density and good coverage of the earth surface for high-resolution images of subsurface targets. Dahlin and Zhou [11] carried out synthetic experiments of ERT using 10 electric arrays and compared their imaging results for four geological models: a buried channel, a narrow dike, dipping blocks and waste ponds. They demonstrated that two three-electrode arrays (pole-dipole and dipole-pole) and three four-electrode arrays (dipole–dipole, Schlumberger and gradient arrays) produce satisfactory images of the subsurface targets. However, due to the use of remote electrodes, the three-electrode arrays are rarely applied for ERT in practice; thus, the four-electrode arrays become popular. Particularly, gradient array [12] is well suited for multichannel data acquisition and can significantly increase the speed of data acquisition in the field, and at the same time, it gives higher data density and lower sensitivity to noise than dipole-dipole array. Figure 3 shows three four-electrode arrays (upper row) and their pseudo-sections of data points $\left(x_{{\mathrm{\rho }}_a},z_{{\mathrm{\rho }}_a}\right)$ (middle row) and their geometry factors (bottom row) for a layout of total 81 electrodes on the earth surface. Figure 3 shows that dipole–dipole array has the biggest range of geometry factor among these three-electrode arrays, and gradient array performs the smallest range, so that it is well suitable for data acquisition with a high data density and good data quality. Zhou and Dahlin [13] based on many sets of field data obtained with different electrode arrays and found that the measured potential errors depend on the reading of potential difference, ΔG = ΔU/I, and the measured datum in field may be expressed by

$$\mathrm{\Delta }{\mathrm{G}}^{\ast }=\mathrm{\Delta G}\left[1+\mathrm{R}{\left(\frac{\alpha }{\mathrm{\Delta G}}\right)}^{\beta },/,100\right],$$E13

where ΔG^* denotes a noise-contaminated datum, R is a random number in (−1,1) and α and β are constants. Due to having a big range of geometry factor (see Figure 3), the dipole-dipole array is more sensitive to noise than Schlumberger and gradient arrays, so that one must pay much attention to the noise contamination when conducting dipole-dipole measurements [13]. According to the relationship between the geometry factor K and reading of potential difference ΔG, Eq. (13) may be applied to predications of noise-contaminated data for different electrode arrays [11].

Figure 3.

Dipole-dipole, Schlumberger and gradient arrays (upper row) for ERT data acquisition and examples of their pseudo-sections of data points ($x_{{\mathrm{\rho }}_a},z_{{\mathrm{\rho }}_a})$(middle row) and geometry factors (bottom row).

Reviewing Figure 3, one may find that gradient array only uses the potential electrodes between C1 and C2 to measure the potential differences. It does not do the measurements at the potential electrodes outside of C1 and C2. It means that gradient array actually misses the outside data for every pair of current electrodes. To obtain a better data coverage and complete subsurface information, the outside potential electrodes of each pair of C1 and C2 should be employed in the gradient measurements. These additional data apparently complement to common gradient measurements and may improve ERT imaging. Accordingly, a new electrode array is naturally formed and called full-range gradient array shown in Figure 4a. The difference from the common gradient array is that the new electrode array uses not only the potential electrodes between C1 and C2 but also the outside potential electrodes of every pair of C1 and C2, if their geometry factors fall in a reasonable range for the data acquisition. Figure 4b and c shows the pseudo-section points and geometry factors for a layout of total 81 electrodes. Comparing with Figure 3, the full-range gradient array does improve the data coverage, and its geometry factors are controlled in a reasonable range.

Figure 4.

Full-range gradient array (a) for ERT data acquisition and examples of its pseudo-section of data points $\left(x_{{\mathrm{\rho }}_a},z_{{\mathrm{\rho }}_a}\right)$ (b) and geometry factors (c) with a layout of total 81 electrodes.

If there are boreholes in a field, cross-hole ERT may be carried out to image the geological structure between the boreholes. Zhou and Greenhalgh [14] investigated all possible electrode arrays for cross-hole ERT data acquisition and found that the electrode arrays of pole-pole (A-M), pole-bipole (A-MN), bipole-pole (AB-M), and bipole-bipole (AM-BN) with their multi-spacing cross-hole profiling and scanning surveys are useful for cross-hole ERT. Here the capital letters A and B stand for two current electrodes, and M and N denote two potential electrodes. These cross-hole electrode arrays are shown in Figure 5 with their sensitivity functions in backgrounds [15]. They also found that the electrode arrays which have either both current electrodes or both potential electrodes in the same borehole, e.g. A-MN, AB-M and AB-MN, have a singularity problem in data acquisition (geometry factor goes to infinite so that apparent resistivity and pseudo-reaction are not applicable), namely, zero readings of the potential or potential difference in cross-hole measurements, so that the potential data are easily obscured by background noise and their images are inferior to those from other cross-hole electrode arrays. The data having the singularity problem may be predicated by zero values of the inverse geometry factors, which should be avoided in cross-hole data acquisition. Therefore, A-M, AM-B, and AM-BN with multi-spaces are recommended for cross-hole ERT.

Figure 5.

Electrode arrays for cross-hole ERT data acquisition. A and B stand for two current electrodes. M and N denote two potential electrodes. The background contours are the sensitivity functions of the electrode array.

4. Numerical modelling

To compute theoretical electric potential U(x,y,z) in a complex 2D or 3D geological model, one has to solve the governing equation (6) or (5) using a numerical approach and a computer. For computational efficiency, an underground artificial boundary Γ₀ is required to truncate an infinite geological model, so that the boundary condition on Γ₀ must be introduced to the governing equations. Applying Eq. (7) to the artificial boundary Γ₀ and zero normal component (J⋅n = 0) of electric current density J on the earth surface, one may find that the electric potential U holds

$$\sigma \nabla U\cdot \mathbf{n}+\mathrm{\nu U}=0,$$E14

where n is the normal vector of the artificial boundary Γ₀ and ν is calculated by ν = (r⋅n)/r². Eq. (14) is the 3D boundary condition on Γ₀. Similarly, one may find the boundary condition for a 2D geological model [16]:

$$\sigma \nabla \stackrel{\sim }{\mathrm{U}}\cdot \mathbf{n}+\mathrm{\nu }\stackrel{\sim }{\mathrm{U}}=0.$$E15

Adding the artificial boundary conditions Eqs. (14) and (15) to Eq. (5) and (6), numerical modelling becomes to solve the following definite governing equation:

$$3\mathrm{D}:\left\{\begin{array}{c}\nabla \cdot \left(\sigma \nabla G\right)=-\delta \left(\mathit{x}-{\mathit{x}}_{\mathit{s}}\right),\mathit{x}\in \mathrm{\Omega },\\ \sigma \nabla G\cdot \mathbf{n}+v\mathrm{G}=0,\mathit{x}\in {\mathrm{\Gamma }}_0\end{array}\right.$$E16

or

$$2\mathrm{D}:\left\{\begin{array}{c}\nabla \cdot \left(\sigma \nabla \tilde{\mathrm{G}}\right)+{\mathrm{k}}_y^2\tilde{\mathrm{G}}=-\delta \left(\mathit{x}-{\mathit{x}}_{\mathit{s}}\right)/2,\mathit{x}\in \mathrm{\Omega },\\ \sigma \nabla \tilde{\mathrm{G}}\cdot \mathbf{n}+v\tilde{\mathrm{G}}=0,\mathit{x}\in {\mathrm{\Gamma }}_0.\end{array}\right.$$E17

The numerical approach to the definite governing equations is called numerical modelling. The most popular numerical approaches are finite-difference and finite-element methods. The former is simple and straightforward, due to directly applying a finite-difference formula to the derivatives of the gradient ∇, so that the definite governing equations at all points in the model domain Ω and on the artificial boundary Γ₀ are discretised and assembled into a linear equation system [16, 17, 18, 19, 20]. The latter converts the definite governing equation into an integral equation by weighting residual principle or variational principle [21] and then carries out numerical volume integration to produce a linear equation system [22, 23, 24, 25, 26]. Therefore, no matter which approach is applied, the numerical modelling becomes to solve a linear equation system. For an instant, the following paragraphs briefly show an advanced 2D version of the finite-element method, called Gaussian quadrature grid (GQG) approach, because of its advantages over finite-difference method and other schemes of finite-element approaches, e.g. high accuracy and easy implementation of discretisation of a geological model having arbitrary free-surface topography. For a 3D numerical modelling, one may follow the 2D procedures.

Weighting residual principle is to calculate the following integral of the governing equation in Eq. (17):

$${\iint }_{\mathrm{\Omega }}\mathrm{W}[\nabla \cdot (\sigma \nabla \stackrel{\sim }{\mathrm{G}})+k_y^2\stackrel{\sim }{\mathrm{G}}+\mathrm{\delta }\left(\mathit{x}\mathbf{-}{\mathit{x}}_{\mathit{s}}\right)/2]d\mathrm{\Omega },$$E18

where W is an arbitrary weighting function. Applying the divergence theorem to the above and then submitting the artificial boundary condition yields

$${\iint }_{\mathrm{\Omega }}\left(\sigma \nabla W\cdot \nabla \stackrel{\sim }{\mathrm{G}}+\mathrm{\sigma }{k}_y^2\mathrm{W}\stackrel{\sim }{\mathrm{G}}\right)d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_0}\mathrm{\nu W}\stackrel{\sim }{\mathrm{G}}d\mathrm{\Gamma }=\mathrm{W}\left({\mathit{x}}_{\mathrm{s}}\right)/2.$$E19

In order to calculate the integrals, the model domain Ω is divided into a set of the no-overlap subdomains {Ω_e, e = 1,2,…,N_e} that matches the free-surface topography and subsurface interfaces. In each subdomain, Gaussian abscissae {(x_α, z_β), α,β = 1,2,…,N_G} are employed (see examples shown in Figure 6) to discretise the subdomain and then Lagrange interpolation:

$$\stackrel{\sim }{\mathrm{G}}={\sum }_{i,j=1}^{N_G}{l}_i\left(x\right),l_j\left(z\right){\stackrel{\sim }{\mathrm{G}}}_{\mathit{ij}}^{\left(e\right)},$$E20

Figure 6.

2D and 3D Gaussian quadrature grids for numerical modelling: (a) a 2D model, 8 × 10 subdomains and 7 × 7 Gaussian abscissae in each subdomain and (b) a 3D model, 10 × 10 × 10 subdomains and 5 × 5 × 5 Gaussian abscissae in each subdomain.

and Gaussian weights {w_αβ} are applied to calculation of the submain integrals. Consequently, Eq. (19) becomes

$$\begin{array}{l}{\displaystyle \sum }_e\{{\displaystyle \sum }_{\text{α},\text{β},i,j}{w}_{\alpha \beta }{\text{σ}}_{\alpha \beta }\left[l_i^{\prime }\left(x_{\text{α}}\right){\text{δ}}_{j\beta }{\partial }_x{\text{W}}_{\alpha \beta }+{\text{δ}}_{i\text{α}}{l}_j^{\prime }\left(z_{\text{β}}\right){\partial }_z{\text{W}}_{\alpha \beta }+{\text{k}}_y^2{\text{W}}_{\alpha \beta }{\text{δ}}_{i\text{α}}{\text{δ}}_{j\beta }\right]{\tilde{\text{G}}}_{\alpha \beta }^{\left(e\right)}\hfill \\ +{\displaystyle \sum }_{\text{α}}{\text{ν}}_{\text{α}}{\text{W}}_{\text{α}}{\tilde{\text{G}}}_{\text{α}}^{\left(e\right)}\}=\text{W}\left(x_s\right)/2.\hfill \end{array}$$E21

Here, ${\stackrel{\sim }{\mathrm{G}}}_{\mathrm{\alpha \beta }}^{\left(e\right)}$ are the discrete values of the electric potential spectra $\stackrel{\sim }{\mathrm{G}}\left(x,{\mathrm{k}}_y,z\right)$ in the subdomain Ω_e. Choosing the Lagrange basis polynomials as the weighting functions, i.e. {${\mathrm{W}}_{\mathit{pq}}=l_p\left(x\right)l_q\left(z\right),\forall p,q\}$, Eq. (21) is changed into

$$\begin{array}{l}{\displaystyle \sum }_e\{{\displaystyle \sum }_{\text{α},\text{β},i,j}{w}_{\alpha \beta }{\text{σ}}_{\alpha \beta }\left[l_i^{\prime }\left(x_{\alpha }\right)l_p^{\prime }\left(x_{\alpha }\right){\text{δ}}_{j\beta }{\text{δ}}_{q\beta }+l_j^{\prime }\left(z_{\beta }\right)l_q^{\prime }\left(z_{\beta }\right){\text{δ}}_{i\alpha }{\text{δ}}_{p\alpha }+{\text{k}}_y^2{\text{δ}}_{p\alpha }{\text{δ}}_{q\beta }{\text{δ}}_{i\alpha }{\text{δ}}_{j\beta }\right]{\tilde{\text{G}}}_{\alpha \beta }^{\left(e\right)}\hfill \\ +{\displaystyle \sum }_{\gamma }{\text{ν}}_{\gamma }{l}_p\left(x_{\gamma }\right)l_q\left(z_{\gamma }\right){\tilde{\text{G}}}_{\gamma }^{\left(e\right)}\}=l_p\left(x_s\right)l_q\left(z_s\right)/2,\hfill \end{array}$$E22

which can be rewritten in the matrix form for all points (p,q) and (i,j) in Ω:

$$\mathbf{M}\stackrel{\sim }{\mathbf{G}}={\mathbf{b}}_s,$$E23

where M is a square symmetric matrix assembled by the coefficients of ${\stackrel{\sim }{\mathrm{G}}}_{\mathrm{\alpha \beta }}^{\left(e\right)}$ in Eq. (22). $\stackrel{\sim }{\mathbf{G}}$ and b_s are two vectors: the former consists of the discrete values of $\stackrel{\sim }{\mathrm{G}}$ at all points of the Gaussian quadrature grid, and the latter is a zero vector except for the component of 1/2 at the electric current source located at x_s. Eq. (22) shows that the matrices M and $\stackrel{\sim }{\mathbf{G}}$ depend on the wavenumber k_y and conductivity model σ_αβ. The wavenumber can be predicated by the minimum and maximum spaces of electrodes [17]. The linear equation system may be efficiently solved for multiple current sources B = {b₁,b₂,…,b_N} by the Cholesky decomposition [27]. After obtaining the spectra $\stackrel{\sim }{\mathbf{G}}$, the theoretical electric potentials are obtained by inverse cosine Fourier transform $\mathbf{U}\left(x,y,z\right)=I\cdot F_c^{-1}\left(\stackrel{\sim }{\mathbf{G}}\right)$.

From Eq. (22), one can find that the model conductivity σ(x,z) is discretised by σ_αβ at the Gaussian abscissae instead of finite elements. The subdomain integrals are calculated by the weights w_αβ to the discrete integrands. The dense abscissae and variable σ_αβ may describe the details of a complex geological model and generate high accurate solutions, so that sizes of the subdomains are not necessarily small but shaped to match the free-surface and subsurface interfaces. Apparently, accuracy of the numerical modelling depends on the number of the Gaussian abscissae in the subdomains. The more abscissae are applied, the more accurate results are yielded but cost more computer time. For efficient and accurate modelling, the minimum electrode space of electrode array may be chosen for the size of the subdomain, to which five Gaussian abscissae is applied to produce satisfactory solutions of numerical modelling. As examples, Figure 7 shows the pseudo-sections of dipole–dipole, Schlumberger and gradient arrays for a layout of total 141 electrodes over an anticline model. These results show that similar pseudo-sections of apparent resistivity are obtained by using the three electrode arrays. From the above analysis, one can see that the finite-element approach does not calculate the high-order derivatives for the governing equation, but any finite-difference scheme does. Therefore, finite-element approach is often called a ‘weak’ solution of the governing equation against a ‘strong’ solution obtained by finite-difference method.

Figure 7.

Numerical modelling for (a) an anticline model and the apparent resistivity pseudo-sections of (b) dipole-dipole, (c) gradient and (d) Schlumberger arrays. The discrete anticline model is also given in the pseudo-sections.

5. Tomographic inversion

Tomographic inversion is to reconstruct the geological model that offers synthetic data matching with observed data. Due to incompleteness and noise contamination of observed data, the model reconstruction is ill posed (multiple solutions). Therefore, tomographic inversion is often defined as an optimisation of data fittingness with regularisation of the model [28], e.g. a generalised objective function is applied:

$$\mathrm{\Phi }\left(\mathbf{m}\right)={‖{\mathbf{d}}_{\mathit{ob}}-{\mathbf{d}}_{\mathit{syn}}\left(\mathbf{m}\right)‖}_{{\mathit{W}}_d}^{l_p}+\mathrm{\lambda }{‖\mathbf{m}‖}_{{\mathit{W}}_{\mathit{m}}}^{l_p},$$E24

where d_ob is a vector of observed data, which are either apparent resistivities or potential differences measured by different electrode arrays. d_syn(m) stands for synthetic data and is calculated by the finite-difference or finite-element method for a guessed geological model m, which consists of discrete conductivities or resistivities. λ is a regularisation parameter that plays a trade-off role between the data fittingness (the first term in Eq. (24)) and model smoothness (the second term in Eq. (24)). ${‖\cdot ‖}_W^{l_p}$ stands for the weighted l_p-norm with a weighting matrix W, e.g. ${‖\mathrm{\epsilon }‖}_W^2={\mathrm{\epsilon }}^T\mathit{W}\mathrm{\epsilon }$ and ${‖\mathrm{\epsilon }‖}_W^1={\sum }_{\mathit{i}}{W}_{\mathit{i}}\mid {\mathrm{\epsilon }}_i\mid$ are the weighted l₂-norm (generalised least square) and the weighted absolute norm. Figure 8 illustrates the difference between a l₂-normed and l₁-normed objective functions. W_d and W_m are weighting matrices to data and model, respectively. It is common to choose a combination of finite-difference operators for W_m, e.g. ${\mathit{W}}_{\mathrm{m}}={\mathrm{\lambda }}_0\mathit{I}\mathit{+}{\mathrm{\lambda }}_1{\mathbf{D}}_x^T{\mathbf{D}}_x\mathit{+}{\mathrm{\lambda }}_2{\mathbf{D}}_y^T{\mathbf{D}}_y\mathit{+}{\mathrm{\lambda }}_3{\mathbf{D}}_z^{\mathit{T}}{\mathbf{D}}_z$ [29]. Here I is a unit matrix, and D_x, D_y and D_z are the finite-difference operators in the x-, y- and z-directions. λ_k (k = 0,1,2,3) are constants and called extensional regularisation parameters used for searching for the smoothest model in three directions. Therefore, tomographic inversion becomes solving the following optimization problem:

$${\mathbf{m}}^{\ast }=min\left\{\mathrm{\Phi }\left(\mathbf{m}\right)\right\}.$$E25

Figure 8.

An l₂-normed (left) and an l₁-normed (right) objective function. Here I is an unit matrix.

To do so, a global or a local search may be applied to Eq. (25) [30], but the global search is extraordinarily computer time-consuming if m has a larger dimension [31]. Therefore, the local search of the standard conjugate gradient method is commonly applied for tomographic inversion. Figure 9 gives a flowchart of the conjugate gradient algorithm and shows that the gradient ∇Φ(m_i) and the Hessian matrix H(m_i) of the objective function are required. Applying the linearised approximation to the synthetic data, d_syn(m_i + 1) ≈ d_syn(m_i) + ($\frac{\partial {\mathbf{d}}_{\mathit{syn}}}{\partial \mathbf{m}}){}^T$(m_i + 1-m_i), and then substituting it for Eq. (24), one can obtain the gradient and the Hessian matrix [28]:

$$\nabla \Phi \left({\mathbf{m}}_i\right)={\left(\frac{\partial \mathbf{d}}{\partial {\mathbf{m}}_i}\right)}^T{\mathbf{W}}_d\left[{\mathbf{d}}_{\mathit{ob}}-{\mathbf{d}}_{\mathit{syn}}\left({\mathbf{m}}_i\right)\right]+\mathrm{\lambda }{\mathbf{W}}_{\mathrm{m}}\left({\mathbf{m}}_i-{\mathbf{m}}_0\right),$$E26

$$\mathbf{H}\left({\mathbf{m}}_i\right)={\left(\frac{\partial \mathbf{d}}{\partial {\mathbf{m}}_i}\right)}^T{\mathbf{W}}_d\frac{\partial \mathbf{d}}{\partial {\mathbf{m}}_i}+\mathrm{\lambda }{\mathbf{W}}_{\mathrm{m}},$$E27

Figure 9.

Flowchart of a standard conjugate gradient algorithm for tomographic inversion. Here m₀ is an initial model, ε is a small value, maxITS is the maximum iterations, and H(m_i) is the Hessian matrix of object function Φ(m_i).

for the l₂-normed objective function.

Note that the gradient ∇Φ(m_i) of the l₁-normed objective function has multiple values at zero misfit Φ(m_i) = 0 (see Figure 8), but for any non-zero misfit Φ(m_i) ≥ ε₀ (ε₀ is a very small value), the computations of the gradient ∇Φ(m_i) and the Hessian matrix H(m_i) become simple, and they are similar to Eqs. (26) and (27) except for the weighting matrices W_d and W_m. Therefore, as Φ(m_i) < ε₀, the l₁-normed gradient ∇Φ(m_i) and the Hessian matrix H(m_i) may be replaced with Eqs. (26) and (27), and as Φ(m_i) ≥ ε₀ the l₁-normed gradient and the Hessian matrix are directly applied to the standard conjugate gradient algorithm. Therefore, the l₁-normed inversion and the l₂-normed inversion are implemented with an almost same algorithm. However, their inversion results may be quite different in the case that the observed data have outliers. Many synthetic and practical experiments have shown that the l₁-normed inversion is less sensitive to the outliers [13]. If the data have no outliers or high qualities, two inversions converge [11].

Figure 10 shows synthetic experiments for imaging the anticline structure shown in Figure 7a with dipole–dipole, Schlumberger and gradient arrays. The l₂-normed inversions were conducted with the apparent resistivities shown in Figure 7b–d, which were computed with the GQG software. The commercial software RES2DINV was applied to the inversions. From these results, one can see that the dipole-dipole and gradient arrays yield competitive images of the anticline structure. Figure 11 gives real applications of surface and cross-hole ERT conducted at a rural site in Australia, where there are many existing drill boreholes which can be used for examinations of both surface and cross-hole ERT experiments. From these boreholes and logging data, the subsurface rocks and main structures were well known (shown in Figure 11b). The main objective of this work was to use the borehole geological and logging information and exam the imaging capabilities of surface and cross-hole ERT in this area, particularly for mapping the base of alluvial overburden and the base of pisolite, as well as predicting clay contamination within pisolite. Surface data acquisition with dipole-dipole and Schlumberger arrays was conducted along a 750 m line, and two pairs of cross-hole ERT on the same line were also carried out for details of the geological structure between the boreholes. Figure 11a shows the integrated cross-hole ERT results with the surface ERT imaging. It gives the resistivity structure along the line from the surface and cross-hole ERT. Figure 11b is the geological section from the existing boreholes and logging data along the same line. Comparing the resistivity imaging results with the geological section, one can see that the integrated ERT results well map the base of alluvial overburden and the base of pisolite. The clay contamination within pisolite is also shown in cross-hole ERT images from the horizontal distance 200–300 m.

Figure 10.

Tomographic inversion of (a) dipole-dipole, (b) Schlumberger and (c) gradient arrays for the anticline model shown in Figure 7.

Figure 11.

Applications of surface and cross-hole electrical resistivity tomography for mapping base of alluvial overburden: (a) integrated resistivity images from surface and cross-hole ERT and (b) geological section from borehole rock samples and logging data.

6. Conclusions

ERT is a useful near-surface imaging technique, which mainly include data acquisition, numerical modelling and tomographic inversion. For surface data acquisition, dipole-dipole, Schlumberger and gradient arrays are applicable for high-resolution image; particularly, the full-range gradient array may complement to gradient array for better data coverage and completeness of data information. For cross-hole data acquisition, pole-pole (A-M), bipole-pole (AM-N) or pole-bipole (N-AM) and bipole-bipole (AM-BN) can be employed, and the geometry factor and numerical modelling may be applied for designing efficient and effective arrays and exam the imaging capability of ERT for specified targets.

GQG approach is a new version of finite-element modelling. It uses Gaussian abscissae to discrete the model domain and Gaussian weights to compute the volume integral. Therefore, it is much easier to match arbitrary free-surface topography and subsurface interface and computation of the subdomain integrals. It does not require a small size of element and complex element mesh generator. The accuracy of modelling depends on the number of Gaussian abscissae in subdomains. The more abscissae are employed in subdomains, the more accurate modelling result is generated but costs more computer time.

Tomographic inversion is generally implemented by a standard conjugate gradient algorithm, which requires to compute the gradient and the Hessian matrix of an objective function. Two types of objective functions can be applied. One is the l₁-norm and the other is the l₂-norm. The former is less sensitive to an outlier in data but hardly computes the gradient at zero misfit. The latter has no problem to compute the gradient and the Hessian matrix, and it is much easier to implement a standard conjugate gradient algorithm. With high-quality data, the two inversions converge. The field experiments show that surface and cross-hole ERT can be applied to map the base of alluvial overburden and the base of pisolite, as well as the clay contamination within pisolite.

Acknowledgments

The author thanks the Abu Dhabi Education Council for the Award for Research Excellence (AARE17-273) to financially support this work and greatly appreciates the ROBE company for sharing the results of ERT imaging experiments.