1.2 Forward Modelling
A typical DC/IP experiment involves inputting a current to the ground and measuring the potential away from the source. In a time-domain system the current alternates the direction and has off-times between the current pulses at which the IP voltages are measured. A typical time-domain signature is shown in Figure 1. In that Figure,  is the potential that is measured in the absence of chargeability effects. This is the "instantaneous" value of the potential measured when the current is turned on. In mathematical terms this potential is related to the electrical conductivity by
 (1)
where the forward mapping operator  dc is define by the equation
 (2)
and appropriate boundary conditions. In equation (2)  is the electrical conductivity in Siemens/metre (S/m),  is the gradient operator, I is the strength of the input current in Amperes, and rs is the location of the current source. For typical earth structures , while positive, can vary over many orders of magnitude. The potential in equation (2) is the potential due to a single current. This is the value that would be measured in a pole-pole experiment. If potentials from pole-dipole or dipole-dipole surveys are to be generated then they can be obtained by using equation (2) and the principle of superposition.
Figure 1. Definition of the three potentials associated with DC/IP experiments.
When the earth material is chargeable the measured voltage will change with time and reach a limit value which is denoted by in Figure 1. There are a multitude of microscopic polarization phenomena which collaborate so that this final value is achieved but all of these effects can be consolidated into a single macroscopic parameter called "chargeability". We denote chargeability by the symbol  . Chargeability is dimensionless, positive, and confined to the region [0,1).
To carry out forward modelling to compute  we adopt Siegel's (1959) formulation which says that the effect of a chargeable ground is modelled by using the DC resistivity forward mapping but with the conductivity replaced by  . Thus
 (3)
or
 (4)
The IP datum can be either the secondary potential or the apparent chargeability. The former is the difference of the forward modelled potentials with, and without, the IP effect:
 (5)
The apparent chargeability is given by the ratio,
 (6)
In this definition the apparent chargeability is dimensionless and, in the case of data acquired over an earth having constant chargeability 0, we have a=0.
Equations (5) and (6) show that the IP data can be computed by carrying out two DC resistivity forward modellings with conductivities and  . The secondary potential is the more general form of IP data and the apparent chargeability is only define when the linear (or polar) arrays are used along a line on the surface or in the same borehole. When the current and potential dipole-electrodes are arranged in 3-D space so they are not aligned, the total potential can take on positive, zero, or negative values. The cross-line experiments on the surface and cross-hole experiment on boreholes are examples of such situations. Because of the zero-crossing in the total potentials, the commonly used apparent chargeability is undefined. In these cases, the appropriate data to measure the IP effect is the secondary potential. Therefore, we will use secondary potential as the basic IP datum except in the case of linear arrays.
The field data from a DC/IP survey are a set of N potentials (ideally but usually ) and a set of N secondary potentials  or a quantity that is related to . The goal of the inversion is to use these data to acquire quantitative information about the distribution of the two physical parameters of interest: conductivity (x,y,z) and chargeability (x,y,z).
The distribution of conductivity and chargeability in the earth can be extremely complicated. Both quantities vary as functions of position in 3-D space. In addition, there is often large topographic relief. In this program library, the 3-D nature of the physical properties and surface topography are fully incorporated. The earth model is divided into cuboidal cells each having a constant value of conductivity and chargeability. The surface topography is approximated by a piecewise constant surface.
Forward modelling example
The forward modelling for the DC potentials and IP apparent chargeabilities is accomplished using a finite volume method (Dey and Morrison, 1979) and a pre-conditioned conjugate gradient technique to solve Equation (2). The program that performs these calculations is DCIPF3D. The DC modelling is performed by a single solution of Equation (2), and the IP modelling is performed by carrying out two DC forward modellings. The IP data are generated according to the operations indicated in Equations (5) or (6).
To illustrate the DC resistivity and IP forward modelling algorithm, we generate synthetic data that would be acquired over the 3-D conductivity structure shown in Fig 2. The model consists of five rectangular blocks buried in a uniform halfspace. Three smaller blocks are placed at the surface while two larger blocks are at depth to simulate the target of the survey. The blocks S1, S2, and B2 are more conductive than the uniform halfspace; and blocks S3 and B1 are more resistive. All five blocks are chargeable. We have placed 7 east-west lines spaced 100–m apart on the surface and there are four vertical boreholes. The surface experiment is carried out using pole-dipole data with a=50m and n=1, 6, while the borehole experiment uses cross-hole pole-dipole configuration with a 50m potential dipole.
| Figure 2. The synthetic model consists of five rectangular blocks in a uniform halfspace. The blocks S1, S2, and B2 are more conductive than the uniform halfspace; and blocks S3 and B1 are more resistive. All five blocks are chargeable. Also shown are the seven east-west lines (spaced 100 m apart), and the four boreholes that extend to a depth of 400 m. |
|
Figure 3 displays the DC resistivity data from three selected lines for the surface experiment. The data are displayed in pseudo-section format. Note the strong responses to the conductivity anomalies on the surface. They appear as pant-legs extending from small n-spacing all the way to the largest n-spacing. The buried blocks are hardly identifiable since their responses have low amplitudes and broad distributions which are masked by the surface anomalies. The apparent chargeability pseudo-sections from the same lines are shown in Figure 4 (note that the apparent chargeability is well-define in this case). The masking effect of the surface blocks are also evident in the IP data. Thus inversion is required.
| Figure 3. Examples of the apparent conductivity pseudo-sections along three east-west traverses. The data are simulated for a pole-dipole array with a=50 m and n=1 to 6 and they have been contaminated by independent Gaussian noise with a standard deviation equal to 2% of the accurate datum magnitude. The pseudo-sections are dominated by the surface responses, but there are some indications of the buried prisms. The grayscale shows the apparent conductivity in mS/m. |
|
| Figure 4. Examples of the apparent chargeability pseudo-sections along three east-west traverses. The data have been contaminated by independent Gaussian noise with a standard deviation equal to 2% of the accurate datum magnitude plus a minimum of 0.001. The same masking effect of near-surface prisms observed in apparent conductivity pseudo-sections is also present here. The grayscale shows the apparent chargeability multiplied by 100. |
 |
