DCIP3D manual: 
1.4 Inversion of DC Resistivity Data

The program library DCIP3D provides two DC resistivity inversion programs. The first, DCINV3D, is based upon a Gauss-Newton method of minimization that requires the linearization of the data equations and the minimization is carried out iteratively. This program can be used to perform rigorous inversion of the DC resistivity data if the interpretation based on the recovered conductivity model is required or if the inversion of IP data requires the best possible background conductivity for calculating the sensitivity matrix. The second program, DCAIM3D, uses a Born approximate inversion and an AIM-MS update to perform the inversion iteratively. This can be used for a first order interpretation of the DC resistivity data, and more importantly, it can be used to inexpensively generate a good approximation to the background conductivity model for IP inversion. We briefly describe both approaches below.

Gauss-Newton method

The inversion of DC resistivity data formulated as the minimization in eq.(10) is nonlinear since the data do not depend linearly upon the conductivity model. We tackle this problem using a Gauss-Newton approach in which the objective function is linearized about a current model, m⁽ⁿ⁾, and a model perturbation is solved for and used to update the current model. Substituting m⁽ⁿ⁺¹⁾ = m⁽ⁿ⁾+m into the objective function in eq.(10),

     (11)

where J is the sensitivity matrix and the element J_ij quantifies the influence of the model change in jth cell on the ith datum,

     (12)

Neglecting the higher order terms and setting to zero the derivative with respect to m yields,

     (13)

Here we assume that the matrix W_d has been absorbed into the sensitivity matrix and data vectors. This is the basic equation that is solved to obtain the model perturbation. The new model is then generated by:

      (14)

where   (0,1] limits the stepsize and is chosen to ensure that the total objective function is reduced.

The major computational effort in this approach includes the calculation of the sensitivity matrix, solution of the basic linearized equation (13), and the choice of regularization parameter . The sensitivity is computed using the standard adjoint equation approach. The equation (13) is solved using a pre-conditioned conjugate gradient technique, in which the sensitivity matrix J is applied to vectors by sparse multiplications in the wavelet domain after it is compressed using fast wavelet transform. Since the wavelet compression is used throughout the DCIP3D Library, we will discuss it in a later section. In the following, we discuss the approaches for choosing the regularization parameter .

Choice of regularization parameter for Gauss-Newton inversion

The choice of the regularization parameter in the DC resistivity inversion ultimately depends upon the magnitude of the error associated with the data. The inversion of noisier data requires heavier regularization, thus a greater value of is required. Since the inversion of DC resistivity data is nonlinear, we also need to avoid the possibility of getting trapped in some local minimum. We have developed two strategies to accomplish this in the program DCINV3D. We discuss their implementations in this section.

If the standard deviation associated with each datum is known, then the data misfit defined by eq.(7) has a known expected value _d^*, which is equal to the number of data when the errors are assumed to be independent Gaussian noise with zero mean. The value of should be such that the expected misfit is achieved. This entails a line search based on the misfit curve as a function of . Because of the computational expense associated with the inversion, we generally cannot afford to perform the line search by carrying out complete solutions for a series of 's. Instead, we start with a sufficientl large value of that reduces the initial misfit by a small fraction. We then reduce by a fixed factor and perform one or two Gauss-Newton updates for each value in the decreasing sequence. Since the sequence starts with a large that leads to stronger regularization, the Gauss-Newton steps bring the model close to the final solution for that . Since the regularization parameter is decreased slowly, it is reasonable to assume that the subsequent models produced by performing one or two Gauss-Newton iterations at the reduced value of would also be close to the corresponding solution. This sequence helps determine the range of the optimum value of . Once the range is established, full minimizations are performed to convergence for a few different values of to determine the optimal value. For each value of in this final stage, the conductivity model from a nearby value is used as initial model so the computational expense is greatly reduced. This strategy is implemented in DCINV3D as the first method for choosing the regularization parameter. The user needs to specify the target misfit value.

In practical applications, the estimate of data error is often not available. Then the degree of regularization, hence the optimal value of , needs to be determined based on other criteria. The L-curve criterion is often used in linear inverse problems to estimate an optimum regularization parameter. It is a heuristic method based upon the relation between the model objective function and data misfit. The underlying premise is the following observation of the Tikhonov curve for the linear inverse problem. As decreases, the data misfit decreases and the model objective function increases. Initially the rapid decrease in misfit is accompanied by a slow increase in model objective function. Beyond a certain point, however, minor reduction in data misfit is achieved only by a large increase in the model objective function. This indicates that the inversion is having difficultie in fitting the noisier features in the data and excessive structure is introduced. When plotted on a log-log scale, the Tikhonov curve appears L-shaped and has a sharp corner. The corner point represents the onset of the rapid increase in model objective function, and it therefore represents the point where the optimum misfit is achieved. We have adopted this criterion in the nonlinear inversion of DC resistivity data. Each iteration of the inversion solves a linearized inverse problem given in eq.(13), and the "error' in the linearized data, d = d⁽ⁿ⁾-d^obs, encapsulates both the actual error in the original data and the error due to the linearization.

