Inversion with an alternative error function in resistivity measurements

Geoexppforation, 23 (1985) 527-536 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

INVERSION WITH AN ALTERNATIVE RESISTIVITY MEASUREMENTS

J. POUS, A. MARCUELLO Departamento de Fisicade 08028 Barcelona (Spain)

ERROR

FUNCTION

527

IiN

and P. QUERALT la Tierra y de1 Cosmos, Universidad de Barcelona, Diagonal 645,

(Received June 13, 1985; accepted September 18,1985)

ABSTRACT

POW, J., Marcuello, A. and Queralt, P., 1985. Inversion with an alternative error function in resistivity measurements. Geoexploration 23: 527-536. An alternative inversion method is proposed, the error function of which evaluates an average between the slopes of apparent resistivity more accurately than the classical one of points. Convergency is achieved by this procedure in cases of initial models for which classical error functions fail to converge. The generation of equivalent models with respect to the slopes is also discussed.

INTRODUCTION

The numerical interpretation of the vertical electrical soundings has been developed by various authors in geophysical literature and different techniques have been applied (Inman et al., 1973; Jupp and Vozoff, 1975; Bichara and Lakshmanan, 1976; Johansen, 1977; Koefoed, 1979). All these techniques basically consist of iterative procedures in which beginning with an initial model (resistivities and thicknesses) successive corrections to this model are applied until the final model is reached. The way to apply the corrections defines each technique. Due to the equivalences between models, there is a strong dependence of the final model with respect to the initial one. The final model is one among the set of equivalent models and it may not be one corresponding to the geological reality. In practice we can be fortunate if the algorithm is able to find.bne between all the equivalent models which corresponds to a decrease of the error function in each iteration. Once the solution is achieved, the generation of its equivalent models provides the set of geoelectrical possibilities and any additional geological information will permit us to choose the real equivalent model. One of the main problems is the choice of the initial model because if a

0016-‘7142/85/$03.30

0 1985

Elsevier Science Publishers B.V.

528

wrong model is chosen, the procedure may not be convergent. This is specially important in techniques that require a linearization and therefore the initial model needs to be close to the solution model. Experience with these techniques demonstrates that the sign of a parameter correction is not the same in all iterations. Moreover, there are incompatibilities in the sign and size of the different components of the parameter correction vector, that is, although the error function decreases in a given iteration, some parameters may be moved away from their correct final values. In this case and if the initial model is correct, the final model can be achieved within some more iterations. If not, the procedure will be divergent. We think this is mainly due to two reasons: (1) The error function commonly used for testing the goodness of a model is an average between the experimental points of the apparent resistivity and the theoretical ones. This implies that minimization of the error function produces a simple smoothing of the theoretical curve with respect to the experimental one and the interpretation may be completely wrong. In this case the error function is within a local minimum, the value of which is too high. Even if the iterative procedure goes on, the algorithm is not able to give up this local minimum. In other words it is not able to find the parameter correction vector that would allow to give up the region around the local minimum and to reach another minimum of smaller value of the error function. (2) Because of the great number of parameters, the hypersurfaces with same value of the error function will have very complex forms. Therefore, the vector correction in each iteration may be very unstable. On the other hand, the properties of the VES curves on the logarithmic scale are well known. Experience with the auxiliary point method (Oreflana and Mooney, 1966) shows that slope adjustments are largely important. We think that an error function that evaluates the average over the apparent experimental and theoretical slopes of the apparent resistivity is better than the classical one averaging the points. Moreover, as slopes do not depend directly on the resistivities of the layers but on their contrast, such a procedure will give a better flexibility in each iteration because the set of possible solutions is increased. This does not happen with error functions such as the ones proposed by Bichara and Lakshmanan (1976) which are more severely constrained. On the other hand an error function with only slopes decreases the incompatibilities mentioned before and concerning techniques that need a linearization, the region around the point where linearity is valid, is increased. In this paper we present inversion over a stratified earth with a new error function that evaluates an average of the slopes. Examples will refer to the Schlumberger array. The analysis of equivalences with respect to the new error function is also discussed here.

529

FORMULATION

OF THE INVERSE

PROBLEM

WITH THE NEW ERROR

FUNCTION

As a new error function that evaluates the slopes, we will use the one proposed by Pous et al. (1984) in the steepest descent method: Q* = l/(M - 1) Mg’