Since we intend to invert these synthetic data, we have added independent Gaussian noise. The standard deviation of the noise is equal to 2% of the datum magnitude plus a small threshold to deal with near zero data. The effect of the added noise can be seen in Figures 3 and 4.
1.3 General Methodology for Inverting DC and IP Data
The computer programs outlined in this manual solve two inverse problems. In the first we invert the DC potentials (or equivalently the data as illustrated in Figure 3) to recover the electrical conductivity (x,y,z). This is a nonlinear inverse problem that requires linearization of the data equations and subsequent iteration. Next we invert the IP data (such as those in Fig ure 4) to recover the chargeability (x,y,z). Because chargeabilities are usually small quantities ( < 0.3 is usual) it is possible to linearize equations (5) and (6) and derive a linear system of equations to be solved. Irrespective of which data set is being inverted however, we basically use the same methodology to carry out the inversions.
To outline our methodology it is convenient to introduce a single notation for the "data" and for the "model". We let d = (d1,d2,...dN) denote the data, where N is the number of data. So di could be the ith potential in a DC resistivity data set, or the secondary potential or apparent chargeability in an IP survey. Let the physical property of interest be denoted by the symbol m. The quantity mi can denote the conductivity or chargeability for the ith cell. For the inversion we choose mi = lni when inverting for conductivities and mi = i when reconstructing the chargeability distribution.
The goal of the inversion is to recover a model vector m = (m1 ,m2 ,...,mM ), where M is the number of model cells that acceptably reproduces the N observations dobs = (d1obs, d2obs,...,dNobs). Importantly, the data are noise contaminated so we don't want to fit them precisely. To do so would ensure that we do not have the correct earth model because some features observed in the constructed model would assuredly be artifacts of the noise. Alternatively, if we fit the data too poorly then information about the conductivity that is coded in the data will not have been recovered. Our objective therefore is neither to underfit nor overfit the data. Rather, we want to find a model that reproduces the data only to within an amount that is justified by the estimated uncertainty in the data. To accomplish this we introduce a global misfit criterion
 (7)
