Direct inversion of the apparent complex-resistivity spectrum

Cole–Cole model parameters are widely used to interpret electrical geophysical methods and are obtained by inverting the induced polarization (IP) spectrum. This paper presents a direct inversion method for parameter estimation based on multifold least-squares estimation. Two algorithms are described that provide optimal parameter estimation in the least-squares sense. Simulations demonstrate that both algorithms can provide direct apparent spectral parameter inversion for complex resistivity data. Moreover, the second algorithm is robust under reasonably high noise.


INTRODUCTION
Spectral induced polarization (IP) is widely used in geological surveys (He et al., 1995;He, 1996;Milson, 1996;Luo and Zhang, 1998).The interpretation of the IP data is often based on the Cole-Cole model (Cole and Cole, 1941;Luo and Zhang, 1998), under which the frequency behaviour of the complex IP impedance (or complex resistivity) is approximated by an equivalent network, as shown in Figure 1 (Major and Silic, 1981).
The impedance of the equivalent network is (1) The complex resistivity expression of the Cole-Cole model (Luo and Zhang, 1998) is written where φ(ω) is complex resistivity, φ 0 is resistivity at zero frequency (ω = 0), m is limited polarizability (or chargeability), φ is a time constant, c is frequency dependence, and ω is angular frequency (rad/s).The chargeability 0 m 1 and the requirement that the complex resistivity amplitude monotoni-Manuscript received by the Editor February 24, 2000; revised manuscript received February 1, 2001.University of Leicester, Department of Engineering, University Road, Leicester LE1 7RH, U.K. E-mail: jx5@le.ac.uk; nbj1@le.ac.uk; fss1@le.ac.uk.‡Institute of Systems Science, Chinese Academy of Sciences, Beijing 100080, P.R. China.E-mail: dcheng@iss03.iss.ac.cn.c 2001 Society of Exploration Geophysicists.All rights reserved.
cally decreases with increasing frequency restricts c to the range 0 c 1 (Major and Silic, 1981).Interpretation based on the Cole-Cole model requires the parameters φ 0 , m, φ, and c to be estimated (Luo and Zhang, 1998).
The existing methods of parameter estimation are generally based on nonlinear iterative inversion (Pelton et al., 1984(Pelton et al., , 1978;;Jaggar and Fell, 1988;Luo and Zhang, 1998).There is a serious shortcoming to this approach: The convergence of an iterative inversion to a global minimum is not ensured because it depends on the initial guess for parameter values.Different starting values may yield different inversion results.
We propose a new direct inversion technique based on multifold least-squares estimation.The basic procedure is described by the following sequence.First, use substitution to eliminate the parameter m from the equations.Second, use least-squares estimation to express the compound parameter X = φ c as a function of φ 0 and c, which is expressed as φ c = X (φ 0 , c).Third, consider the reduced system of equations for φ 0 and c, which is linear with respect to φ 0 .Then the secondfold least-squares estimation can be used to get a solution for φ 0 as a function of c.
The final step is to substitute all the estimated parameters into the equations to get a set of linear equations of X .Both the coefficients and the unknowns, X , of the linear equations are functions of c.A square error can be constructed for the leastsquares solution of X as a function of c.Since the parameter c ∈ (0, 1), the golden section technique is used to search for the minimum solution c. Finally, the least-squares technique is used to estimate m.
This inversion technique is based on a straightforward computation that identifies the four parameters simultaneously and guarantees a unique solution.Moreover, the solution is optimal in the least-squares sense, and the method can be applied easily to the case of multiple, linearly additive Cole-Cole models.

DIRECT INVERSION OF THE COLE-COLE MODEL
Consider inverting the ].That is, from experimental data estimate φ 0 , m, φ, and c.
Assume N + 1 data are obtained as where φ k is complex resistivity at ω = ω k .
To reduce the number of parameters, we divide the kth equation by the (k + 1)th equation for k = 1, 2, . . ., N .After some algebra we get N equations as Taking the reciprocal, equation ( 4) is equivalent to where Note that data with repeated φ k should be deleted to avoid ballooning to infinity.
Setting φ c = X and separating the real and imaginary parts of equation ( 5) yields the following equations: where and