((1 n Pai+

- ln Pai) - [ln Pa(si+ l,P) - In PatSiS)

1 2 (1)

i=l

where pai is the measured

value of the apparent resistivity at abscissae si, the components of which are the logarithms of the parameters, resistivities and thicknesses, p = (In p1 In pz , . . ., In pi, ln tl, . . . In tN_l), N being the number of layers. pa(si,p) is the theoretical apparent resistivity corresponding to model p in abscissae Si. We look for a p vector that minimizes Q*. We assume that the experimental and theoretical slopes are i = 1, . . . M. p is a vector

??li = In Pai+ 1 - In Pai mI = ln Pa(Si+l,P)

i=l,...M-1

(2)

- In Pa(si,P)

where the abscissae are supposed to follow a geometrical progression of a given ratio. Otherwise, Q* would not represent an error function for the slopes. By expansion of rni in a Taylor series in the parameter space around the point p” and discarding all terms higher than the first order we get: 2N-1 mi z

mi”

+

c

(amilapj)po

* AP~

j=l

where mi” corresponds top’ and Ap = p - p”. In matrix notation eq. 1 will be:

Q* = Ilb - AxI1 = (b - ii~)~(b - Ax) where vectors

b and x have the components:

bi = mi _ mi” Xi = Api and the matrix 7i has the elements: Aij= am{/apj

i = 1, . . . . M - 1

j=1,...,2N-1 It is clear that: Ai,j = Ai+ 12 - Ai,j

(3)

530

where Ai,j are the elements: Ai,i = a ln P,(si,P)/aPj Matrix A is the one involved in the classical error function the apparent resistivity points. Minimization of Q* is achieved by the generalized inverse:

by averaging

--

of Lanczos where A+ = Vq A,-- 1 $ from the singular value decomposition (1961). _& is a diagonal matrix with the r~ non-zero singular values of x. Matrix U, is formed with the q orthonormalised eigenvectors of dimension (M - 1) associated with the column vectors of x and the matrix v, of q orthonormalised eigenvectors of dimension (2N - 1) associated with the row vector of A. The new parameter model will be: p=po

+;

After various iterations the solution is achieved. Then, the corresponding apparent resistivity curve will be vertically displaced (in the bilogarithmic representation) from the experimental one. This displacement can be found by calculating the vertical distance between the experimental points and the theoretical ones which correspond to model j. Therefore the displacement will be: M

K = l/M

C

P,ilP,(SivP)

(4)

i=l

thus if j = (In pl, . . . In PN, In tl,. . . In tN-1 ), the final shifted model will be: p* = (In Kp,,

. . . In K/IN, In tl, . . . In to-

1)

As it is known the solution i given by A’ may lie far outside the environment where linearity is valid. This is due to the existence of singular values close to zero. If there is one zero singular value, it will not appear in matrix A+ (in this case the solution is not unique) and therefore there will not be any convergency problem. Several strategies have been proposed in the literature to attenuate matrix K+. We have chosen to give up (assume zero) the smallest singular values until the Q* convergency is achieved. Matrix x is the same one for all models having the same resistivity contrast between the layers. There will be a strictly zero singular value and its parameter eigenvector will be proportional to vector (l,, l,, . . . lo, 0,, . . . ON_ 1 ). Such a singular value will not appear in matrix li’. In Table I we show an example of inversion with the new error function. The experimental curve is for the model pi = 10, pz = 100, p3 = 2, p4 = 10,

531 TABLE

I

Successive

models

Iteration

p,

0

10.25 10.47 9.39 8.92 a.54 8.38 8.38 10.00

1 2 3 4 5 6 Shift

in an iterative

procedure

p3

4.20 2.81 2.16 1.85 1.72 1.68 1.68 2.01

25.02 37.80 53.29 67.68 79.40 84.54 84.56 100.89

9.28 8.98 9.27 8.93 8.55 8.38 8.38 10.00

of inversion

with the new error function’

t,

t*

t,

Error max. (%)

Q*

Q

1.00 2.30 1.79 1.97 1.99 2.00 2.00 2.00

10.00 7.20 5.10 3.95 3.22 2.97 2.97 2.97

68.13 31.08 29.17 26.04 25.26 25.14 25.11 25.11

0.1

0.008857 0.004419 0.000091 0.000031 0.000018 0.000024 0.000021 -

10V6

*I Iteration 0 corresponds to the initial model. ,Error max. is the maximum tween the experimental apparent resistivity and the theoretical one.

relative

error

3, t3 = 25. As initial model in the iterative procedure we have the one that corresponds to the maxima and minima of the experimental curve (iteration 0 in Table I). After six iterations, the process has been stopped because Q* does not decrease significantly and a very good adjustment is reached. The apparent resistivity curve corresponding to the final model has the same shape as the experimental one but vertically displaced in the bilogarithmic representation. Once the displacement of eq. 4 has been carried out, the shift model of Table I is obtained. The Q error function averaging the points in Table I is calculated according to: ti

= 2, t2 =

chosen

&=1/M

g i=

[l n Pai - in Fa(si$)l* 1

It must be noticed that although the solution obtained is one that corresponds to the experimental curve, this is purely a coincidence. As we will see later on, many equivalent models do exist for this case. In Tables II and III we present another example of inversion performed TABLE Example points

II of non-convergency

in inversion

with the classical

error function

averaging

Iteration

p1

pz

pa

pq

t,

t2

t,

Error max. (%)

Q

0 1 2 3 4

9.90 8.45 8.23 8.25 8.24

4.40 3.36 3.55 3.53 3.54

8.20 5.95 7.28 7.36 7.39

1.01 1.14 0.95 0.93 0.94

1.00 2.88 2.94 2.94 2.94

10.00 1.36 1.32 1.32 1.32

68.00 58.99 47.45 49.17 48.68

51 47 32 33 33

0.114294 0.052575 0.032354 0.032011 0.032001

the

be-

532 TABLE

III

Same example as of Table II but computed Iteration

0 1 2 3 4 5 6 7 8 Shift

P,

9.90

8.90 15.17 10.82 11.55 10.67 11.01 10.35 9.94 10.01

with the new error function

pi

p3

P*

t,

t,

t,

Error max. (%)

Q*

Q

4.40 3.85 2.30 2.40 1.63 1.54 1.06 0.99 0.93 0.94

8.20 6.71 9.28 11.99 16.04 20.27 28.03 34.03 39.07 39.34

1.01 1.57 1.11 1.16 1.20 1.08 1.10 1.04 1.00 1.00

1.00 2.47 1.24 2.06 2.74 2.39 2.53 2.50 2.50 2.50

10.00 1.43 1.69 1.77 2.54 2.56 1.81 1.86 1.84 1.84

68.00 65.41 33.41 37.78 26.60 18.75 14.19 10.94 9.15 9.15

0.2

0.042470 0.027353 0.019119 0.004737 0.002686 0.000081 0.000034 0.000006 0.000002 -

3*10--6

with the classical eq. 5 and the new error functions. The experimental curve corresponds to the model p1 = 10, pz = 1, p3 = 60, p4 = 1, t, = 2.5, t2 = 2, t3 = 6. As can be seen from Table II, the classical procedure is not convergent, whereas in Table III after 8 iterations a very fair solution is reached. After applying the displacement of eq. 4 the shift model is obtained, the apparent resistivities of which have differences smaller than 0.2% with respect to the experimental ones. Therefore, this model is equivalent to that corresponding to the experimental curve. GENERATING

EQUIVALENT

MODELS

As is known, the observation of the parameter eigenvectors can give information about the equivalence of a model (Johansen, 1977). With this new error function, we assume that two models p and p* are equivalent when the components mi@) and mi(p*) differ from each other by less than a preassigned quantity, for instance between the 3% and the 5% limits usually established for the errors on data. Let us consider a model p*. We can find another p-equivalent model considering an expansion in a Taylor series:

m’(p)=

m’@*)+&?_p*)