Figure 5 Flow chart illustrating the solution of the 3D DC resistivity inversion by DCINV3D using different strategies for choosing the regularization parameter.

Applying the L-curve criterion to the linearized problem at each iteration is expected to yield an optimum value of the regularization parameter for that particular iteration, and the resulting model perturbation m is limited in structural complexity and gives a good direction to update the model. As the solution approaches convergence, the linearization error becomes small and less influential, and the estimated should converge to the optimum value that is appropriate for the actual noise in the data. This approach works well when the data have a reasonably good coverage over the model area. In which case the data exert sufficient constraints on the model so that the noise in the data is not easily reproduced by a simple structure. However, the L-curve criterion can greatly under-estimate the noise when very sparsely distributed data, such as those collected along one or several widely separated lines, are inverted. Nevertheless, it offers an effective approach for determining the degree of regularization in most cases for which the DCIP3D library is intended. Therefore, we have implemented it as one of three methods for choosing the regularization parameter in DCINV3D.

Figure 5 illustrates the structure of the program DCINV3D. It has three options for determining the regularization parameters. The controlling parameter is mode. When mode=1, the
line search based on a known target value of data misfit is used. The regularization parameter is reduced gradually and several solutions for different values of are tested to obtain one that produces the target misfit. When mode=2, the user specifies a regularization parameter and a single solution is produced. When mode=3, the program uses the L-curve criterion to estimate the optimum regularization parameter.

AIM inversion

Given an approximate method by which to invert DC resistivity data, it is possible to iteratively update the conductivity model by using the approximate inverse mapping (AIM) formalism (Oldenburg and Ellis, 1991) such that a final model reproducing the 3D observations is constructed. The AIM inversion requires only one forward modelling and one linear inverse at each iteration and therefore is computationally very efficient. Furthermore, we have observed that the most misfit reduction is achieved within the first two to three iterations (Li and Oldenburg, 1994). Thus, by performing only a limited number of AIM updates, we can obtain a good conductivity model for use in preliminary interpretation or in the IP inversion.

Our approximate inverse mapping is a linearized inversion which uses sensitivity from a conductivity. This is identical to carrying out the first iteration of the Gauss-Newton inversion (see preceeding section) and the equation to be solved is

     (15)

where m^b is the background conductivity, d^b is the predicted data, and d is the data to be inverted. The sensitivity J is calculated using m^b and is stored for use throughout the inversion (which can also be compressed using the wavelet transform for numerical efficiency). Note that unlike Gauss-Newton iterations, the sensitivity matrix is not updated. The regularization parameter is chosen according to the errors in the observed data by using a L-curve criterion. The approximate model is given by =m^b+m. This can be symbolically written as the approximate inverse mapping ^-1[d]=.

The AIM iterations proceed as follows. Let d^obs be the observed data and let m⁽ⁿ⁾ and d⁽ⁿ⁾ be the model and predicted data at the n'th iteration. The application of the approximate inverse mapping to the observed and predicted data yields respectively

     (16)

The updated model m⁽ⁿ⁺¹⁾ is given by

     (17)

The iteration starts with an initial model which can be supplied by ^obs.

Figure 6 The conductivity model recovered from inversion of surface data using a Gauss-Newton method. The model is shown in one cross-section and two plan-sections. The positions of the true prisms are indicated by the white lines.

The AIM algorithm implemented in DCAIM3D uses either a user-supplied model or, by default, the best-fitting constant model as m^b. The user can allow the program to iterate to convergence or limit the number of iterations so an intermediate approximation is obtained. For the purpose of generating a conductivity model for IP inversion, three to five iterations are usually sufficient.

Example: Inversion of synthetic DC Resistivity data generated earlier

The inversion of the apparent conductivity data including those shown in Figure 3 are carried out using the program DCINV3D. The data are collected using pole-dipole arrays with a=50m and n=1 to 6. A total of 1,386 observations are inverted. The same mesh for forward modelling is used in the inversion, and the conductivity value in 15,548 cells is found so that the potential data are adequately fit. Since we added the noise to the data, we use the known values of the standard deviations in the inversion. The appropriate value for the expected data misfit is _d^*. Thus we have run the program using mode=1.

The model objective function is that given in eq.(8). The reference model has a constant conductivity of 1mS/m, and the coefficients are set to _s=0.0001, and _x=_y=_z=1.0. This corresponds to a length scale of 100m in all three directions. The initial model is identical to the reference model. We have also used Daubechies-4 wavelet to perform the compression of sensitivity matrix with a relative threshold of 0.002. An average compression ratio of about 20 was achieved, which leads to the same factor of reduction in CPU time required to solve the linearized equations.

The inversion took 10 iterations to converge. The recovered model is shown in two plan-sections and one cross-section in Figure 6. The recovered model compares favourably with the true model. The surface blocks are well resolved, and there is little excessive structure away from the these blocks. The buried conductive prism is clearly visible and is placed at the correct location. There is also indication of the presence of the buried resistive prism, but it appears to be smoothed out much more. The decrease in resolution at depth is expected since we only have data of limited array separation from the surface.