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:

$$\Phi(\mathbf{m}) = \left\| \mathbf{d}\_{ab} - \mathbf{d}\_{\rm sym}(\mathbf{m}) \right\|\_{\mathbf{W}\_d}^{l\_p} + \lambda \left\| \mathbf{m} \right\|\_{\mathbf{W}\_m}^{l\_p} \tag{24}$$

where dob is a vector of observed data, which are either apparent resistivities or potential differences measured by different electrode arrays. dsyn(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 lp in Eq. (24)). k k� <sup>W</sup> stands for the weighted lp-norm with <sup>a</sup> weighting matrix <sup>W</sup>, e.g. <sup>2</sup> <sup>1</sup> k k<sup>ε</sup> <sup>W</sup> <sup>¼</sup> <sup>ε</sup><sup>T</sup>W<sup>ε</sup> and k k<sup>ε</sup> <sup>W</sup> <sup>¼</sup> <sup>∑</sup><sup>i</sup> Wi∣εi∣ are the weighted l2-norm (generalised least

Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

Figure 8. An l2-normed (left) and an l1-normed (right) objective function. Here I is an unit matrix.

square) and the weighted absolute norm. Figure 8 illustrates the difference between a l2-normed and l1-normed objective functions. W<sup>d</sup> and W<sup>m</sup> are weighting matrices to data and model, respectively. It is common to choose a combination of finitedifference operators for <sup>W</sup>m, e.g. <sup>W</sup><sup>m</sup> <sup>¼</sup> <sup>λ</sup>0Iþλ1D<sup>T</sup>Dxþλ2D<sup>T</sup>Dyþλ3DTD<sup>z</sup> [29]. <sup>x</sup> <sup>y</sup> <sup>z</sup> Here I is a unit matrix, and Dx, D<sup>y</sup> and D<sup>z</sup> are the finite-difference operators in the x-, y- and z-directions. λ<sup>k</sup> (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}^\* = \min \{ \Phi(\mathbf{m}) \}. \tag{25}$$

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 ∇Φ(mi) and the Hessian matrix H(mi) of the objective function are required. Applying the linearised approximation to the synthetic data, <sup>d</sup>syn(mi+1) <sup>≈</sup> <sup>d</sup>syn(mi)+(<sup>∂</sup>dsynÞ<sup>T</sup>(mi+1-mi), <sup>∂</sup><sup>m</sup> and then substituting it for Eq. (24), one can obtain the gradient and the Hessian matrix [28]:

#### Figure 9.

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

$$\nabla \Phi(\mathbf{m}\_i) = \left(\frac{\partial \mathbf{d}}{\partial \mathbf{m}\_i}\right)^T \mathbf{W}\_d \left[\mathbf{d}\_{vb} - \mathbf{d}\_{\partial n}(\mathbf{m}\_i)\right] + \lambda \mathbf{W}\_m (\mathbf{m}\_i - \mathbf{m}\_0), \tag{26}$$

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

for the l2-normed objective function.

Note that the gradient ∇Φ(mi) of the l1-normed objective function has multiple values at zero misfit Φ(mi) = 0 (see Figure 8), but for any non-zero misfit Φ(mi) ≥ ε<sup>0</sup> (ε<sup>0</sup> is a very small value), the computations of the gradient ∇Φ(mi) and the Hessian matrix H(mi) become simple, and they are similar to Eqs. (26) and (27) except for the weighting matrices W<sup>d</sup> and Wm. Therefore, as Φ(mi) < ε0, the l1 normed gradient ∇Φ(mi) and the Hessian matrix H(mi) may be replaced with Eqs. (26) and (27), and as Φ(mi) ≥ ε<sup>0</sup> the l1-normed gradient and the Hessian matrix are directly applied to the standard conjugate gradient algorithm. Therefore, the l1-normed inversion and the l2-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 l1-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 l2 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

#### Figure 10.

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

Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

#### 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.

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.