where Wd is a datum weighting matrix. In this work we shall assume that the noise contaminating the jth observation is an uncorrelated Gaussian random variable having zero mean and standard deviation  j. As such, an appropriate form for the matrix is Wd = diag{1/1,...,1/N}. With this choice,  d is the random variable distributed as chi-squared with N degrees of freedom. Its expected value is approximately equal to N and accordingly, d*, the target misfit for the inversion, should be about this value.
Earth conductivity (and chargeability) distributions are complicated. To allow maximum flexibility to produce a model of arbitrary shape it is important that M, the number of cells representing the model, is large. In our inversions M will almost always be greater than N, the number of data. The inverse problem therefore reduces to finding a set of M parameters using only N data constraints under the condition that M>N. Therefore the solution is nonunique and this nonuniqueness represents the principle obstacle for obtaining unambiguous information about earth structure from the observations.
Any inversion algorithm (if it works) will produce a model which reproduces the data. But there are infinitely many possible models. So which one does the algorithm produce? It is not good practise to let the program make a random selection. Rather, a responsible approach is to direct the inversion algorithm to produce a model that is geologically reasonable and is constrained by additional information if that information is available. This can be implemented by formulating a "model objective function" which, when minimized, produces a model with desirable characteristics. The critical aspect of the inversion is therefore to form the model objective function which we characterize by m . To do this, the inversionist must ask the question "what type of model is desired?". Should the model be smooth, should it be blocky? Is there a reference or background model that the constructed model should emulate? If there is such a reference model, is it better known in some places than others? In other words, should the constructed model be close to the reference model in certain locations while being allowed to depart from our preconceived ideas in other areas? Whatever the answer to these questions, a guiding philosophy should always be to find a model which (in some sense) is "as simple as possible". One important guiding principle is that the nonuniqueness inherent in the inversion generally means we can generate models which are arbitrarily complicated, yet we cannot make models that are arbitrarily simple. For example a halfspace will generally not reproduce data acquired from a geophysical survey.
In the inversion algorithms in DCIP3D our choice for the objective function m is guided by a desire to find a model which has (i) minimum structure in the vertical and horizontal directions and at the same time is (ii) close to a reference model m0. To accomplish this we minimize a discretized approximation to
 (8)
In equation (8) the functions are specified by the user and the constant  s controls the importance of closeness of the constructed model to the reference model m0. The constants x, y, z control the roughness of the model in the three directions. We can de\fn a length scale in each direction as  The greater the length scale
is in each direction, the smoother is the model. Varying these scales allows the construction of models that are smoother, thus more elongated, in either x-, y-, or z-direction. To obtain a reasonably smooth model, the length scale should be no less than two cell widths. Given that we always work with a finite model domain, the length scales should be smaller than the respective dimension of the model region.
The discrete form of equation (8) is
 (9)
The matrices Ws, Wx, Wy, Wz, are formed by finite difference approximation of the integrals in eq.(8).
The inverse problem is now properly formulated as an optimization problem:
 (10)
Appropriate techniques can be employed to carry out the minimization and the minimizer m yields the model we are seeking. For DC resistivity and IP inversions we use different minimization techniques that will be discussed in the following sections. In addition, we also impose positivity in the IP inversion to ensure that the recovered chargeability is positive.
References
Dey, A. and Morrison, H. F., 1979. Resistivity modelling structures for arbitrarily shaped three-dimensional, Geophysics, 44, 753–780.
Li, Y., and Oldenburg, D. W., 1994, 3–D inversion of inverse, DC resistivity data usingapproximate, Geophys. J. Int. 109, 343–362.
Li, Y. and Oldenburg, D. W., 1997, Fast inversion of large scale magnetic data using wavelet transform, JACI Annual Report.
Li, Y. and Oldenburg, D. W., 2000, 3-D inversion of induced polarization data, Geophysics, 65 , #6, pp 1931-1945.
Oldenburg, D. W. and Li, Y, 1996, Inversion of induced polarization data, Geophysics, 59, 1327–1341.
Oldenburg, D. W. and Ellis, R. G, 1991, Inversion of inverse mapping, Geophys. J. Int. 105, 325–353.
Siegel, H. O., 1959. Mathematical formulation and geophysical data using an type curves for inducedapproximate polarization, Geophysics, 24, 547–565.
|