A base of the null space of matrix x will be a generator of increments 0, p*) such that the difference m$$) - RI@*) will be zero. The null space is generated by the parameter eigenvectors associated with the zero singular values. But since the classical forward problem for a stratified earth has a unique solution, there will be only one zero singular value which results from the fact that models with the same layer-resistivity contrast have the same value of slopes. Therefore the equivalence in classical

533

terms will correspond to the eigenvectors associated with the small singular values. We call pseudo null space the subspace generated by these eigenvectors. Any element of it will be a @ - p”) such that differences ml(p) rni(p’) be small. Thus this element will generate a p model equivalent top’. For instance, if u is a parameter eigenvector with a quasi-zero singular value and we want to find out a model p which is equivalent to the model p”, we may choose Ap = p - p” parallel to V. That is Ap = C-Gand from eq. 6 m’(p) = m!CpO) f 7iffF =m’Cp”)+ffXii = m’@O) for any real IYsuch that linearity of m’(p) between p and p” is valid. H is the data eigenvector associated with the x singular value (Lanczos, 1961). In general any combination of the eigenvectors belonging to the pseudo null space will generate equivalent models, as long as the linearity assumption is valid. From the A matrix decomposition a similar treatment is applied to generate equivalent models with the same apparent resistivities. f: f2 R p4 t1 tz

t3

h2=0.7020

h&=0.2973

x5=

0.0576

hs=o.o097

x7-0.0

Fig. 1. Singular values and parameter eigenvectors from the A matrix corresponding to iteration 6 of Tabie I.

534

