Geosoft has added induced polarization (IP) and resistivity data inversion to its VOXI Earth Modelling 3D inversion software service. The new toolset is available as part of Geosoft’s 9.1 software update released this month. ERTLab™ is the 3D resistivity and chargeability inversion software that has radically changed the approach to electrical resistivity tomography (ERT) surveys, allowing full flexibility in the three-dimensional arrangement of the electrodes and creating the market of resistivity surveys around buildings.
- Resistivity Inversion Software
- 3d Resistivity Inversion Software Updates
- 3d Resistivity Inversion Software Updates
- 3d Resistivity Inversion Software Update Software
- 3d Resistivity Inversion Software Update
- Software Update on Finite Element Simulation of. 3D Resistivity hp-Finite Element Software — α version— (D. Inversion of Borehole Measurements.
- RES3DINVx64 - 3D RESISTIVITY & IP INVERSION software for Windows XP/Vista/7/8/10 Uses the smoothness constrained Gauss-Newton least-squares optimization method Supports on land and underwater surveys.
Summary
We present a novel technique for the determination of resistivity structures associated with arbitrary surface topography. The approach represents a triple-grid inversion technique that is based on unstructured tetrahedral meshes and finite-element forward calculation. The three grids are characterized as follows: A relatively coarse parameter grid defines the elements whose resistivities are to be determined. On the secondary field grid the forward calculations in each inversion step are carried out using a secondary potential (SP) approach. The primary fields are provided by a one-time simulation on the highly refined primary field grid at the beginning of the inversion process.
We use a Gauss—Newton method with inexact line search to fit the data within error bounds. A global regularization scheme using special smoothness constraints is applied. The regularization parameter compromising data misfit and model roughness is determined by an L-curve method and finally evaluated by the discrepancy principle. To solve the inverse subproblem efficiently, a least-squares solver is presented.
We apply our technique to synthetic data from a burial mound to demonstrate its effectiveness. A resolution-dependent parametrization helps to keep the inverse problem small to cope with memory limitations of today's standard PCs. Furthermore, the SP calculation reduces the computation time significantly. This is a crucial issue since the forward calculation is generally very time consuming. Thus, the approach can be applied to large-scale 3-D problems as encountered in practice, which is finally proved on field data.
As a by-product of the primary potential calculation we obtain a quantification of the topography effect and the corresponding geometric factors. The latter are used for calculation of apparent resistivities to prevent the reconstruction process from topography induced artefacts.
electrical resistivity, finite-element method, inversion, tomography, topography
1 Introduction
The direct current (dc) resistivity method gains information about subsurface conductivity structures by injecting electric currents into the ground and measuring electric voltages at different locations. Since the measured data, however, do not provide the desired information directly, inversion techniques have to be formulated to reconstruct the spatial distribution of conductivity. Although uniqueness has been proved (Druskin 1998), the inverse dc resistivity problem is generally ill-posed with respect to data errors and incomplete data sets.
While initially resistivity soundings used to be applied to determine horizontal layers of varying conductivity, nowadays profile measurements using pre-installed multielectrode lines are state of the art to investigate 2-D structures. However, many objects of investigation exhibit 3-D features, which limit a 2-D interpretation.
As a result of the development of modern field equipment it is now possible to obtain several thousands of data per day. Thus, a target area can be covered by large electrode arrays or a large number of parallel profile lines to make feasible 3-D investigations. Fortunately, computer facilities and numerical techniques are developing rapidly, too, allowing for the development and application of 3-D simulation and inversion software.
3-D inversion of resistivity data is non-linear and usually solved in an iterative process that applies a forward modelling routine for nearly arbitrary resistivity distributions in every inversion step. The forward operator is generally obtained by finite-difference (FD) or finite-element (FE) methods. Since this paper focuses on the inverse procedure we have discussed the forward calculation more rigorously in Rücker (2006, this issue). There we concentrate on accuracy and refinement strategies of a FE formulation operating on unstructured tetrahedral meshes and investigate the solution of the system of equations. Multifrontal direct equation solvers are proposed in conjunction with reordering strategies for the stiffness matrix. The calculation of the secondary potential (SP) may be carried out on coarse grids, which is particularly advantageous for the inversion. However, the determination of the reference primary potential is necessary once. This is obtained most efficiently on locally refined unstructured meshes using quadratic shape functions, which represents the initial point of the inverse approach described here.
The first 3-D inversion scheme was published by Park & Van (1991) followed by Ellis & Oldenburg (1994), Zhang (1995) and Loke & Barker (1996). Most of them are based on Gauss—Newton techniques, some avoid the storage of the Jacobian matrix (Zhang 1995). In contrast, Ellis & Oldenburg (1994) introduced a non-linear conjugate gradient method, see also Zhdanov & Keller (1994). In all these approaches the model parameters are represented by orthogonal grids and the forward calculation is carried out using FDs.
The first inversion scheme using a FE forward calculation was presented by Sasaki (1989, 1994). However, like La Brecque (1995), Yi (2001) and Pain (2003), they use structured grids of rectangular block orientation even if blocks are decomposed into tetrahedrons. There are very few references to unstructured tetrahedral grids (Sugimoto 1999). Zhdanov (2002) gives an overview on existing approaches.
Here we present a full 3-D dc resistivity inversion scheme that involves three different unstructured tetrahedral grids and is able to deal with arbitrary surface topography. The forward calculation is done in two steps, the calculation of the primary potential is approximated by quadratic shape functions and the secondary field by linear ones according to Rücker (2006). We show that the parametrization can be adopted to the resolution properties of the inverse problem. Thus, the number of model parameters can be reduced significantly to save computing time and memory.
2 Triple-Grid Inversion Approach
2.1 Motivation
One of the main problems to deal with is the grid transformation between the inverse and the forward procedure. On the one hand, the parametrization grid has to be relatively crude to limit the degrees of freedom and the ill-posedness of the inverse problem. On the other hand, the forward calculation requires a very fine grid to provide sufficiently accurate results. In order to save computing time, we want to use the SP technique (Rücker 2006) for which we need the primary potentials. If topography is present, those can only be assessed by numerical modelling on a highly refined grid. However this work has to be done only once in the beginning of the inversion. These considerations lead to the triple-grid technique as presented in the following.
The formulation of boundary conditions for the solution of the partial differential equations requires boundaries at the modelling domain that are generally far from the sources and parameter contrasts. This is often accomplished by adding cells of increasing size toward the boundaries. In contrast to the forward modelling grid, the model parameters are usually represented in a parameter grid restricted to the neighbourhood of the sensors. Its number of cells, that is, the degrees of freedom must not be too large due to limitations in run time and computer memory since the Jacobian matrix is dense. At best, the appropriate element size should be defined by the physical resolution. Resolution measures can be obtained by resolution analysis as pointed out by Friedel (2003). In summary, the grids of parameter representation and forward calculation are connected by a refinement and prolongation operation (Figure 1). The primary field mesh is independent of these but locally refined.
The three meshes of inversion for a 2-D discretization: parameter mesh (left), secondary field mesh (centre) and primary field mesh (right).
The three meshes of inversion for a 2-D discretization: parameter mesh (left), secondary field mesh (centre) and primary field mesh (right).
2.2 The Three Meshes
The three meshes needed in the inversion procedure are constructed with respect to the input parameters. The mesh generation itself is a very complicated procedure. We, therefore, use TetGen (Si 2003), a non-commercial and flexible software for the generation of tetrahedral meshes.
The parameter mesh describes the cells, whose resistivities are to be determined and includes the topography of the target area. The lateral model extension is determined by the electrode layout plus a surrounding region in the dimension of several electrode spacings. To determine the model depth we extend the ideas of Roy & Apparao (1971) and Barker (1989) to multielectrode data. A 1-D-sensitivity function for a homogeneous half-space is calculated and summed up for all data. A cumulative value of 90 per cent depicts a reasonable model depth and can be applied to arbitrary electrode configurations (Günther 2004). As an advantage of unstructured meshes their element sizes can be varied flexibly. An a priori resolution analysis (Friedel 2003) helps to get an idea of resolution radii as a function of space. Thus, the smallest elements appear near the electrodes and have typical edge lengths of 0.5—1 times the electrode distance. The element size increases with depth such that the largest elements are created at depth or near the boundaries and have extensions of several electrode distances.
On the secondary field mesh the forward calculations are carried out for the SP. It is obtained through global refinement of the parameter mesh. To avoid effects of the boundary conditions, a mesh prolongation has to be applied. Additional boundaries in the size of a few times the electrode layout are suggested.
The primary field mesh is used to calculate the primary potentials as needed for the SP procedure. We assume a constant conductivity of σ= 1 S/m and scale the potential with the individual value of σ. To cope with the large potential gradients, the mesh needs to be highly refined near the electrode positions to obtain accurate results. Additionally, we use higher-order polynomial shape functions. This also increases the degrees of freedom but turns out to be more effective than purely spatial refinement (Rücker 2006). Consequently, the primary field mesh usually has a huge number of nodal points but the associated system of equations has to be solved only once at the beginning of the inversion process.
2.3 Inversion Scheme
2.3.1 The objective function
A Gauss—Newton scheme with global regularization is used as described in the following. Let m be the model vector of M single model parameters m= (m1, m2, … , mM)T. The individual mj denote the physical properties of the individual tetrahedrons or a function of it. To ensure positivity of the resistivities ρj we use the logarithms mj= log ρj.
In the same manner the logarithms of the measured apparent resistivities log(ρai) =di are used to build up the data vector of N single data d= (d1, d2, … , dN)T. Each data is associated with an error εi, which is used for weighting. Using an ℓ2-norm of the weighted residual between data and model response f (m), the data functional Φd to be minimized is defined by 1
If the errors are not measured, they have to be estimated appropriately. Following Friedel (2003) we assume the errors to be composed of a percentage error of several (p) per cent and a voltage error δU2
Note that the use of logarithmized apparent resistivities leads to the error estimate εi= log (1 +δρai/ρai) (Friedel 2003).
Since the minimization of Φd represents a highly ill-posed problem, we impose regularization by introducing an additional model functional Φm (Tikhonov & Arsenin 1977). It is weighted by the regularization parameter λ yielding the functional to be minimized 3
Eq. (3) can also be obtained by the minimization of Φm while fitting the data down to the error level Φ*d=N. Thus, the regularization parameter has to be chosen to satisfy Φd=Φ*d.
Φm is a squared norm of a product of a constraint matrix C and the difference between the model m and a reference model m04
From here two ways are possible: If m0 is a model containing a priori information and C is the identity or a diagonal weighting matrix, the model is kept close to m0. The other way is to treat m0 as a constant vector and to use C to control the model characteristics. Since the problem is highly underdetermined and the measurements are usually carried out at the surface, the application of smoothness constraints (Constable 1987; de Groot-Hedlin & Constable 1990) seems to be the method of choice. Both methods may also be combined, for example, if a layered half-space is given as reference model m0 and the inhomogeneities are expected to be smooth.
For regular grids the smoothness matrix C represents a discrete approximation of a partial differential operator of first or second order. Here, special smoothness constraints for unstructured meshes have to be defined taking into account the neighbouring relations. Two tetrahedrons i and j are referred to being neighbours, if they share a common triangle face. For each face f we set Cf,i=−1 and Cf,j=+1, instead of 1 the face area size may be inserted. Thus, C∈RB×M is a sparse matrix with 2M entries, where B is the number of boundaries.
In an iterative scheme the model vector is updated 5
with the line search parameter τk, the superscript denotes the iteration number.The application of the Gauss—Newton scheme on minimizing Φ leads to (Park & Van 1991) 6
where S is the Jacobian matrix.2.3.2 Sensitivity Calculation
The sensitivity or Jacobian matrix S∈RN×M contains the partial derivatives of the model responses fi with respect to the model parameters mj7
There exist several methods to obtain sensitivities or Frechet derivatives (McGillivray & Oldenburg 1990). Sensitivity studies have been extensively carried out to obtain concepts of detection and resolution properties. Spitzer (1998) described sensitivities for various configurations and gave an overview on methods for the numerical calculation.
From the reciprocity theorem Geselowitz (1971) derived an analytic expression for the sensitivity of the impedance Z with respect to a conductivity change δσ8
where Ωi is the anomalous region. uS is the potential caused by the source electrodes and uR is the one obtained from interchanging source and receiver electrodes. The corresponding currents are denoted by IS and IR.The potential gradients are known from the FE forward simulation and constant within each element for linear shape functions. In the following we apply the formulation of Kemna (2000), going back to Sasaki (1989), in which the sensitivity is obtained by summation over the potentials and the entries in the element matrix A(e). For details of the FE simulation and particularly the composition of the matrix A(e) see Rücker (2006).
However, problems arise for elements adjacent to electrode positions. The infinite potential on this node has to be replaced by a finite one to restrict the sensitivity values. Note also that the formulation approximates the potential by piecewise linear functions, which is satisfying for the secondary, but not for the total potential on the moderate secondary field grid.
In the course of iterations, the sensitivity matrix can be recalculated. Alternatively, quasi-Newton methods may be applied that use a rank-one update mechanism (Broyden 1972). Since the difference is negligible we use Broyden's method in the following examples to save computing time.
2.3.3 Solution of the Inverse Subproblem
The system of equations (6) has to be solved in every iteration step. In contrast to ray tomography the matrix S is dense. For large-scale problems it is sometimes useful to neglect small absolute values to represent S as a sparse matrix.
To avoid the storage of the left-hand side matrix and the transpose of the Jacobian we use conjugate gradient techniques for solving eq. (6) approximately. The presented algorithm denoted as CGLSCD is an adaptation of the CGLS algorithm (Press 1988) introducing data weighting and model constraints (Günther 2004). For the solution of it reads in the initialization (x0 - initial guess) and computes in every iteration k+ 1 until the residual vector rk is small enough. In every occurrence is replaced by DS, the vector b is initialized with D(d−f (mk)), δx=mk−m0 and L=CTC. The result is the sought model update Δmk. Note that all scalars and vectors may be overwritten in each iteration. The only exception is to save ||rk||2 before rk+1 is computed. A further improvement may be achieved by preconditioning (Haber 2005).
2.3.4 Choice of the Regularization Parameter
The regularization parameter λ is a trade-off between data fit and model roughness and has, therefore, to be chosen appropriately. Small values lead to highly rough-textured models with good data fit, whereas large values correspond to smooth models with weak data fit (Tikhonov & Arsenin 1977). There exist several methods to compromise data fit and model constraints (Farquharson & Oldenburg 2004). A very simple but not always stable method is referred to as L-curve method (Hansen & O'Leary 1993, see Vogel 1996 for instabilities). By plotting the data misfit Φd against the model roughness Φm for a series of λi a curve is constructed, which has often an L-shaped appearance. Figure 2 shows such an L-curve as arising from the first synthetic example of Section 3.1. A reasonable trade-off between Φd and Φm occurs near the ‘corner’. The location may be determined, for example, by estimating the curvature (Li & Oldenburg 1999).
L-curve (data fit against model roughness) derived from the synthetic example of Section 3.1. The optimum regularization strength may be determined searching the maximum of the curvature.
L-curve (data fit against model roughness) derived from the synthetic example of Section 3.1. The optimum regularization strength may be determined searching the maximum of the curvature.
Since the system is better conditioned for strong regularization, we start solving eq. (6) using a relatively large λ and decrease it successively by a fixed factor until the curvature of the L-curve begins to decrease again. This point indicates the optimum λ. Since solutions Δm corresponding to neighbouring λi yield similar results, the result for one value of λ is used as starting vector x0 for the next smaller value. This reduces the number of iterations in the CGLSCD algorithm drastically (Frommer & Maass 1999). In order to construct the global L-curve, we need the forward responses for all values of λ. To avoid this we investigate the L-curve of the linearized inverse problem.
Note that finally the discrepancy principle (fitting data within error bounds) dictates the regularization strength. However, in non-linear inversion many runs would be necessary to iterate until Φd=Φ*d. On the other hand, data errors are seldom known accurately and have to be estimated. In practice, values of 1 to 5 for the chi-squared misfit χ2=Φd/N show reliable results without overfitting or underfitting the data.
2.3.5 Forward Calculation
In every iteration step k, the model response f (mk) has to be calculated. For details we refer to Rücker (2006). We use the singularity removal technique to approximate the SP us, which is the deviation from a background potential up. In matrix notation it is written as9
where Aσ and Aσp−σ are the stiffness matrices for the present conductivity and the conductivity difference σp−σ.Since the calculation is done on a relatively moderate grid we use the multifrontal Cholesky decomposition (Toledo 2001) with approximate minimum degree reordering (Amestoy 1996) to solve the system of equations. The main computational effort bears the Cholesky decomposition, which, however, is done only once. The subsequent back-substitutions are carried out quite rapidly.
The reciprocity measure as proposed by Rücker (2006) is used to control the accuracy of the forward computation. It can be plotted as pseudo-section to visualize the distribution of errors. By observing the mean and maximum values we decide if the secondary field grid has to be refined to increase the accuracy of the forward calculation. In the current stage the refinement is global throughout the modelling domain.
2.3.6 Line search
In each iteration the step length parameter τk has to be determined to prevent the model from overshooting due to non-linearity. We first compute the forward response of the full step length f(mk+Δmk). To avoid the forward calculation for many values τki, we use a linear interpolation between the old and the new model response which is a poor approximation but sufficient to obtain an appropriate step length (Günther 2004). τk is determined such that is minimized. Then, the model is updated mk+1=mk+τkΔmk and the forward response f (mk+1) is calculated exactly. Thus at most two forward runs are needed per iteration step.
2.4 Inversion Scheme
Figure 3 displays the inversion scheme. The following steps are passed through:
Inversion scheme. The topography is used to build the meshes. The primary potential defines the geometric factors and thus the apparent resistivities. The forward calculation is carried out on the secondary mesh. The frame indicates the central loop of the forward and inverse step until the stopping criterion is reached and a solution is obtained.
Inversion scheme. The topography is used to build the meshes. The primary potential defines the geometric factors and thus the apparent resistivities. The forward calculation is carried out on the secondary mesh. The frame indicates the central loop of the forward and inverse step until the stopping criterion is reached and a solution is obtained.
- Read data & topography, estimate errors
- (ii)Triangulate topographical surface
- Generate parameter mesh by resolution estimates
- (iv)Generate secondary field mesh by prolongation and global refinement
- Generate primary field mesh
- (vi)Calculate primary potentials for σ= 1 and interpolation to the secondary mesh nodes
- Determine geometric factors and apparent resistivities
- (viii)Choose starting model and constraints
- Approximate Jacobian matrix from potential
- (x)Calculate model update vector by solution of eq. (6), use L-curve criterion to determine λ in first step
- Determine step length by inexact line search
- (xii)Update model and calculate forward response
- Repeat last four steps until convergence
The stopping criterion applies if (a) the functional Φ stagnates or (b) the data are fitted within errors, that is, χ2=Φd/N≈ 1. If necessary, the regularization parameter λ may be adjusted and the inversion is continued.
3 A Synthetic Example
3.1 Synthetic Model and Data
The following example is descended from the investigation of a burial mound. It illustrates the flexibility of handling a complicated geometry. The topography is built up by intersection of an ellipsoid with a plane. Its footprint is an ellipse of 20 × 15 m extent, the top of the mound is located 5 m above the ground level. In the following, we show that the steep topography, particularly at the feet of the slopes, affects the dc measurements severely and cannot be neglected in the inversion process.
Figure 4 shows the topography and the electrode layout. 9 profiles are placed in x-direction and 13 profiles in y-direction. The electrode distance is 1 m on all profiles, the profile distance is 2 m. In summary 487 electrodes are used.
Topography and electrode layout of the synthetic example. The two red profiles are regarded in the following.
Topography and electrode layout of the synthetic example. The two red profiles are regarded in the following.
The mound itself is supposed to have a resistivity of 200 Ωm whereas for the underground 100 Ωm are assumed. Additionally, a cavity inside the mound is considered. It has an ellipsoidal shape with half-axes of 1.5, 1.0 and 0.5 m in x, y and z direction, respectively, and is described by Neumann boundary conditions. Its centre is located 2.5 m below the top of the mound. Figure 5 shows the synthetic model and the grid used for calculating the synthetic data. The electrode positions are clearly indicated by the high refinement. According to Rücker (2006) we used quadratic shape functions on a mesh with 634 242 nodes.
Synthetic model of a burial mound of 200 Ωm over a homogeneous half-space of 100 Ωm. An ellipsoidal cavity is located in the centre of the mound.
Synthetic model of a burial mound of 200 Ωm over a homogeneous half-space of 100 Ωm. An ellipsoidal cavity is located in the centre of the mound.
Dipole—dipole configurations are considered on all profiles. To counteract the growing geometric factor and redundancy, the dipole lengths are enlarged with increasing separation factor n. In fact, the dipole length is calculated by Δe= (n+ 3 mod 4)a, where a is the initial dipole length and n= 1, 2, 3…. A total of 3439 single data is simulated. The synthetic data are noisified with Gaussian random values. According to eq. (2) we use a standard deviation of 3 per cent plus a voltage-dependent noise of 100 μV at an assumed current of 100 mA for both noise and weighting. In the following we exemplify the two profiles at x= 0 m and x=−10 m.
3.2 Topography Effect
In Figure 6, the influence of the topography on the response along the example profiles are shown. Note that this is done only for illustration, since the geometric factors are numerically known. We (wrongly) use the flat-earth geometric factors for the calculation of the apparent resistivities ρa (Figure 6a). Both in the profile aside the mound (top) as well as in the central profile (bottom) we see strong anomalies, most of which are due to the topography as we will see right away.
Topography effect for two example profiles x=−10 m (top) and x= 0 m (bottom) of the synthetic model from Figure 5: ρa based on flat-earth geometric factors (a), topography effect (b) and ρa using simulated geometric factors (c).
Topography effect for two example profiles x=−10 m (top) and x= 0 m (bottom) of the synthetic model from Figure 5: ρa based on flat-earth geometric factors (a), topography effect (b) and ρa using simulated geometric factors (c).
To eliminate the topographical effects, the real geometric factors were derived from the primary potentials. Referring to Rücker (2006) we define the ratio of flat-earth and simulated geometric factor as topography effect (Figure 6b). It is very strong and produces an increase as well as a decrease in the apparent resistivity of a factor of 3 to 4. The mound shows up with increased apparent resistivities at the central part of the profile, whereas decreased values are indicated at the slopes. Compared to the raw data in (Figure 6a) we see nearly identical structures. Hence, the data are dominated by the topography effect.
Using the simulated geometric factor the data can then be transformed such that the topographical effect is removed (Figure 6c). The remaining anomalies can be clearly associated with the parameter structure. In the central part of the central profile (Figure 6c bottom) the pseudo-section shows the resistivity of the mound (200 Ωm), whereas the data at the profile ends represent the background resistivity of 100 Ωm. Additionally, the high-resistive cavity can be clearly observed in the data by isolated high resistivities in the central part.
3.3 Inversion
We use an automated mesh generation merely based on the topographical information. The grids created in the following are independent of the one used for producing the synthetic data set in Section 3.1. The parameter mesh contains 23 109 parameter cells whose resistivities are to be determined. The secondary field mesh obtains 112 040 nodes and the primary field mesh 421 157 nodes so that the first one may be solved directly by Cholesky decomposition and the latter iteratively by conjugate gradients. Figure 7 shows the parameter mesh and the secondary field mesh. The primary field mesh is similar to the one shown in Figure 5 and, therefore, not shown.
Meshes used for the inversion. The parameter mesh (a) is embedded into the secondary field mesh (b) which is globally refined and prolonged.
Meshes used for the inversion. The parameter mesh (a) is embedded into the secondary field mesh (b) which is globally refined and prolonged.
Starting model is a homogeneous resistivity of 200 Ωm. The regularization parameter has been chosen according to the L-curve criterion in the first iteration (see Figure 2) to be λ= 55.
Table 1 shows the inversion facts for each iteration. The chi-squared misfit (weighted data functional per data) was decreased from 117 to 1.1 in 5 steps. By the line search procedure an essential non-linearity is indicated. The main reduction of the misfit is achieved by just one iteration. Due to the global regularization technique the maximum model parameter range (columns 3 and 4 of Table 1) does not increase in the further course of inversion. Finally, the inversion stagnates at the expected level of 1, that is, the data are fitted within the (known) noise level.
Line search parameter τ, model range and chi-squared data fit in the course of iterations.
Line search parameter τ, model range and chi-squared data fit in the course of iterations.
Figure 8 displays the inversion result. An isosurface with the value of 300 Ωm clearly depicts the location of the cavity. Also, the resistivity contrast between the mound and the underground is clearly visible.
Inversion result of the synthetic model data.
Inversion result of the synthetic model data.
3.4 Comparison with Even Surface
The question arises if the data might be explained without incorporation of topography. Therefore, the locations of the electrodes are projected onto a horizontal surface at z= 0 m (ground level) and the same approach as before is used except that the primary potentials were calculated analytically. Due to the simpler geometry the number of parameter cells reduces to 7448. We inverted the apparent resistivities based on flat-earth geometric factors. However, the inversion process became instable and parameter values exceeded the range of 1–1e5 Ωm without reaching a data fit of less than 20 per cent. We, therefore, abstain from showing these results.
Thereafter we used this parametrization to invert the apparent resistivities based on the simulated geometric factors since they already include the topography effect. Consequently, we were able to fit the data approximately (χ2= 2.5). Figure 9 shows a section of the corresponding inversion result.
Inversion result using the calculated apparent resistivities on even parametrization.
Inversion result using the calculated apparent resistivities on even parametrization.
It shows the main features of the resistive mound in conductive background and also the cavity can be seen. However, the latter anomaly is weaker and smeared over a larger region. Moreover, stronger artefacts at the surface arise. To summarize, although the simulated geometric factors allow for easy inversion methods, the topography needs to be involved into the parametrization to obtain a geometrically reliable model. Note also that the determination of geometric factors by calculating the primary potentials is already the main portion of runtime (see Table 2).
Runtime needed for the inversion of the field data.
Runtime needed for the inversion of the field data.
Assume the topography is sufficiently smooth so that their effects can be neglected. In order to compare with a block-oriented discretization, we created an equidistant grid with dx=dy= 1 m. Using eight layers we obtain a number of 34 × 28 × 8 = 7616 parameters, which is comparable to the unstructured parameter mesh. However, with the latter the model can be described more accurately since the elements are smaller at the surface where the physical resolution is better. Hence, the unstructured meshes seem more efficient even for flat surface, particularly if refinement and prolongation of the forward mesh is considered.
4 Application to Field Data
As a field example, we apply our technique to the investigation of an abandoned mining dump which contains slag material originating from steel production. The deposition has been ceased in 1995. Hardpans formed when the dump material was exposed to precipitation and evaporation thus leading to alteration of the material and the development of silicate gels. dc resistivity measurements have been carried out by the Federal Institute for Geosciences and Natural Resources, Hannover. The aim of the survey was to delineate the hardpan thickness, its function as sealing agent, and, furthermore, to determine internal structures of the dump.
The dump is about 200 m long, 70 m wide and the maximum height is 14 m above the ground level. Figure 10 shows the topography and the electrode layout. The top profile has an electrode separation of 1 m, whereas the 13 cross profiles have electrode separations of 2 m. The profile distance is 10 m, which is five times the electrode separation and, thus, just sufficient to try a 3-D reconstruction of the dump's interior.
Topography and electrode layout of the mining dump. Electrodes were located along 13 cross profiles and one top profile.
Topography and electrode layout of the mining dump. Electrodes were located along 13 cross profiles and one top profile.
In summary 4245 data have been collected with the Wenner array using in total 577 different electrode positions. Figure 11 shows the raw data (a), the topography effect (b) and the apparent resistivities (c) obtained by the simulated geometric factors for the profile at x= 220 m. Some of the main features in the raw data can be recognized in the topography effect which shows deviations from the flat earth case of up to 30 per cent in both directions. For example, the low values of the topography effect along the lower left-hand side flank of the pseudo-section in Figure 11(b) show up as high values in the flat earth apparent resistivities of Figure 11(a). They do not appear in the apparent resistivities with simulated geometric factors (Figure 11c). Although the topography effect is not very pronounced it will account for possible misinterpretations. Hence, we include the topography in the inversion.
Sample profile of the field data in terms of flat earth apparent resistivity data (a), topography effect (b) and real apparent resistivities (c) at x= 220 m.
Sample profile of the field data in terms of flat earth apparent resistivity data (a), topography effect (b) and real apparent resistivities (c) at x= 220 m.
The parameter mesh contains a number of 62 045 tetrahedrons. This is far more than the number of data but it is needed to describe the topography sufficiently. For the primary mesh with quadratic shape functions we obtained 422 980 nodes. Such a large number requires an approximate solution of the system of equations for which we choose the ICCG method (Kershaw 1978; Rücker 2006). The primary field calculation needs about 18 seconds for each electrode on a 2.4 GHz PC such that the whole calculation is done in 3 hr and 18 min including the pre-conditioner.
In contrast, the secondary mesh contains only 111 096 nodes. The resulting system of equations can thus be solved directly. The multifrontal Cholesky decomposition (Rücker 2006) takes about 21 s. One back-substitution for a single source is carried out in less than a second such that the total forward modelling procedure needs about 8 min.
For the error weighting according to eq. (2) we used p= 5% and a voltage uncertainty of δU= 50 μV. By the application of the L-curve criterion we obtained a regularization strength of λ= 22.5. After 7 inversion steps the chi-squared misfit was decreased from 1433 to 4.7 which corresponds to a relative rms of 7 per cent. Table 2 assembles the time needed for the individual inversion parts. The number of 14 forward calculations is due to the inexact line-search procedure in every iteration. Note that for any further run, for example, with changed regularization type or strength, only the last two tasks in Table 2 have to be carried out.
Figure 12 shows the inversion result. Two isosurfaces delineate the conductive dump material (<6 Ωm) and a resistive dump cover (>20 Ωm), possibly due to hardpan development. In regions of missing or thinner hardpan, bulges in the low-resistivity body are noticeable.
Inversion result of the dump data. An isosurface of 20 Ωm is used to delineate the hardpan formation. The dump body is indicated by resistivities below 6 Ωm.
Inversion result of the dump data. An isosurface of 20 Ωm is used to delineate the hardpan formation. The dump body is indicated by resistivities below 6 Ωm.
No obvious internal structure of the dump is revealed which at this stage is interpreted as hint for a homogeneous dump body as far as resistivity is concerned. Laboratory experiments are carried out in order to relate the resistivity to the water content of the material. The resistive dump cover is heterogeneously developed. Measurements to estimate the water retention capability of different hardpan types and thicknesses, as revealed by the resistivity interpretation, are currently ongoing. The interpretation, however, is preliminary since the project is in progress.
With the large amount of data and unknowns the storage of the sensitivity matrix reaches beyond a standard pc's capabilities. Since many of the values are expected to be very small, for example, small cells far from the used four electrodes, we try to save memory by neglecting absolute values below a certain threshold. Consequently, we store S as a sparse matrix. Note that 12 Bytes are used for the indices and the value of each non-zero element in contrast to 8 Byte for a full matrix element.
We decrease the threshold value logarithmically from 10−2 down to 10−6. In order to avoid effects of the forward calculation we compute the difference between the first iteration model (without line search) with and without sensitivity sparsing. Table 3 shows the non-zero values and the model difference as relative rms.
Effect of sparsing the sensitivity matrix on the result of the inverse subproblem for different threshold values, the number of non-zero elements is also given as fraction of the full matrix and required memory, model rms is the relative rms of the first iteration model with and without sparsing.
Effect of sparsing the sensitivity matrix on the result of the inverse subproblem for different threshold values, the number of non-zero elements is also given as fraction of the full matrix and required memory, model rms is the relative rms of the first iteration model with and without sparsing.
We observe a very good accordance from tol = 10−5 or below. Thus, we are able to obtain accurate results with 10 to 20 per cent of the memory requirements, which is supported by comparing the final models for various data sets.
5 Discussion and Conclusions
An inversion strategy for the reconstruction of conductivity from dc measurements for arbitrary topography was presented. It is an enhancement of existing techniques and combines the fast convergence of regularized Gauss—Newton methods with accurate FE forward calculations. Unstructured tetrahedral meshes are used in order to describe the topography of the measurement area with high accuracy. Special attention is paid to an efficient forward simulation using a SP approach. Therefore, we split up the forward process into two parts: the time-intensive calculation of the primary potential and the fast calculation on the moderate secondary field mesh.
We use a Gauss—Newton approach with inexact line search for solving the inverse problem. Thus a solution is guaranteed after a small number of iterations. To stabilize the solution we apply first-order smoothness constraints. However, other regularization methods have already been implemented such as second-order smoothness constraints, the introduction of a priori information as well as weighting mechanisms for different model regions. However, their detailed description goes beyond the scope of this paper. A further enhancement may be obtained by active constraint balancing (Yi 2003). Also, some kind of focused inversion (Portniaguine & Zhdanov 1999) may be introduced to obtain models with sharp boundaries.
The strategy was applied to noisified synthetic data as well as to field data. The run time is strongly decreased by the singularity removal technique and the use of fast direct equation solvers. In addition, the use of unstructured grids provides a resolution-dependent parametrization, which finally saves computing time and memory. Thus, the approach is especially interesting for surveys with steep topography or large geometry contrasts. Also, the approach suits well for complicated model geometries, for example, underwater surveys, model tank experiments or measurements on cylindrical geometries.
In cases of strong conductivity contrasts, particularly near the electrodes, errors arise in the forward calculation that can only be diminished by further refinement. In order to achieve high speed and low memory consumption, the refinement should be carried out adaptively such that only regions of weak approximation are locally refined. Another possibility of improving the forward calculation is a grid hierarchy used for multigrid solvers or pre-conditioners. Since in the presented technique the Jacobian matrix needs to be stored, the memory capacities of standard PCs may represent a limit for many large-scale problems. Thus methods are sought that avoid the explicit storage of the Jacobian matrix. One solution may be the technique of Zhang (1995) or the non-linear conjugate gradients (NLCG) method (Ellis & Oldenburg 1994). Due to the slower convergence of NLCG much more iterations are needed to find the minimum of the objective function. However, since the forward calculation is very fast once the primary potential has been calculated, this technique is particularly promising for large-scale problems. Other memory-saving techniques are quasi-Newton methods as proposed by Haber (2005) where the Jacobian is approximated by a sum of rank-two update matrices each stored by two vectors. Thus the method starts with the slow convergence of gradient methods and ends up with the quadratic convergence of Newton-type methods.
A further development can be achieved by incorporating induced polarization effects using complex conductivities (Kemna 2000). Since the imaginary part is generally small, the primary potentials may remain real, only the SPs have to be obtained complex. Moreover, the use of electric in-hole, hole-to-surface and cross-hole measurements can restrict the ambiguity of the surface data. Finally, appraisal techniques as presented by Alumbaugh & Newman (2000) or Friedel (2003) may help to find the reliability of the models.
Acknowledgments
The research was partially supported by the German Research Foundation (Ja 590/18-1). We wish to thank U. Noell, M. Furche and C. Grissemann from the Federal Institute of Geosciences and Natural Resources, Hannover, for providing the high-quality field data within a project supported by the Federal Ministry of Education and Research (BMBF, FKZ 0330523).
References
D.L.
G.A.
, . Image appraisal for 2-d and 3-d electromagnetic inversion
, , 65
(), 1455
–. P.R.
T.A.
I.S.
, . An approximate minimum degree ordering algorithm
, , 17
(), 886
–. R.D.
, . Depth of investigation of collinear symmetrical four-electrode arrays
, , 54
(), 1031
–. C.G.
, . Quasi-Newton methods
, in , ed. W.
, , London.
S.C.
R.L.
C.G.
, . Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data
, , 52
, –300
. C.
S.C.
, . Occams inversion: to generate smooth two-dimensional models from magnetotelluric data
, , 55
, –1624
. V.
, . On the uniqueness of inverse problems from incomplete boundary data
, , 58
(), 1591
–. R.G.
D.W.
, . The pole-pole 3-d dc-resistivity inverse problem: a conjugate gradient approach
, , 119
, –194
. C.G.
D.W.
, . A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems
, , 156
, –425
. S.
, . Resolution, stability and efficiency of resistivity tomography estimated from a generalized inverse approach
, , 153
, –316
. A.
Resistivity Inversion Software
P.
, . Fast cg-based methods for Tikhonov-Phillips regularization
, , 20
(), 1831
–. D.B.
, . An application of electrocardiographic lead theory to impedance plethysmography
, , BME-18
(), 38
–. T.
, . Inversion Methods and Resolution Analysis for the 2D/3D Reconstruction of Resistivity Structures from DC Measurements
, , Freiberg University of Mining and Technology, available from http://fridolin.tu-freiberg.deE.
, . Quasi-newton methods for large-scale electromagnetic inversion problems
, , 21
, –323
. P.C.
D.P.
, . The use of the L-curve in the regularization of discrete ill-posed problems
, , 114
, –1503
. A.
, . Tomographic inversion of complex resistivity
, Dissertation, Der andere Verlag, Ruhr-UniversitätBochum.
D.S.
, . The incomplete cholesky-conjugate gradient method for the iterative solution of systems of linear equations
, , 26
, –65
. D.
G.
A.
P.
, . Occam's inversion of 3D ERT data
, in Proceedings I. International Symposium on Three-Dimensional Electromagnetics in Ridgefield, CT, Oct 4–6, 1995
, pp. –477
, . 3d Resistivity Inversion Software Updates
Y.
D.W.
, . 3-d inversion of dc resistivity data using an L-curve criterion
, in 69th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts
, pp. –254
. M.H.
R.D.
, . Practical techniques for 3d resistivity surveys and data inversion
, , 44
, –523
. P.R.
D.W.
, . Methods for calculating frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study
, , 38
(), 499
–. C.C.
J.V.
J.H.
M.H.
C.R.E.
, . Anisotropic resistivity inversion
, , 19
, –1111
. S.K.
G.P.
, . Inversion of pole-pole data for 3-d resistivity structure beneath arrays of electrodes
, , 56
, –960
. O.
M.S.
, . Focusing geophysical inversion images
, , 64
(), 874
–. W.H.
S.A.
W.T.
B.P.
, 1992. Numerical Recipes in C: The art of scientific computing
, , Cambridge
. A.
A.
, . Depth of investigation in direct current methods
, , 36
(), 943
–. C.
T.
K.
, . Three-dimensional modeling and inversion of dc resistivity data incorporating topography — I. Modelling
, , doi: (this issue). Y.
, . Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data
, , 54
, –262
. Y.
, . 3-d resistivity inversion using the finite-element method
, , 59
(), 1839
–. H.
, . TetGen: A 3d delaunay tetrahedral mesh generator
, K.
, . The three-dimensional dc sensitivity for surface and subsurface sources
, , 134
, –746
. Y.
, . Shallow high-resolution 2–d and 3–d electrical crosshole imaging
, , 18
, pp. –1428
. A.N.
V.Y.
, . Solution of Ill-Posed Problems
, , Washington DC
. S.
D.
V.
, . Taucs—a library of sparse linear solvers
. C.R.
, . Non-convergence of the L-curve regularization parameter selection method
, , 12
, –547
. M.-J.
J.-H.
Y.
3d Resistivity Inversion Software Updates
S.-J.
S.-H.
J.-H.
, . Three-dimensional imaging of subsurface structures using resistivity data
, , 49
(), 483
–. M.-J.
J.-H.
S.-H.
, . Enhancing the resolving power of least-squares inversion with active constraint balancing
, , 68
(), 931
–. J.
R.L.
T.R.
, . 3-d resistivity forward modeling and inversion using conjugate gradients
, , 60
(), 3d Resistivity Inversion Software Update Software
1313
3d Resistivity Inversion Software Update
–.M.S.
, . Geophysical Inverse Theory and Regularization Problems
, , Elsevier
, . M.S.
G.V.
, . The Geoelectrical Methods in Geophysical Exploration, Methods in Geochemistry and Geophysics
, , Amsterdam, London, New York, Tokyo
. © 2006 The Authors Journal compilation © 2006 RAS