FIG. 1. Equivalent circuit of IP geophysical system
To get the least-squares approximate solution of equation (7), we have to minimize the square error form S: Setting dS/dX = 0, the minimum point of the strictly positive S is obtained as Substituting the expressions for A k , B k , and C k into equation (11), after some simplification, yields where X (c, φ 0 ) emphasizes that the least-squares approximate solution for X depends on both c and φ 0 .
Next, to estimate φ 0 some new notation is convenient.Denote Then, from equation ( 6) the following expressions can be derived: and Substituting equations ( 15) and ( 16) into equation ( 12) provides where and Applying equations ( 17)-( 20) to equation ( 7) produces Equation ( 21) can be expressed in matrix form as where and Here we use M(c) and L(c) to emphasize they depend on c only.Now the least-squares solution of equation ( 22) is Using equation ( 25), the square error of equation ( 7) becomes In practice, we can assume that φ 0 is a real number.That is, φ 0 = R 0 and I 0 = 0. Then equation ( 23) can be replaced by and To determine c, we minimize S(c).That is, Since c ∈ (0, 1), it lies on a restricted range.So we use the golden section method to search for the optimal solution (Press et al., 1986).
After the optimal c (denoted c ) is obtained, equation ( 25) [or equation (28) in the real case] provides the best estimate of φ 0 ; X = φ c can be obtained by equation ( 11).The estimation of φ follows easily.
Finally, equation ( 3) can be used to estimate m.Since m is real and m > 0, equation ( 3) can be rewritten as Now, for this set of real equations the technique used in equations ( 7)-( 11) can be used again for obtaining a least-squares solution as
The sixth and final step in this algorithm is to use equation (26) to calculate S(c 1 ) and S(c 2 ) and to compare them to decide if a convergence criterion for the solution has been satisfied.If |c 1 c 2 | < ε for a given 0 < ε 1 , stop.Otherwise, use the golden section algorithm to choose a new point and set new c 1 , c 2 .Then go back to step 3.
We have written a Matlab program that performs this algorithm.The following example illustrates the accuracy of the approximation.Take the sampling frequencies as that is, from ω = 2 12 to ω = 2 +7 rad/s.Use equation ( 2) to get synthetic data values of φ k = φ( jω k ).Then invert the data set.
The results using the above algorithm are shown in Tables 1  and 2 (where c is the estimated value and E c % is the relative error, etc.).
McInnes (personal communication, 2000) points out that the above algorithm is sensitive to noise in the observed apparent complex resistivity data.Motivated by this observation, we modified the performance criterion S to optimize c.Let S be the real square error of the estimation Algorithm 2 We modified algorithm 1 as follows.Using c, φ 0 (c), andφ(c), we estimated m(c) first.Then we used the estimated parameters and Cole-Cole model to calculate the model response φ c (ω k ).Equation ( 32) can be used to calculate the model response misfit.Finally, a single-variable minimization algorithm can be used to find the best estimation of c.
To perform this refined new algorithm, you repeat the first five steps for algorithm 1.As the sixth step, use equations ( 28), (17), and (31) to find φ 0 , φ, and m, respectively.For the supposed c, the corresponding model is built.As the seventh and final step, use equation ( 32) to calculate S(c 1 ) and S(c 2 ).Compare them to decide if a convergence criterion for the solution has been satisfied.That is, if |c 1 c 2 | < ε for a given 0 < ε 1, stop.Otherwise use the golden section algorithm to choose a new point and set new c 1 , c 2 .Then go back to step 3.
It is likely that the direct Cole-Cole inversion function S(c) will have multiple local minima with some data sets.Then the interval (0, 1] can be divided into several subintervals, and the algorithm can be used over each subinterval to find the real minimum. The algorithm is described by Figure 2, which consists of two essential parts.The first is the multifold, least-squares estimation, which uses the assumed parameter c and the field data Table 1.Estimation for parameters with different accuracy.to estimate φ 0 (c), m(c), and φ(c).Part 2 is the golden section algorithm.
To represent the accuracy of the approximation, we use a square error criterion.That is, let p(n) and q(n), n = 1, 2, . . ., N , be the field data and the theoretical data obtained from the Cole-Cole model, respectively.Then the square error is defined as , which is the same error criterion minimized by Jaggar and Fell (1988) to estimate their model parameters.
Using algorithm 2, we revisit example 1.When there is no measurement noise, the parameters can be recovered precisely without visible error (within Matlab's precision for long variables, that is, ε < 0.000015).When white noise is added to the complex resistivities, the values of estimated parameters will change.The algorithm is said to be robust if this sensitivity of the parameters to measurement noise is small.Table 2 illustrates the robustness of the algorithm.
A relative noise level of 10 2 % represents a 0.1 mrad IP phase error for the first harmonics, which is good-quality data for typical dipole-dipole surveys (S.McInnes, personal communication, 2000).From Table 2 one can see that the robustness of algorithm 2 is adequate for practical use.
The range of frequencies used affects the accuracy of the approximation.In the example above, if we narrow the range of frequencies to N = 12 and ω = 2 k , k = 6, 5, . . ., 5, the range of the frequency, in hertz, is 2.49 10 3 < f < 5.1.Then the related estimation is obtained as in Table 3, which shows the accuracy is still acceptable.Jaggar and Fell (1988) report N = 19 data values, but highfrequency data values have large uncertainty.Table 4 shows the results estimated via algorithm 2 by using the first N = 12, FIG. 2. Flow chart of algorithm 2. 13, . . ., 17 data values.The best-fit model suggested by Jaggar and Fell (1988), data size 19 on Table 4, is compared with the models found with our methods in Table 5.

FIELD DATA INVERSION
Another advantage of our algorithm is that it does not need an initial guess.Table 6 shows that the ridge regression method adopted by Jaggar and Fell (1988) clearly depends on the initial guess.It might only find a local minimum.Table 6 shows three models (N -13, N -16, N -17) provide better estimation than the J-F model, even over all 19 sampling data points.
If we ignore the last two outlier data values, then the model errors over the set of N = 12-17 reliable data values are shown in Table 7.So from both Tables 6 and 7 we can see that model N-16 seems to best fit the field data.We therefore can choose N -16 as our finally accepted model.Figures 3 and 4 compare the field data with the theoretical data.

CONCLUSION
A direct inversion technique for the Cole-Cole model was presented.The method is based on a multifold least-squares estimation combined with an optimal searching technique.Two step-by-step algorithms were derived and described via several numerical examples.The advantages of the final new approach are 1) the estimation is optimal in the least-squares sense, the existence of a solution is guaranteed, 3) the algorithm is straightforward and simple, 4) no iterative operator intervention is required, and 5) amplitude and phase can be fit simultaneously.
Using the field data provided by Jaggar and Fell (1988), our algorithm has been verified and a detailed comparison with their existing algorithm has been completed.The results show the advantages of the new algorithm.

Table 4 . Estimated parameters via field data.
us the lack of robustness in the previous version of the algorithm.That challenged us to improve the technique presented here.