- Email: [email protected]
A generalized method to calculate standard curves for geoelectrical soundings
Geoexploration, Elsevier Science
21 (1983) 91-103 Publishers B.V., Amsterdam
A GENERALIZED GEOELECTRICAL
R. SANTINI
91 - Printed
in The Netherlands
METHOD TO CALCULATE SOUNDINGS
STANDARD
CURVES FOR
and R. ZAMBRANO
Flitel S.p.A., Via Julia, Capriccio di Vigonm, Padova (Italy) Istituto di Geofisica, Via Rudena, 3, Padova (Italy) (Received
December
30, 1981; accepted
after revision
July 30, 1982)
ABSTRACT Santini, R. and Zambrano, geoelectrical soundings.
R., 1983. A generalized method Geoexploration, 21: 91-103.
to calculate
standard
curves for
A generalized numerical method to calculate apparent resistivity curves has been developed. It consists in approximating the resistivity transform of the geoelectrical model by a linear combination of a low number of real exponentials. The apparent resistivity is then computed by means of a simple linear combination of elementary functions of which the mathematical expression depends on the particular electrode configuration considered. The method has been applied to Schlumberger, Wenner and Two-electrode arrays, and several examples are shown, which confirm the reliability of the proposed method for any geoelectrical model. Comparison with master curves denotes also a very good accuracy of the results, with low computing time.
INTRODUCTION
For many years several methods and computer programs to calculate apparent resistivity curves have been set up. Some of these methods have been conceived to be used in the field (Flathe, 1955; Van Dam, 1965); others require the use of a digital computer (Mooney and Wetzel, 1956; Onodera, 1963; Orellana and Mooney, 1966; Mooney et al., 1966; Nabighian, 1966; Argelo, 1967; Van Dam, 1967; Ghosh, 1971; Baranov, 1976; Brizzolari, 1977; Das Gupta et al., 1979; Das and Verma, 1980). The method shown here deals with the calculation of the resistivity curves in the case of horizontally stratified layers; it is first developed for the Schlumberger array and later is extended to the Wenner and Two-electrode arrays, considering also the case of insulating substratum. PROPOSED
METHOD
As is known,
0016-7141/83/$03.00
- SCHLUMBERGER
ARRAY
in the case of a horizontally
0 1983 Elsevier
Science
stratified
Publishers
medium
B.V.
consisting
of n
92
laterally homogeneous and isotropic layers having resistivity the Schlumberger apparent resistivity is given by:
pi and thickness di,
m
P&)
= s*
s
T(pi,di;X)Jl
(sX)XdX
0
where s is half the distance between the current electrodes, h the integration variable, J1 the Bessel function of order one, and T(pi,di;h) the resistivity transform (Koefoed, 1970). This function depends only on the parameters ~i,di and can be easily calculated by means of recurrence formulas, for example:
1 + Ti+l(h)pT’
(2)
tanh(dih)
starting from the last layer n, where T,,(h) = ,o~. Other recurrence formulas referring to T(X) or to similar functions closely linked to the resistivity transform, the so-called “kernel” functions introduced by Stefanescu, Slichter, Koefoed, etc. can be used (e.g., Sunde, 1949; Flathe, 1955; Van’Jan, 1962; Koefoed, 1965, 1970; Strakhov, 1966; Johansen, 1975). The proposed method consists of the following three steps. (a) Once the parameters pi and di have been assigned, the numerical values Tj of the transform at pre-established abscissae hj (j = 1, 2, . . . k) are computed by means of one of the above mentioned recurrence formulae. (b) The discrete function Tj = T(A.) is approximated by the linear combination of real exponentials: T,(h)
= p1 + g
ai exp(-qh)
(3)
i= 1
using the least square method. The si values (E > 0) are suitably preset so that only the coefficients ai must be calculated. The number m of exponent& depends on the required precision. (c) By substituting eq. 3 into 1 and considering that:
s
e-‘*hJ,(sX)dh
= S(E’ + s*)-
3/*
0 we
get an approximated apparent resistivity :
Pm,S(S)=Pl
+
E
Ui[l +
(Ej/S)2]-3’2
i= 1
where ai and ei have the same values assumed in eq. 3.
(4)
93
In Fig. 1 the approximating
functions:
(5)
g(f, A) = exp(4 and its transform
in the s-domain,
(6)
Fig. 1. The fitting function g( E,h) of the resistivity tion f( 6,s).
transform and the corresponding
func-
are plotted (versus A-’ ). In using log abscissae the shape of the above functions does not change when varying E: the curves are only shifted along the abscissae axis. With the present method both the calculation of the Bessel function and the numerical integration of eq. 1 are avoided, because only the numerical solution of a linear set of equations (resulting from the least square approximation of T(X)) is necessary. APPROXIMATION
OF THE RESISTIVITY
TRANSFORM
Let us suppose that both the k abscissae A. where T(A) will be sampled and the number m of fitting functions 5 have already been established. Then we want to determine the function T, (A) which approximates T(X) according to the least square criterion.
94
Therefore J$
we must have:
[Tm(hj)-
T(Xj)]* = min.
(7)
j=l
or: k
m
pi + C =[ j=l
i=
ai exp(-cjh)
-
T(Aj)
*
=
I
1
min.
(8)
The condition 8 would lead to a non-linear system in the 2m unknowns ai, I, which poses great difficulties for its numerical solution, and is very timeconsuming. Therefore we have preferred to pre-establish the ei values, so that only the coefficients ai are the unknowns. In so doing the following symmetric normal linear system of m equations for m unknowns is obtained: UI@I,I
+$01,2
+
. . . . +amol,rn
=771
a101,2
+a202,2
+
. . . . +%o*,rn
=772
(9)
+amum,m = 7)m
+a2U2,m + ....
alul,m
where: Uq,r
=
h
exp[-(E,
+
4)xjl
(10)
exp(-Eq$)
(11)
j=l
vq
=
5
[Wji)
-AI
j=l
Once the ai values are calculated by solving system 9, both the approximated resistivity transform T, (A) and the apparent resistivity function Pm, s (s) according to eq. 1 is entirely determined. Sampling of the resistivity transform Before the numerical computation of T(h) at a discrete number of points can be carried out, the interval and the law of the sampling must be taken into account. Theoretically we should consider an unlimited sampling interval as the resistivity transform extends on an infinite number of log cycles; yet the most significant portion of the transform function appears on a limited interval. The authors have taken as lower (left) extreme of the sampling interval: Xl
= 0.5 d,
95
putting x = A-‘. The value xk_1 of the upper (right) extreme of the above interval has been related to the pseudo-thickness d, and to the pseudo-resistivity peq of the layer equivalent to the first n - 1 layers of the considered geoelectrical model. Thus:
20 deqpnlpeq= 20 CP,
for Pn > Peq
20 &qPnlPeq
for
Xk_l=
= 20
RIP,
Pn < Peq
where d eq = m
;
Peq = JW
and R = ‘x1 pidi ; C = ‘x1 i=l i= 1
di/pi
(12)
indicate
the transverse .Resistance and the longitudinal Conductance of the first n - 1 layers respectively. However, it is useful for a better approximation of the resistivity transform at the right of the sampling interval to consider also the abscissa Xk = m, where T(A) is equal to pn. Towards this end, in doing the numerical computation it is enough to assign to Xk a very large value, for example lOlO, whereas it is not necessary to consider the abscissa at the origin, where the transform function is equal to pl, because this would be equivalent to adding an identically null equation in system 9. As far as the distribution of the sampling points is concerned it is preferable to take them approximately equally spaced along the log axis, similar to what is done as a rule in the field measurements of the apparent resistivity. This is equivalent to setting the values xi in a geometrical progression. With the proposed method, however, a rigorous equal spacing is not necessary, provided that the sampled points are reasonably distributed along the interval (XI, xk-1). Then it is enough to take 6 points for each log cycle, independently of the geoelectrical layering. Nevertheless a greater number of points can be considered without increasing the computer time. Number of fitting functions and their placing The determination of the number m of fitting functions 5, and consequently of the corresponding functions 6, is conditioned by two opposing requirements: on the one hand we want to obtain a good approximation to the resistivity transform (which generally improves as m increases) and on the other hand limit the number of eq. 9 on which depends the computing time. The experience acquired on different stratigraphic models suggests utilizing
96
normally 4-5 functions for each log cycle of the significant portion of T(h). It may be necessary to slightly increase this number in the case of very strong contrasts (>104) in the apparent resistivity curve or, which is the same, in the Dar-Zarrouk curve. Regarding the values E1, E*, , . E, , it is convenient for a better approximation that they also are taken in a geometrical progression. By so doing the values of functions 5 and 6 become “equally spaced” if plotted versus log abscissae (Fig. 1). On the limits ei, E, depends the width of the interval in which the above functions are placed: the greater the number of log cycles required to describe completely the resistivity transform, the larger must be the interval (ei, em ). Whatever be the parameters pi and di, the authors have found sufficient to take : =
fl
em=
2x, for h
;“k-1
> Peq
for Pn < Peq
i xk-lPn/Peq
The fitting functions
placed too far outside the above interval are wasted.
CASE OF SUBSTRATUM WITH INFINITE RESISTIVITY
When pn = m the point function T(&) cannot be approximated according to eq. 3 by a limited number of functions 5. Therefore other fitting functions have been used, which become infinite to order 1, like T(h) for large values of the abscissa x = h-l. It is known in fact that when pn = 00 and A-’ is sufficiently large we have: T(X) = x-‘/c C being the longitudinal conductance ing substratum. Thus the fitting function:
of the n - 1 layers lying on the insulat-
g(e, X) = exp(-eh)/eX has been considered, f(q,s) =
2.[l
(13) to which corresponds,
in the s-domain,
the function:
- (1 + s+2)--1’2]
(14)
E
calculated
by means of the transform
J e-%,(sh)dh
formula
= [(e’ + s~)“~ - e]/s(e’
2 and observing that:
+ s’)~‘~
0
The function
14 is an infinity
of the first order for large values of the abscissa s,
97
as is necessary
since
lim p&s(s) -s/C s+Also the functions 13 and 14, and the similar ones 5 and 6, maintain their shape when E varies if they are plotted versus log abscissae. To sum up, in the case pn = 00 the resistivity transform is approximated the following linear combination of real exponentials: Ui eXp(-fzih)/qh
Tm(N=P1 + g i=
(15)
1
The corresponding becomes: p,,Js)
by
= p, + g i= 1
approximated
ui;[l
expression
for the apparent resistiuity
(16)
- (1 +s2/+-1’21
i
The coefficients ci are calculated however, we have:
again by means of system 9 where now,
k
C expH~q + ~r)~l/~q+$
‘4,r =
(17)
j=l
rlq =
5
[T@j) -PI
1 eXP(-QXj)/eqAj
j=l
About the sampling interval and the placing of the fitting functions 13, considerations similar to those made for pn finite led us to take in every case the following standard values: Xl
= 0.5 dl
Cl
Xk-1
=
Em = 0.1 Xk-1
1
ACCURACY
20d,,
=2x,
= dl = 2 d,,
OF THE METHOD AND EXAMPLES
The procedure here proposed for calculating the apparent resistivity essentially consists of two steps: (1) to determine the function Tm (A) that approximates the sampled points T(A.); (2) to get, the apparent resistivity Pm,S(s) by means of a purely analytic operation expressed by the Hankel transform 2. Only in the first step errors (of approximation) can be introduced, which are transmitted to the next step, which is not affected by errors. However the authors have verified that in the second step the relative error in Tm (A) is usually transferred to pm,s(s) with a sensible amplification that sometimes can reach the value of several tens.
98
That is just the opposite of what is observed in the inverse problem of computing the resistivity transform from apparent resistivity curves (Santini and Zambrano, 1981) where an attenuation of the error occurs. Nevertheless if 5 fitting functions per log cycle are taken, a relative error in the apparent resistivity lower than 10m3 results even in the most critical case which can be encountered in practice. In order to check the accuracy of the present method several resistivity curves for 3 and 4 layers have been calculated and then compared to those given in the Orellana-Mooney collection (1966). In every test we observed that the accuracy of the ,LJ~,s(s) values was, point by point, not lower than the accuracy of the master values. For example, in the layering: p = 1-20-0.1-1
;
d = l-5-25
and p = l-0,025-
;
d = l-25,
we had a mean relative error of 5.2 10W4 in the first case, and 5.1 l.0F3 in the second case. It is interesting to note that the mean relative error in the approximation of T(h) was 2.1*10-* and 3.7. 10W4, which signifies an error amplification of 25 and 14, respectively. In the above examples 20 and 11 fitting functions have been utilized consuming 0.8 and 0.6 s on an IBM 370/158 computer. l
l
EXTENSION OF THE METHOD TO OTHER ELE~RODE
The proposed method can be extended in theory Towards this end we shall follow again the criterion resistivity transform function by the combination:
~ONFI~URA~ONS
to any configuration. of approximating the
i= 1
to which every time a specific expression resistivity:
for the (approximated)
apparent
m
will correspond by virtue of the Hankel transform. As far as the approximation of T(X) is concerned it is possible to use the functions 5 when ~a is finite and the functions 13 when pn is infinite, since the resistivity transform only depends on the parameters of the geoelectrical model. For the same reason the criteria about the sampling of T(X) and the posi-
99
tions of the fitting functions do not change, and the authors have verified this for the Wenner, two-electrode and dipolar arrays. Nevertheless, they noted that among these configurations the dipole arrays (except obviously the azimuthal one) present the largest amplification of the error and hence the smallest final accuracy of apparent resistivity. The writers intend to study an improvement of the procedure for the last arrays in another article. WENNERARRAY
Case of substratum
with p,, finite
For the Wenner configuration P&W(S) =
2s
I
T(X)[J,(Xs)
the apparent
resistivity
is:
- J,,(2Xs)] dX
(19)
0
where s denotes the distance Bessel function of first kind Substituting in eq. 19 the of T(h) and keeping in mind
between two adjacent electrodes and JO the and order zero. approximation !I’, (A), again given by 3, in place the known relation:
OD
s
e-‘“J&h)dh
= (E’ + s*)-~/*
(20)
0
Wenner apparent
we have for the (approximated) P m, w(s) =
p, + 5
aj[2(1
+ ef/s*)-‘I*
resistiuity :
- (1 + ef/4s*)-‘/*]
(21)
i=l
Using the latter formula p = 100-5000-2-5000-2
the curve of Fig. 2, corresponding
to the layering:
d = 1-5.3-15.8-60
has been calculated. The Dar-Zarrouk curve has been drawn as a useful reference and comparison, but it is not required by the method. In spite of the great “extension” of the curve - 8 log cycles - and the very strong resistivity contrasts, the computing time was 1.7 s, and 24 fitting functions were used. This example shows also that the present method does not need the thicknesses to be integers, and multiples of a common thickness. Case of substratum
with pn = m
In this case the resistivity of eq. 13.
transform
will be approximated
by the function
100 *
0 -i-
i-l
t
---_
r-
Rtsistivity
t
Dar
transform
Zarrouk
Apparent True
curve
resistivity
resistivity
meters
Fig. 2. Wenner apparent resistivity curves: above, p = 100-5000-2-5000-2, 1-5.3-15.8-60; below, p = 0.025-lO-0.025-9 d = 1-194-5000.
d =
In virtue of eq. 20 we can write: m
s
eMf”[&@A) - J,,(Bsh)] dh = (e’ + s~)-~/~ - (e’ + 4s2)-‘12
0
Integrating both members of the above eq. with respect to the parameter E we get: 91 s E
m s 0
eRfh [JO(sX) - J0(2sh)] dXde =
s
[(r? + sy--1/2 - (8 + 4s 2) -‘/2]dE
c
Then, inverting the order of integration in the first member and solving the
101
integral with respect to E, we have:
s
oI F
[Jo(sh) - J,(2sh)]dX
= log[(s + (e* + 4s’)“*)f(e
+ (e2 + s’)‘/~)]
0
Therefore, considering eqs. 15 and 19 we easily attain at least:
P,,w(s)=P1 +2
5 (li+ log[(l i= 1
+ (1 + 4s2/e:)“*)/(l
f (1 + 5*/e: )“2)l
(22)
i
An application example appears in Fig. 2 for the layering: p = 0.025-10-0.025-
d = 1-199-5000.
Twenty-four functions were used, consuming 1.6 s. TWO-ELECTRODE
ARRAY
In such a confi~rati~n, where one current electrode and one potenti~ electrode are put at “infinity”, the apparent resistivity is given by the integral relation:
PAT
=s j
V)JofsWX
(23)
0
being the distance between the active current electrode and the active potential one. A procedure similar to that carried out in the Wenner case leads to the following approximated expression for the apparent resistivity:
s
p,,$s)
= p, +
5 up
f qs*,+‘*
(24)
i=l
On the contrary it is impossible to treat in the same way the case pn = 00, because we arrive at the integral:
r= eeh J,(sh)dh ji
0
x
which does not converge. This unfavourable result does not depend on the particular choice of function 13, but rather on the intrinsic ch~ac~~stics of this array. In fact the resistivity transform becomes infinite of order 1 at the origin, as already seen, while J,,(O) = 1, and thus the eq. 23 diverges.
102
In Fig. 3 two examples
corresponding
p = l-0.2~-5-0.01
a! = l-3-10,
p = l-0.2--1000
d = l-25
to the layering:
are shown. In the first example 15 fitting functions were used, and 33 in the second one, consuming 0.7 and 2.2 s, respectively. The latter example shows that, if we assign to pn a high value, the impossibility of dealing with the theoretical case pn = 03 does not imply a practical limitation.
1-l
0.1
1
10
102
103
104
105
106
s
Cr
107
meters
Fig. 3. Two-electrode apparent resistivity curves: above, p = 1-0.2-5-0.01, below, p = 1-0.2-1000, d = l-25. (For the legend see Fig. 2.)
d = l- .3-10;
REFERENCES Argelo, S.M., 1967. Two computer programs for the calculation of standard graphs for resistivity prospecting. Geophys. Prospect., 15: 71-91. Baranov, W., 1976. Calcul des courbes de sondages Blectriques a l’aide des fonctions d’dchantillonage. Geophys. Prospect., 24: 222-232. Brizzolari, E., 1977. Calcolo di curve teoriche di resistivita apparente. Boll. Geofis. Teor. Appl., 20(73-74): 27-32.
103 Das Gupta, S.P., Barun Prosad, P. and Buddhadeb Banerjee, 1979. Rapid computation of resistivity sounding master curves. Geoexploration, 17: 223-228. Das, U.C. and Verma, S.K., 1980. Digital linear filter for computing type curves for the two-electrode system of resistivity sounding. Geophys. Prospect., 28: 610-619. Flathe, H., 1955. A practical method of calculating geoelectrical model graphs for horizontally stratified media. Geophys. Prospect., 3: 268-294. Ghosh, D.P., 1971. Inverse filter coefficients for the computation of apparent resistivity standard curves for a horizontally stratified earth. Geophys. Prospect., 19: 769-775. Johansen, H.K., 1975. An interactive computer graphic-display-terminal system for interpretation of resistivity soundings. Geophys. Prospect., 23: 449-458. Koefoed, O., 1965. Direct methods of interpreting resistivity observations. Geophys. Prospect., 13: 568-592. Koefoed, O., 1970. A fast method for determining the layer distribution from the raised kernel function in geoelectrical sounding. Geophys. Prospect., 18: 564-570. Mooney, H.M. and Wetzel, W.W., 1956. The Potentials about a Point Electrode and Apparent Resistivity Cu::ves for a Two, Three, and Four Layered Earth. Univ. Minn. Press, Minneapolis, 146 pp. Mooney, H.M., Orellana, E., Pickett, H. and Tornheim, L., 1966. A resistivity computation method for layered earth model. Geophysics, 31: 192-203. Nabighian, M.N., 1966. The application of finite forward differences in the resistivity computations over a layered earth. Geophysics, 31: 971-980. Onodera, S., 1963. Numerical analysis of relative resistivity for a horizontally layered earth. Geophysics, 28: 222-231. Orellana, E. and Mooney, H.M., 1966. Master Tables and Curves for Vertical Electrical Sounding over Layered Structures. Interciencia, Madrid, 34 pp. Santini, R. and Zambrano, R., 1981. A numerical method of calculating the kernel function from Schlumberger apparent resistivity data. Geophys. Prospect., 29: 108-127. Strakhov, V.N., 1966. The analytic determination of the parameters of a horizontally layered medium from data obtained by vertical electrical prospecting. Izv. Akad. Nauk, USSR, Earth Phys., 4: 237-242. Sunde, E.D., 1949. Earth Conductive Effects in Transmission Systems. Van Nostrand, New York, NY. Van Dam, J.C., 1965. A simple method for the calculation of standard-graphs to be used in geoelectrical prospecting. Geophys. Prospect., 13: 37-65. Van Dam, J.C., 1967. Mathematical denotation of standard-graphs for resistivity prospecting in view of their calculation by means of a digital computer. Geophys. Prospect., 15: 57-70. Van’Jan, L.I., 1962. On the calculation of theoretical electrical sounding curve. Prikl. Geofiz., 34: 135-144. (In Russian).
21 (1983) 91-103 Publishers B.V., Amsterdam
A GENERALIZED GEOELECTRICAL
R. SANTINI
91 - Printed
in The Netherlands
METHOD TO CALCULATE SOUNDINGS
STANDARD
CURVES FOR
and R. ZAMBRANO
Flitel S.p.A., Via Julia, Capriccio di Vigonm, Padova (Italy) Istituto di Geofisica, Via Rudena, 3, Padova (Italy) (Received
December
30, 1981; accepted
after revision
July 30, 1982)
ABSTRACT Santini, R. and Zambrano, geoelectrical soundings.
R., 1983. A generalized method Geoexploration, 21: 91-103.
to calculate
standard
curves for
A generalized numerical method to calculate apparent resistivity curves has been developed. It consists in approximating the resistivity transform of the geoelectrical model by a linear combination of a low number of real exponentials. The apparent resistivity is then computed by means of a simple linear combination of elementary functions of which the mathematical expression depends on the particular electrode configuration considered. The method has been applied to Schlumberger, Wenner and Two-electrode arrays, and several examples are shown, which confirm the reliability of the proposed method for any geoelectrical model. Comparison with master curves denotes also a very good accuracy of the results, with low computing time.
INTRODUCTION
For many years several methods and computer programs to calculate apparent resistivity curves have been set up. Some of these methods have been conceived to be used in the field (Flathe, 1955; Van Dam, 1965); others require the use of a digital computer (Mooney and Wetzel, 1956; Onodera, 1963; Orellana and Mooney, 1966; Mooney et al., 1966; Nabighian, 1966; Argelo, 1967; Van Dam, 1967; Ghosh, 1971; Baranov, 1976; Brizzolari, 1977; Das Gupta et al., 1979; Das and Verma, 1980). The method shown here deals with the calculation of the resistivity curves in the case of horizontally stratified layers; it is first developed for the Schlumberger array and later is extended to the Wenner and Two-electrode arrays, considering also the case of insulating substratum. PROPOSED
METHOD
As is known,
0016-7141/83/$03.00
- SCHLUMBERGER
ARRAY
in the case of a horizontally
0 1983 Elsevier
Science
stratified
Publishers
medium
B.V.
consisting
of n
92
laterally homogeneous and isotropic layers having resistivity the Schlumberger apparent resistivity is given by:
pi and thickness di,
m
P&)
= s*
s
T(pi,di;X)Jl
(sX)XdX
0
where s is half the distance between the current electrodes, h the integration variable, J1 the Bessel function of order one, and T(pi,di;h) the resistivity transform (Koefoed, 1970). This function depends only on the parameters ~i,di and can be easily calculated by means of recurrence formulas, for example:
1 + Ti+l(h)pT’
(2)
tanh(dih)
starting from the last layer n, where T,,(h) = ,o~. Other recurrence formulas referring to T(X) or to similar functions closely linked to the resistivity transform, the so-called “kernel” functions introduced by Stefanescu, Slichter, Koefoed, etc. can be used (e.g., Sunde, 1949; Flathe, 1955; Van’Jan, 1962; Koefoed, 1965, 1970; Strakhov, 1966; Johansen, 1975). The proposed method consists of the following three steps. (a) Once the parameters pi and di have been assigned, the numerical values Tj of the transform at pre-established abscissae hj (j = 1, 2, . . . k) are computed by means of one of the above mentioned recurrence formulae. (b) The discrete function Tj = T(A.) is approximated by the linear combination of real exponentials: T,(h)
= p1 + g
ai exp(-qh)
(3)
i= 1
using the least square method. The si values (E > 0) are suitably preset so that only the coefficients ai must be calculated. The number m of exponent& depends on the required precision. (c) By substituting eq. 3 into 1 and considering that:
s
e-‘*hJ,(sX)dh
= S(E’ + s*)-
3/*
0 we
get an approximated apparent resistivity :
Pm,S(S)=Pl
+
E
Ui[l +
(Ej/S)2]-3’2
i= 1
where ai and ei have the same values assumed in eq. 3.
(4)
93
In Fig. 1 the approximating
functions:
(5)
g(f, A) = exp(4 and its transform
in the s-domain,
(6)
Fig. 1. The fitting function g( E,h) of the resistivity tion f( 6,s).
transform and the corresponding
func-
are plotted (versus A-’ ). In using log abscissae the shape of the above functions does not change when varying E: the curves are only shifted along the abscissae axis. With the present method both the calculation of the Bessel function and the numerical integration of eq. 1 are avoided, because only the numerical solution of a linear set of equations (resulting from the least square approximation of T(X)) is necessary. APPROXIMATION
OF THE RESISTIVITY
TRANSFORM
Let us suppose that both the k abscissae A. where T(A) will be sampled and the number m of fitting functions 5 have already been established. Then we want to determine the function T, (A) which approximates T(X) according to the least square criterion.
94
Therefore J$
we must have:
[Tm(hj)-
T(Xj)]* = min.
(7)
j=l
or: k
m
pi + C =[ j=l
i=
ai exp(-cjh)
-
T(Aj)
*
=
I
1
min.
(8)
The condition 8 would lead to a non-linear system in the 2m unknowns ai, I, which poses great difficulties for its numerical solution, and is very timeconsuming. Therefore we have preferred to pre-establish the ei values, so that only the coefficients ai are the unknowns. In so doing the following symmetric normal linear system of m equations for m unknowns is obtained: UI@I,I
+$01,2
+
. . . . +amol,rn
=771
a101,2
+a202,2
+
. . . . +%o*,rn
=772
(9)
+amum,m = 7)m
+a2U2,m + ....
alul,m
where: Uq,r
=
h
exp[-(E,
+
4)xjl
(10)
exp(-Eq$)
(11)
j=l
vq
=
5
[Wji)
-AI
j=l
Once the ai values are calculated by solving system 9, both the approximated resistivity transform T, (A) and the apparent resistivity function Pm, s (s) according to eq. 1 is entirely determined. Sampling of the resistivity transform Before the numerical computation of T(h) at a discrete number of points can be carried out, the interval and the law of the sampling must be taken into account. Theoretically we should consider an unlimited sampling interval as the resistivity transform extends on an infinite number of log cycles; yet the most significant portion of the transform function appears on a limited interval. The authors have taken as lower (left) extreme of the sampling interval: Xl
= 0.5 d,
95
putting x = A-‘. The value xk_1 of the upper (right) extreme of the above interval has been related to the pseudo-thickness d, and to the pseudo-resistivity peq of the layer equivalent to the first n - 1 layers of the considered geoelectrical model. Thus:
20 deqpnlpeq= 20 CP,
for Pn > Peq
20 &qPnlPeq
for
Xk_l=
= 20
RIP,
Pn < Peq
where d eq = m
;
Peq = JW
and R = ‘x1 pidi ; C = ‘x1 i=l i= 1
di/pi
(12)
indicate
the transverse .Resistance and the longitudinal Conductance of the first n - 1 layers respectively. However, it is useful for a better approximation of the resistivity transform at the right of the sampling interval to consider also the abscissa Xk = m, where T(A) is equal to pn. Towards this end, in doing the numerical computation it is enough to assign to Xk a very large value, for example lOlO, whereas it is not necessary to consider the abscissa at the origin, where the transform function is equal to pl, because this would be equivalent to adding an identically null equation in system 9. As far as the distribution of the sampling points is concerned it is preferable to take them approximately equally spaced along the log axis, similar to what is done as a rule in the field measurements of the apparent resistivity. This is equivalent to setting the values xi in a geometrical progression. With the proposed method, however, a rigorous equal spacing is not necessary, provided that the sampled points are reasonably distributed along the interval (XI, xk-1). Then it is enough to take 6 points for each log cycle, independently of the geoelectrical layering. Nevertheless a greater number of points can be considered without increasing the computer time. Number of fitting functions and their placing The determination of the number m of fitting functions 5, and consequently of the corresponding functions 6, is conditioned by two opposing requirements: on the one hand we want to obtain a good approximation to the resistivity transform (which generally improves as m increases) and on the other hand limit the number of eq. 9 on which depends the computing time. The experience acquired on different stratigraphic models suggests utilizing
96
normally 4-5 functions for each log cycle of the significant portion of T(h). It may be necessary to slightly increase this number in the case of very strong contrasts (>104) in the apparent resistivity curve or, which is the same, in the Dar-Zarrouk curve. Regarding the values E1, E*, , . E, , it is convenient for a better approximation that they also are taken in a geometrical progression. By so doing the values of functions 5 and 6 become “equally spaced” if plotted versus log abscissae (Fig. 1). On the limits ei, E, depends the width of the interval in which the above functions are placed: the greater the number of log cycles required to describe completely the resistivity transform, the larger must be the interval (ei, em ). Whatever be the parameters pi and di, the authors have found sufficient to take : =
fl
em=
2x, for h
;“k-1
> Peq
for Pn < Peq
i xk-lPn/Peq
The fitting functions
placed too far outside the above interval are wasted.
CASE OF SUBSTRATUM WITH INFINITE RESISTIVITY
When pn = m the point function T(&) cannot be approximated according to eq. 3 by a limited number of functions 5. Therefore other fitting functions have been used, which become infinite to order 1, like T(h) for large values of the abscissa x = h-l. It is known in fact that when pn = 00 and A-’ is sufficiently large we have: T(X) = x-‘/c C being the longitudinal conductance ing substratum. Thus the fitting function:
of the n - 1 layers lying on the insulat-
g(e, X) = exp(-eh)/eX has been considered, f(q,s) =
2.[l
(13) to which corresponds,
in the s-domain,
the function:
- (1 + s+2)--1’2]
(14)
E
calculated
by means of the transform
J e-%,(sh)dh
formula
= [(e’ + s~)“~ - e]/s(e’
2 and observing that:
+ s’)~‘~
0
The function
14 is an infinity
of the first order for large values of the abscissa s,
97
as is necessary
since
lim p&s(s) -s/C s+Also the functions 13 and 14, and the similar ones 5 and 6, maintain their shape when E varies if they are plotted versus log abscissae. To sum up, in the case pn = 00 the resistivity transform is approximated the following linear combination of real exponentials: Ui eXp(-fzih)/qh
Tm(N=P1 + g i=
(15)
1
The corresponding becomes: p,,Js)
by
= p, + g i= 1
approximated
ui;[l
expression
for the apparent resistiuity
(16)
- (1 +s2/+-1’21
i
The coefficients ci are calculated however, we have:
again by means of system 9 where now,
k
C expH~q + ~r)~l/~q+$
‘4,r =
(17)
j=l
rlq =
5
[T@j) -PI
1 eXP(-QXj)/eqAj
j=l
About the sampling interval and the placing of the fitting functions 13, considerations similar to those made for pn finite led us to take in every case the following standard values: Xl
= 0.5 dl
Cl
Xk-1
=
Em = 0.1 Xk-1
1
ACCURACY
20d,,
=2x,
= dl = 2 d,,
OF THE METHOD AND EXAMPLES
The procedure here proposed for calculating the apparent resistivity essentially consists of two steps: (1) to determine the function Tm (A) that approximates the sampled points T(A.); (2) to get, the apparent resistivity Pm,S(s) by means of a purely analytic operation expressed by the Hankel transform 2. Only in the first step errors (of approximation) can be introduced, which are transmitted to the next step, which is not affected by errors. However the authors have verified that in the second step the relative error in Tm (A) is usually transferred to pm,s(s) with a sensible amplification that sometimes can reach the value of several tens.
98
That is just the opposite of what is observed in the inverse problem of computing the resistivity transform from apparent resistivity curves (Santini and Zambrano, 1981) where an attenuation of the error occurs. Nevertheless if 5 fitting functions per log cycle are taken, a relative error in the apparent resistivity lower than 10m3 results even in the most critical case which can be encountered in practice. In order to check the accuracy of the present method several resistivity curves for 3 and 4 layers have been calculated and then compared to those given in the Orellana-Mooney collection (1966). In every test we observed that the accuracy of the ,LJ~,s(s) values was, point by point, not lower than the accuracy of the master values. For example, in the layering: p = 1-20-0.1-1
;
d = l-5-25
and p = l-0,025-
;
d = l-25,
we had a mean relative error of 5.2 10W4 in the first case, and 5.1 l.0F3 in the second case. It is interesting to note that the mean relative error in the approximation of T(h) was 2.1*10-* and 3.7. 10W4, which signifies an error amplification of 25 and 14, respectively. In the above examples 20 and 11 fitting functions have been utilized consuming 0.8 and 0.6 s on an IBM 370/158 computer. l
l
EXTENSION OF THE METHOD TO OTHER ELE~RODE
The proposed method can be extended in theory Towards this end we shall follow again the criterion resistivity transform function by the combination:
~ONFI~URA~ONS
to any configuration. of approximating the
i= 1
to which every time a specific expression resistivity:
for the (approximated)
apparent
m
will correspond by virtue of the Hankel transform. As far as the approximation of T(X) is concerned it is possible to use the functions 5 when ~a is finite and the functions 13 when pn is infinite, since the resistivity transform only depends on the parameters of the geoelectrical model. For the same reason the criteria about the sampling of T(X) and the posi-
99
tions of the fitting functions do not change, and the authors have verified this for the Wenner, two-electrode and dipolar arrays. Nevertheless, they noted that among these configurations the dipole arrays (except obviously the azimuthal one) present the largest amplification of the error and hence the smallest final accuracy of apparent resistivity. The writers intend to study an improvement of the procedure for the last arrays in another article. WENNERARRAY
Case of substratum
with p,, finite
For the Wenner configuration P&W(S) =
2s
I
T(X)[J,(Xs)
the apparent
resistivity
is:
- J,,(2Xs)] dX
(19)
0
where s denotes the distance Bessel function of first kind Substituting in eq. 19 the of T(h) and keeping in mind
between two adjacent electrodes and JO the and order zero. approximation !I’, (A), again given by 3, in place the known relation:
OD
s
e-‘“J&h)dh
= (E’ + s*)-~/*
(20)
0
Wenner apparent
we have for the (approximated) P m, w(s) =
p, + 5
aj[2(1
+ ef/s*)-‘I*
resistiuity :
- (1 + ef/4s*)-‘/*]
(21)
i=l
Using the latter formula p = 100-5000-2-5000-2
the curve of Fig. 2, corresponding
to the layering:
d = 1-5.3-15.8-60
has been calculated. The Dar-Zarrouk curve has been drawn as a useful reference and comparison, but it is not required by the method. In spite of the great “extension” of the curve - 8 log cycles - and the very strong resistivity contrasts, the computing time was 1.7 s, and 24 fitting functions were used. This example shows also that the present method does not need the thicknesses to be integers, and multiples of a common thickness. Case of substratum
with pn = m
In this case the resistivity of eq. 13.
transform
will be approximated
by the function
100 *
0 -i-
i-l
t
---_
r-
Rtsistivity
t
Dar
transform
Zarrouk
Apparent True
curve
resistivity
resistivity
meters
Fig. 2. Wenner apparent resistivity curves: above, p = 100-5000-2-5000-2, 1-5.3-15.8-60; below, p = 0.025-lO-0.025-9 d = 1-194-5000.
d =
In virtue of eq. 20 we can write: m
s
eMf”[&@A) - J,,(Bsh)] dh = (e’ + s~)-~/~ - (e’ + 4s2)-‘12
0
Integrating both members of the above eq. with respect to the parameter E we get: 91 s E
m s 0
eRfh [JO(sX) - J0(2sh)] dXde =
s
[(r? + sy--1/2 - (8 + 4s 2) -‘/2]dE
c
Then, inverting the order of integration in the first member and solving the
101
integral with respect to E, we have:
s
oI F
[Jo(sh) - J,(2sh)]dX
= log[(s + (e* + 4s’)“*)f(e
+ (e2 + s’)‘/~)]
0
Therefore, considering eqs. 15 and 19 we easily attain at least:
P,,w(s)=P1 +2
5 (li+ log[(l i= 1
+ (1 + 4s2/e:)“*)/(l
f (1 + 5*/e: )“2)l
(22)
i
An application example appears in Fig. 2 for the layering: p = 0.025-10-0.025-
d = 1-199-5000.
Twenty-four functions were used, consuming 1.6 s. TWO-ELECTRODE
ARRAY
In such a confi~rati~n, where one current electrode and one potenti~ electrode are put at “infinity”, the apparent resistivity is given by the integral relation:
PAT
=s j
V)JofsWX
(23)
0
being the distance between the active current electrode and the active potential one. A procedure similar to that carried out in the Wenner case leads to the following approximated expression for the apparent resistivity:
s
p,,$s)
= p, +
5 up
f qs*,+‘*
(24)
i=l
On the contrary it is impossible to treat in the same way the case pn = 00, because we arrive at the integral:
r= eeh J,(sh)dh ji
0
x
which does not converge. This unfavourable result does not depend on the particular choice of function 13, but rather on the intrinsic ch~ac~~stics of this array. In fact the resistivity transform becomes infinite of order 1 at the origin, as already seen, while J,,(O) = 1, and thus the eq. 23 diverges.
102
In Fig. 3 two examples
corresponding
p = l-0.2~-5-0.01
a! = l-3-10,
p = l-0.2--1000
d = l-25
to the layering:
are shown. In the first example 15 fitting functions were used, and 33 in the second one, consuming 0.7 and 2.2 s, respectively. The latter example shows that, if we assign to pn a high value, the impossibility of dealing with the theoretical case pn = 03 does not imply a practical limitation.
1-l
0.1
1
10
102
103
104
105
106
s
Cr
107
meters
Fig. 3. Two-electrode apparent resistivity curves: above, p = 1-0.2-5-0.01, below, p = 1-0.2-1000, d = l-25. (For the legend see Fig. 2.)
d = l- .3-10;
REFERENCES Argelo, S.M., 1967. Two computer programs for the calculation of standard graphs for resistivity prospecting. Geophys. Prospect., 15: 71-91. Baranov, W., 1976. Calcul des courbes de sondages Blectriques a l’aide des fonctions d’dchantillonage. Geophys. Prospect., 24: 222-232. Brizzolari, E., 1977. Calcolo di curve teoriche di resistivita apparente. Boll. Geofis. Teor. Appl., 20(73-74): 27-32.
103 Das Gupta, S.P., Barun Prosad, P. and Buddhadeb Banerjee, 1979. Rapid computation of resistivity sounding master curves. Geoexploration, 17: 223-228. Das, U.C. and Verma, S.K., 1980. Digital linear filter for computing type curves for the two-electrode system of resistivity sounding. Geophys. Prospect., 28: 610-619. Flathe, H., 1955. A practical method of calculating geoelectrical model graphs for horizontally stratified media. Geophys. Prospect., 3: 268-294. Ghosh, D.P., 1971. Inverse filter coefficients for the computation of apparent resistivity standard curves for a horizontally stratified earth. Geophys. Prospect., 19: 769-775. Johansen, H.K., 1975. An interactive computer graphic-display-terminal system for interpretation of resistivity soundings. Geophys. Prospect., 23: 449-458. Koefoed, O., 1965. Direct methods of interpreting resistivity observations. Geophys. Prospect., 13: 568-592. Koefoed, O., 1970. A fast method for determining the layer distribution from the raised kernel function in geoelectrical sounding. Geophys. Prospect., 18: 564-570. Mooney, H.M. and Wetzel, W.W., 1956. The Potentials about a Point Electrode and Apparent Resistivity Cu::ves for a Two, Three, and Four Layered Earth. Univ. Minn. Press, Minneapolis, 146 pp. Mooney, H.M., Orellana, E., Pickett, H. and Tornheim, L., 1966. A resistivity computation method for layered earth model. Geophysics, 31: 192-203. Nabighian, M.N., 1966. The application of finite forward differences in the resistivity computations over a layered earth. Geophysics, 31: 971-980. Onodera, S., 1963. Numerical analysis of relative resistivity for a horizontally layered earth. Geophysics, 28: 222-231. Orellana, E. and Mooney, H.M., 1966. Master Tables and Curves for Vertical Electrical Sounding over Layered Structures. Interciencia, Madrid, 34 pp. Santini, R. and Zambrano, R., 1981. A numerical method of calculating the kernel function from Schlumberger apparent resistivity data. Geophys. Prospect., 29: 108-127. Strakhov, V.N., 1966. The analytic determination of the parameters of a horizontally layered medium from data obtained by vertical electrical prospecting. Izv. Akad. Nauk, USSR, Earth Phys., 4: 237-242. Sunde, E.D., 1949. Earth Conductive Effects in Transmission Systems. Van Nostrand, New York, NY. Van Dam, J.C., 1965. A simple method for the calculation of standard-graphs to be used in geoelectrical prospecting. Geophys. Prospect., 13: 37-65. Van Dam, J.C., 1967. Mathematical denotation of standard-graphs for resistivity prospecting in view of their calculation by means of a digital computer. Geophys. Prospect., 15: 57-70. Van’Jan, L.I., 1962. On the calculation of theoretical electrical sounding curve. Prikl. Geofiz., 34: 135-144. (In Russian).