Fig. 1 shows the singular values and parameter eigenvectors of the matrix A decomposition corresponding to the final model of iteration (eq. 6) on the example of Table I before the displacement is carried out. In Fig. 2 we present the singular values and parameter eigenvectors corresponding to the same model but calculated from matrix A. Comparison of the two figures is worthwhile. In matrix A decomposition (Fig. 2) we observed that there are two important equivalences defined by the two parameter eigenvectors u6 and u7 associated with the two smallest singular values. We observe that the main components of vector u7 are those corresponding to pz and tz which have the same value but with opposite sign, while other components are practically zero. In this case the so-called classic equivalence in T (Keller and Frischkneeht, 1966) is obtained, because an equivalent model will be generated by: p2 = P: exp (4) ta = t: exp (MY”,) and by multiplying

(7) both equations

we obtain the classical form:

Concerning vg, the main components correspond to p3 and t3 with a small contribution from p2 also. By modifying the modulus of the increment in the direction of ug we obtain the different equivalent models generated by this vector. In this case the equivalences obtained are not the classical ones, due to the small contribution of pz and due to the components p3 and t3 which have not the samevalue. From matrix A decomposition three equivalences are found. One corresponds to the zero singular value and the other two correspond to the very small eigenvalues. These singular values are associated respectively with eigenvectors U7, ii, and GS. In Fig. 1 it can be seen that v7 eigenvector is proportional to the (1,1,&l, O,O,O) vector. This means that any model with the same resistivity contrast between the layers has the same m’Cp) function. On the other hand E6 and G5 generate the same equivalence as the u7 and u6 eigenvectors of the A matrix. V6 can be decomposed into two vectors of different behaviour. One of them corresponds to the equivalence in T as in the case of eigenvector u7 and the other to the vertical displacement as u., : v, = uk f TU:

where the real 7 can be chosen such that z& be parallel to v7. In the same way the G5 eigenvector can be decomposed into two components, showing the same equivalences as u6 and &, This decomposition is always possible and the same equivalences as with matrix A decomposition are obtained.

535

hz= 2.6520

h3=2.0470

h=o

0864

h5=07061

‘1

V6

*-

1

I

-1

h=0.0577

11

0

V?

I

I

A7 =

0.0098

-1

Fig, 2. Singular values and parameter eigenvectors from the A matrix corresponding to iteration 6 of Table I.

In this example we have C6 = (0.2,-0.5,0.2,0.2,0,0.7,0) O,O,-O.?,O). Therefore T = 0.2 and the same equivalences vertical displacement, are found.

and u7 = (0,0.7,0, as in eq. 7 plus the

DISCUSSION

Experience has proved that convergency towards an acceptable solution is achieved with this new error function in cases where a classical error function evaluating the points of the apparent resistivity fails to convergence. Besides, convergency is generally faster, although it is not assured for all the initial models. Given an experimental curve, if convergency is not reached it may be due to two facts: the initial model is not adequate, or the experimental curve does not result from a layered earth model. The latter case has no solution unless the departure from the stratified model be small and can be introduced as a random error in the data. The former case can be solved simply by a better choice of the initial model. If not, an alternative inversion method can be tried. In fact, the new error function may be considered simply as an alternative when the classical error function gives rise to a non-acceptable local Q minimum.

636

Experience demonstrates also that when the final model is achieved, that is when Q* does not decrease any further, the solution can be refined through iterations with the classical error function. This is specially important when Q* is within a local minimum but has a too high value. In this case changing the error function may permit us to give up the local minimum and to fall into another minimum with a smaller Q* value. The combined use of the two error functions always gives a better adjustment. REFERENCES Bichara, M. and Lakshmanan, J., 1976. Fast automatic processing of resistivity soundings. Geophys. Prospect., 24: 354-370. inversion. Geophysics, 38: Inman, J.R., Ryu, J. and Ward, S.H., 1973. Resistivity 1088-1108. Johansen, H.K., 1977. A man/computer interpretation system for a resistivity sounding over a horizontally stratified earth. Geophys. Prospect., 25: 667-691. Jupp, D.L.3. and Vozoff, K., 1975. Stable iterative methods for the inversion of geophysical data. Geophys. J. R. Astron. Sot., 42: 952-976. Keller, G.V. and Frischknecht, F.C., 1966. Electrical Methods in Geophysical Prospecting. Pergamon Press, Oxford, 519 pp. Sounding Measurements, Geosounding Principles 1. Koefoed, O., 1979. Resistivity Elsevier, Amsterdam, 276 pp. Lanczos, C., 1961. Linear Differential Operators. Van Nostrand, London, 564 pp. Orellana, E. and Mooney, H.M., 1966. Master Tables and Curves for Vertical Electrical Sounding over Layered Structures. Interciencia, Madrid, 150 pp. Pow, J., Badiella, P. and Nifierola, J.M., 1984. Nuevo metodo de1 gradiente para la interpretaci6n de sondeos elkctricos. Rev. Geofis., 40: 147-158.