- Email: [email protected]
The versatility of digital linear filters used in computing resistivity and em sounding curves
Geoexploration, 18 (1981) 297-310 Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands
297
THE VERSATILITY OF DIGITAL LINEAR FILTERS USED IN COMPUTING RESISTIVITY AND EM SOUNDING CURVES
U.C. DAS and SK. VERMA National Geophysical Research Institute, Hyderabad (India) (Received May 7, 1980; accepted October 5, 1980)
ABSTRACT Das, UC. and Verma, SK., 1981. The versatility of digital linear filters used in computing resistivity and EM sounding curves. Geoexploration, 18: 297-310. The applicability of EM filters to compute resistivity sounding curves and that of resistivity filters in computing EM sounding curves is shown. This has been possible mainly because of the flexible nature of the filtering process involved which permits us to redefine an input function in order to suit a particular filter function. The desired output can then be obtained by convolving the modified input with the chosen filter. This flexibility can be used advantageously to compute the electrical and EM response of a layered earth to a number of source-receiver configurations by using the resistivity and the EM filters individually or jointly, as required. The study leads us to establish the relationships among various coil configurations used in EM sounding. A great computational advantage accruing from the present study is that the nine point filter to compute Schlumberger resistivity curves can also be used to compute EM sounding curves for perpendicular and vertical coplanar loop systems on programmable pocket calculators such as the HP-67 or TI-59.
INTRODUCTION
The extensive use of electric filters in communization engineering in order to clean signals contaminated with noise is well known. In geophysical operations the same technique is utilized for recording the seismic signals. However, the application of electric filter theory to geophysical data processing and inte~re~tion was possible only after Swartz (1954) indicated the similarity between the electric filter and the transformation of geophysical quantities. It was then realized that the operators used to convert residual gravity maps into their derivatives, or to continue the residual maps upward or downward were in fact linear filtering operations. Dean (1958) put this thinking on a sound basis. The most ou~t~~ng applications of linear filter theory have been made in the field of seismic data processing and the credit goes jointly to Robinson and Treital (1973). However, in the field of resistivity interpretation, it is quite surprising that 0016-71~2~81/0000-OOOOf$O2.50
01981
Etsevier Scientific Publishing Company
298
although active research has been going on since as early as the 1930’s, it was as late as the 1970’s (Ghosh, 1970; 1971a; 1971b) that the principle of linear filtering was applied to the problem of direct interpretation of Schlumberger, Wenner and dipole-dipole sounding data (Ghosh, 1971a; Das et al., 1974) and later for the computation of theoretical resistivity and EM sounding curves (Ghosh, 1971b; Das and Ghosh, 1974; Kumar and Das, personal communication, 1975; Johansen, 1975; Nyman and Landisman, 1977; Koefoed et al., 1972). Expressions for various quantities such as apparent resistivity, EM fields or mutual impedance ratio etc., in d.c. and EM prospecting involve Bessel functions of orders 0 and 1. A suitable change of variables in these Hankel transforms results in convolution integrals whereby filters are developed to evaluate these integrals. Computing these Bessel function filters (J, and J,), not associated with any multiplying factor as in resistivity and EM filters referred earlier, at identical abscissae, Anderson (1979) has shown how to compute the related Hankel transforms in which the discrete values of the initial kernels are stored for repeated use in the evaluation of other integrals. Avoidance of re-evaluation of kernels thus saves substantial computer time. The availability of these Bessel function filters (Jo and J,) will lead to the computation of any Hankel transform. In this paper, however, we show a new aspect of the available resistivity and EM filters - their versatility. By simple manipulation of the kernels we can use EM filters to compute resistivity curves or resistivity filters to compute EM sounding curves. The advantages of this versatility will be discussed later. BASIC EQUATIONS
FOR DESIGNING
THE FILTERS
The basis of designing the resistivity and the EM filters is the linear relationship between the kernels and the output functions. The output functions are the apparent resistivity (pa) and the mutual coupling ratios (Z/Z,) in the case of the resistivity and the EM sounding respectively. The expression for Schlumberger apparent resistivity, pas (s) is given by: Pas(S)
=
s*_I-hT(X) J,(hs)dh 0
where s is half the current electrode separation and the resistivity transform T(h) is given by T(h) = T,. The recurrence relation for T, is: Ti-1 =
Ti + pi-1 tanh(Xdi-1) 1 +Ti tanh (Adi- lpi-1
(pi and di being the resistivity and the thickness of the ith layer). Applying Hankel transformation (Fourier-Bessel integral (see e.g., Watson, 1962, section 14.3, eq. 3) to eq. 1 we have (Ghosh, 1971a):
299
T(h)
=
s
PUS J, (As) ds/s
(2)
0
With the change of variables given by: In(s) ln(l/h)
X=
Y=
(3)
eq. 2 assumes the following form: T(Y)
= s
-m
adX)
J, ‘&d-b+X)1
(4)
IdX
Eq. 4 is the convolution integral through which Ghosh (1971a, p. 195, eq. 9) established the linearity between T(y) and pas(x). Here, T(y) = resistivity transform, the output function; ,o~s(x) = Schlumberger apparent resistivity, the input function; and J,{exp[-(y -x)1} = the filter function. Expressions for the mutual coupling ratios (Z/Z,) for various coil configurations are (Ward, 1967; Koefoed et al., 1972): (a) for horizontal coplanar loops: (2/2,)1
= 1 - s3 7
A’R(h) J,(hs)dh
(5)
0
(b) for perpendicular loops: (Z/Zo)11
=
-s3 j=- h2R(h)J,(hs)dX
(6)
0 (c)
for vertical coplanar loops:
(Z/Z,)m = 1 - s2 7
XR(h)J,(hs)dh
(7)
0
and (d) for vertical coaxial loops: (2/2,)1,
i
1-Z
j&(h)J,(Xs)dh 1
-s
0
7
A2R(h)Jo(Xs)dh
(8)
0
where R (A) is the kernel function given by, R (A) = Ro,n(h), in which the first suffix represents the consideration of the impedance in air and the second suffix is the number of subsurface layers. The recurrence relation for Ro,n(h) is:
V(i-l),i + Ri,rz(X) hVW%)l R(i-l)‘n(h)
= 1 + V(i_l),i{Ri,n(h)[eXp(-2diVi)]
}
300
with R,,n(X)
= 0 and;
vj = (X2 + ,;)l’* VT
=
iWpoUi
vi, k = (vi
- vk)/(vi
+ vk)
and di being the conductivity and thickness of the ith layer). With the change of variables given in eq. 3 and rearranging, the above equations assume the following forms: (oi
7
{e-2Y[R(y)]
(Z/Z,)1
= 1 - ezx
(Z/Z&I
= -ezx
(Z/Z,)m
= 1 - J {R(y)} -cc
(Z/Z,,),,
= l-
-cc
7
--Lo
1
{R(Y)){-
[Jo(e(“-Y))]
{-e”-YIJl(e”-Y)]
{R(y)e-2Y}
}dy
{-e ‘(“-Y)[J1(ex-Y)]}dy
e3(x-Y)[J,,(ex-Y)]
}dy]
}dy
(9)
(19)
{-e2(x--Y)[ J, (eX-Y)] }dy
7 {R(y)} c -m
cc s
}{-e-(x--Y)
(11)
-
(12)
-cc
The above equations are convolution integrals. The term within the first set of brackets in the integrands above represents the input function, and that in the second set, the filter function, while various integrals are output functions, The filter functions in eq. 9 and eq. 10 were developed by Koefoed et al. (1972). VERSATILITY
OF THE FILTERING
OPERATIONS
In this section we shall demonstrate the versatile nature of various filtering operations involved in the computation of resistivity and EM sounding curves. It will be shown that by some simple mathematical manipulations we can use the electromagnetic filters to compute the resistivity sounding curves and the resistivity filters to compute the EM sounding curves. Use of EM filters in resistivity
sounding
Schlumberger sounding Once the linearity between T(y) and pas(x) was established through eq. 4,
301
reversing the concept of the input and the output functions, Ghosh (1971b) designed the filter for computing Schlumberger apparent resistivity curves. Combining eq. 1 and eq. 3 we write the Schlumberger apparent resistivity as: PUS
=
ezx_m I- e-YT(y)J,
(eX-Y) (-e-Y)dy
(13)
Rearranging:
J
pas(x) = ex
-m
e-YT(y){-eX-Y[J,(eX-Y)]
}dy
(14)
Comparing this convolution integral with eq. 10 we find that the term within .brackets in the integrand is the same filter function as that used in eq. 10 for the perpendicular loop system. This would mean that convolution of exp(-y)T(y) with this filter set (Koefoed et al., 1972, p.418, table 2) and multiplying the result by exp(x) could yield the Schlumberger apparent resistivity. Let us compare the filter function used in eq. 14 with the Schlumberger resistivity filter (Ghosh, 1971b). For this we rewrite eq. 13 in the following form : pas(x)
=
7
-m
5!‘(y){-e’(“-Y)[Jl(eX-y)]}dy
(15)
When T(y) is the input to the filter, the filter function is given by {-exp[ 2(x-y)] [J, exp(z-y )] } which has been derived by many workers (see introduction). Instead of this resistivity filter, if we use the EM filter as in equation (14), the input to the former [i.e., T(y)] changes to exp(x)exp(--y)?‘(y). This implies that the action of the resistivity filter can be derived from the EM filter by introducing appropriate changes in the input function.
Two-electrode sounding The expression for the two-electrode apparent resistivity, pat(s) is given by Kumar (1974): Pat(s) =
s
s
T(h) J,(hs)dh
(16)
0
where s = AM, the separatibn between the two active electrodes. Introducing the change of variables given by eq. 3, eq. 16 assumes the form: Pat(x) =
ex s -cc
W)Jo(e”-Y)(-e-Y)dy
302
or: pat(x)
=
s
~(y){-eX-YJo(eX-Y)}dy
(17)
-m
Comparing this convolution integral with eq. 9 we find that the term within the brackets in the integrand is the same filter function as that used in eq. 11 for the horizontal coplanar loop system. Convolution of T(y) with this filter set (Koefoed et al., 1972, p.415, table 1) would yield the apparent resistivities for the two-electrode system. In this case the action of the required resistivity filter can be derived from the EM filter without changing the input T(y). In other words, the resistivity filter for the two electrode system is the same as the EM filter for the coplanar horizontal loop system. Wenner sounding The expression for the Wenner apparent resistivity, paw(s) is given by: paw(s) = 2s j: T(h)[ J,(h)
- Jo
Ws)l dh
(18)
0
where s = a, the inter-electrode spacing. Let us write eq. 18 in the following form: Paw(s) = 2&Z&) - &(2s)
(19)
where:
P&(S)= s j”
T(h)Jo(hs)dX
(20)
= 2s 7 T(X)J, (X2s)dh
(21)
0
and: ,032s)
0
Comparing eq. 20 and eq. 21 with eq. 16, we can write eq. 19 as: Paw(s) = 2&t(s) - Pat@) Js’+Q
(22)
Eq. 22 suggests that the Wenner apparent resistivity paw(s) at s could be obtained from the two electrode apparent resistivity values at s and 2s computed by using the EM filter as mentioned in the preceding section. Using the above mentioned techniques, apparent resistivity curves for various electrode configurations have been computed. These are shown in Fig. 1 by the solid lines. Models for which these curves are computed are given in the figure. Blank circles on these curves are the apparent resistivity values from the tables given by Orellana and Mooney (1966) for the Schlumberger
303
-
USING
EM FILTERS
TABULATED
VALUES
TWO-ELECTRODE 121
Fig, 1. Apparent
resistivity
curves
for various
electrode
configurations
computed
by EM
filters.
and the Wenner system and by Mooney and Wetzel(1956) for the two-electrode system. The largest discrepancy between the computed and the tabulated values is found to be one point in the third decimal place. Use of resistivity
filters in EM sounding
Perpendicular loop system Let us consider eq. 10 for the perpendicular loop system which can be rearranged into the following form: (z/Z,),
= -ex
7
e-YR(y){-e2(X-Y)
Jl(ex-Y)}dy
(23)
_m
Comparing the above convolution integral with eq. 15 we find that the term within the brackets in the integrand is the same filter function as that used in eq. 16 for computing the Schlumberger apparent resistivities. The input to this filter is exp(x)exp(-y)R (y) instead of exp(2x)exp(--2y)R(y) as used in eq. 10. Introducing this modification in the input, (Z/Z,)11 can be computed with the help of the Schlumberger filter. Vertical coplanar loop system Comparing the convolution integral for this system given by eq. 11 with eq. 15 we find that the term within the brackets in the integrand of eq. 11 is the Schlumberger resistivity filter and its convolution with the EM kernel, R(y) would yield the output, the integral itself.
304
Vertical
coaxial loop system
Let us consider the dipole--dipole apparent resistivity, POD (Das and Ghosh, 1973, p.397, eq. A(6)): PaD(S) = s2 (1-b)
hT(h)J,(hs)dX-bs
j-
s h’T(h)J,,(As)dh 0
0
given by
1
where b is a constant depending upon the type of dipole-dipole configuration used. For radial dipole-dipole configuration b = 0.5 and the apparent resistivity for this system is: s2 par(s) = 2
[
m s hT(h)J,(hs)dh
-s s
h2T(h)Jo(hs)dX
0
0
(25)
3
For computing par(s) filters are available (Das and Ghosh, 1974, p. 772, table 1; Nyman and Landisman, 1977, p. 1043, table 2), which convolve with T(h) to yield par(s). Comparing eq. 25 with eq. 8 we observe that: (Z/Z,),,
= l-
par(s) I T(A)&(A)
Thus we can use a radial dipole--dipole filter which convolves with R (A) to yield (Z/Z,), . Using the above techniques we have computed the absolute values of the mutual coupling ratios (represented by the solid line) for various coil configurations for the models shown in Fig.?. The blank circles on the curves are
COIL
SEPARATION
= 1000
?,?c
MODELS
L+
‘.?
,3l
-
P
‘31 58
,/
t
\ \-
/
iz/z,!
\
PEiPENDlCULAR (I 1
32 3.3
L\
“ERTKPC COAXIAL
CC-
I
lo-5
m
T’ -4 IO
Fig. 2. Mutual
I 10
0
“SING
RESlSTlVlT\
“SlNG
EM FILTERS
L_
-3
coupling
IO
-2 ratios
FILTERS
I 16’
to”
I [SEC
for various
1
coil configurations
computed
by resistivity
filters.
305
the values computed by using the EM filters (Koefoed et al., 1972; Verma, 1973). The agreement between the solid lines and the corresponding blank circles is good enough except for the small differences at the lower periods in the case of the perpendicular loop system. These differences are due to the use of Ghosh’s filter. An extensive study of the frequency content of the resistivity kernels led Ghosh to choose a large sampling interval of (In lo)/3 which he considered to be the optimum from the computational speed and accuracy point of view. However, this large sampling resulted in a smaller frequency domain which can not accomodate the frequency content of all the electromagnetic sounding curves and thus deviations of the computed values with Ghosh’s filter from the actual values may occur in some cases having high frequency contents. In those cases, resistivity filters with smaller intervals (Das and Kumar; 1975; Johansen, 1975; Nyman and Landisman, 1977) could have been used. However, if the slight inaccuracy involved in the computation on the high frequency side is accepted, Ghosh’s filter could be used advantageously by a field party to compute an EM curve on programmable pocket calculators such as the HP-67 or TI-59. The results obtained in the preceding two sub-sections are summarized in Table I. The table clearly illustrates how various available resistivity and EMfilters can be used to compute responses for other systems. The need to redefine the input function for this purpose is also indicated in the last column. RELATIONSHIPS
Vertical
AMONG
VARIOUS
COIL CONFIGURATIONS
coaxial and vertical coplanar
loop systems
We have seen in the above analysis that the Schlumberger and the radial dipole-dipole filters could be used to convolve the EM kernel R (h) to compute the mutual coupling ratio for the vertical coplanar and the vertical coaxial loop systems, respectively. This leads us to study the relationship between these two systems analogous to the following relation that exists between the Schlumberger and the dipole-dipole apparent resistivities (Al’pin, 1950; Das and Ghosh, 1973): &D(s)
=
hS(S)
6Pas@) - beti6s
Differentiating eq. 7 with respect to s we have:
;
(Z/Zo),, = s
J hR(h)J,(hs)dX
+ s2
s
h’R(X)J,(hs)dX
(27)
0
0
Multiplying eq. 27 by (s/2) and adding it to eq. 7 we have:
(Z/Z,),,
+ ; ; (Z/Zo),, = l-
;
hR(h)J,(Xs)dX
--s j- h’R(h)J,(hs)dh 0
(28)
1
I
Jo + J,
Jl
Jo
Bessel function
Pas(S) =
-s j: hZT(h)JO(hs)dh] 0
0
$ [I hT( s.)J, fhs)dh
Vertical coplanar loops -J, (exp(=-Y) Irexp2(x-_y )_ Resistivity filter Radial dipole-dipole array
EM filters Perpendicular loops -J, (exp(n-y)}exp(x-y)
Resistivity filter Schlumberger array 4, {exp(x-y)]exp2(x-y)
EM filter Horizontal coplanar loops -JO {expk-y))exp(x-y)
Resistivity filter Wenner array --.[J,{exp(x-y)Tt--J,{2exp(x-y))l2exp(x-y)
Filters for different configurations
Summary of results
TABLE
array
Schlumberger
--s
h2R(h)J,(hs)dh
s0
Vertical coaxial loops
array
Schlumberger
Perpendicular loops Vertical coplanar loops
Two-electrode array Wenner array (at two different spacings)
Configuration using the same filter
1
no
n0
yes
no
yes
no no
Need to redefine input function
301
The right hand side of eq. 28 is equal to (Z/Z,),,. (-WG.I),,
Thus:
+ 2s -C 6s (W-O),
= (%%)I,
(29)
Eq. 29 is similar to eq. 26 revealing that: (1) when (Z/Z,)Irr approaches a constant value, (Z/Z,)Iv approaches (Z/Z,JID; (2) (Z/Z,),v passes through the minima and maxima of (Z/Z,)In as the gradient of the latter is zero at these points; and (3) for the rising part of (Z/Z,)III the values of (Z/Z,)Iv are greater than (Z/Z,,),, while in the descending part (Z/Z,)Iv is less than W-W,. COIL
0.0l-q ICY5 Fig. 3. Mutual earth model.
SEPARATION
1
/
1o-4 coupling
I
= 1000m
I,
I
1
10-3[SEC] 10* ratios
for vertical
I
I
/
lo-’ coaxial
103 and vertical
coplanar
loops
over the same
These statements are illustrated in Fig. 3. Through eq. 29 the mutual coupling ratio of the vertical coaxial loop system could be derived from that of the vertical coplanar loop system. This requires the computation of the gradient of the latter system at a point where this transformation is desired. Al’pin (1958) has given the formulae for computing the gradient of the Schlumberger curve in evaluating eq. 26. These formulae can be directly used to obtain the gradient of (Z/Z,), which would yield (Z/Z,)Iv when substituted in eq. 29. Horizontal
coplanar and vertical coplanar loop systems
Multiplying eq. 27 by s we obtain: hR(h)J1 0
(hs)dh + s3
s
0
h2R (X) J,(hs)dh
(39)
308
From eq. 5, 7 and 30 we have: s ;
(Z/Z,),
= l-
P-/w,,
+
1-
~mb),
or: (31) Thus (Z/Z& eq. 31.
can also be computed from (Z/Z,),,
and its gradient through
CONCLUSIONS
We have shown the versatility of resistivity and EM filters in computing the response of a layered earth. The concept could be used not only for the electrode and coil configurations mentioned in this paper but also for computing the components of the EM fields due to a number of other sourcereceiver configurations (Wait, 1972; Ward, 1967). In fact, computation of the components of EM fields due to a number of source-receiver configurations requires the evaluation of Hankel transform integrals which could be obtained by redefining appropriately the input functions and using the resistivity and the EM filters individually or jointly, as required. The versatility shown here should be useful for field geophysicists in computing approximate EM sounding curves on a programmable pocket calculator. However, for precise computations we do not recommend the use of a single filter. “For each purpose the most desirable filter must be sought, and in this search the possibilities implied in the definition of the input function may well be utilized. Not, however, without due consideration of the effect of the input function upon the sampling interval and the termination point of the filter” (Koefoed, personal communication, 1978). ACKNOWLEDGEMENT
The study was supported by UNDP under project No. IND/74/012.
309
NOTATION List of symbols and suffixes Symbol
Explanation
s
electrode spacing in Schlumberger, Wenner, two-electrode and dipoledipole configurations i.e. AB/2, AB/3, AM (two active electrodes) and OQ (centres between two dipoles) respectively. In EM sounding, it is the separation between the two coils (source and receiver) located at the surface of a horizontally stratified earth. variable of integration resistivity transform, an associated kernel EM kernel apparent resistivity Bessel function of order n ln (s) ln (l/h) mutual coupling ratio constant
A
T(A) R(h) Pa
Jfl X Y
(Z/Z, ) b
Suffixes S W t D ; II III IV
-
Schlumberger configuration Wenner configuration two-electrode configuration dipole-dipole configuration radial dipole-dipole configuration horizontal coplanar loop system perpendicular loop system vertical coplanar loop system vertical coaxial loop system
REFERENCES Al’pin, L.M., 1950. The theory of dipole sounding. Geostoptekhizdat, 1950. (Translation in: Dipole methods for measuring earth conductivity. Plenum Press, New York, N.Y., pp. l-60, 1966.) Al’pin, L.M., 1958. Transformation of sounding curves. Prikl. Geofiz., 19: 23-46. (Translation in: Dipole methods for measuring earth conductivity. Plenum Press, New York, N.Y., pp. 61-78, 1966.) Anderson, W.L., 1975. Improved digital filters for evaluating Fourier and Hankel transform integral. U.S. Geol. Surv., Rept., USGS-GD-75-012: 223 pp. Anderson, W.L., 1976. An optimal method for evaluating a class of convolution integrals with related kernels. U.S. Geol. Surv. Rept., USGS-GD-76-003. Das, U.C. and Ghosh, D.P., 1973. A study on the direct interpretation of dipole sounding resistivity measurements over layered earth. Geophys. Prospect., 21: 379-400. Das, U.C. and Ghosh, D.P., 1974. The determination of filter coefficients for the computation of standard curves for dipole resistivity sounding over layered earth by linear digital filtering. Geophys. Prospect., 22: 765-780. Das, U.C., Ghosh, D.P. and Biewinga, D.T., 1974. Transformation of dipole resistivity sounding measurements over layered earth by linear digital filtering. Geophys. Prospect., 22: 416-489.
310 Dean, W.C., 1958. Frequency analysis for gravity and magnetic interpretation. Geophysics, 23: 97-127. Ghosh, D.P., 1970. The application of linear filter theory to the direct interpretation of geoelectric resistivity measurements. Thesis, Technological University, Delft. Ghosh, D.P., 1971a. The application of linear filter theory to the direct interpretation of geoelectric resistivity sounding measurements. Geophys. Prospect., 19: 192-217. Ghosh, D.P., 1971b. 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 sounding. Geophys. Prospect., 23: 449-458. Johansen; H.K. and Srensen, K., 1977. Fast Hankel transforms. E.A.E.G., Zagreb. Koefoed, O., Ghosh, D.P. and Polman, G.J., 1972. Computation of type curves for electromagnetic depth sounding with a horizontal transmitting coil by means of a digital linear filter. Geophys. Prospect., 20: 406-420. Kumar, R., 1974. Direct interpretation of two-electrode resistivity soundings. Geophys. Prospect., 22: 224-237. Mooney, M.M. and Wetzel, W.W., 1956. The potentials about point electrode and apparent resistivity curves for a two, three, and four-layered earth. University of Minnesota Press, Minneapolis, Minn. Nyman, D.C. and Landisman, M., 1977. VES dipole-dipole filter coefficients. Geophysics, 42: 1037-1044. Orellana, E. and Mooney, H.M., 1966. Master tables and curves for vertical electric sounding over layered structures. Intercientia, Madrid. Robinson, E.A. and Treital, S., 1973. The Robinson and Trital Reader. Seismograph Service Corp., Tulsa, Okl., 3rd ed. Swarts, C.A., 1954. Some geometrical properties of residual maps. Geophysics, 15: 46-70. Verma, R.K., 1973. A feasibility study of electromagnetic depth sounding methods. Thesis, Technological University, Delft. Wait, J.R., 1970. Electromagnetic waves in stratified media. Macmillan, New York, N.Y., 2nd ed. Ward, S.H., 1967. Electromagnetic theory for geophysical methods. Mining Geophys., vol. 2, SEG, Tulsa, U.S.A. Watson, G.N., 1966. A treatise on the theory of Bessel functions. Cambridge University Press, Camb., 2nd ed.
297
THE VERSATILITY OF DIGITAL LINEAR FILTERS USED IN COMPUTING RESISTIVITY AND EM SOUNDING CURVES
U.C. DAS and SK. VERMA National Geophysical Research Institute, Hyderabad (India) (Received May 7, 1980; accepted October 5, 1980)
ABSTRACT Das, UC. and Verma, SK., 1981. The versatility of digital linear filters used in computing resistivity and EM sounding curves. Geoexploration, 18: 297-310. The applicability of EM filters to compute resistivity sounding curves and that of resistivity filters in computing EM sounding curves is shown. This has been possible mainly because of the flexible nature of the filtering process involved which permits us to redefine an input function in order to suit a particular filter function. The desired output can then be obtained by convolving the modified input with the chosen filter. This flexibility can be used advantageously to compute the electrical and EM response of a layered earth to a number of source-receiver configurations by using the resistivity and the EM filters individually or jointly, as required. The study leads us to establish the relationships among various coil configurations used in EM sounding. A great computational advantage accruing from the present study is that the nine point filter to compute Schlumberger resistivity curves can also be used to compute EM sounding curves for perpendicular and vertical coplanar loop systems on programmable pocket calculators such as the HP-67 or TI-59.
INTRODUCTION
The extensive use of electric filters in communization engineering in order to clean signals contaminated with noise is well known. In geophysical operations the same technique is utilized for recording the seismic signals. However, the application of electric filter theory to geophysical data processing and inte~re~tion was possible only after Swartz (1954) indicated the similarity between the electric filter and the transformation of geophysical quantities. It was then realized that the operators used to convert residual gravity maps into their derivatives, or to continue the residual maps upward or downward were in fact linear filtering operations. Dean (1958) put this thinking on a sound basis. The most ou~t~~ng applications of linear filter theory have been made in the field of seismic data processing and the credit goes jointly to Robinson and Treital (1973). However, in the field of resistivity interpretation, it is quite surprising that 0016-71~2~81/0000-OOOOf$O2.50
01981
Etsevier Scientific Publishing Company
298
although active research has been going on since as early as the 1930’s, it was as late as the 1970’s (Ghosh, 1970; 1971a; 1971b) that the principle of linear filtering was applied to the problem of direct interpretation of Schlumberger, Wenner and dipole-dipole sounding data (Ghosh, 1971a; Das et al., 1974) and later for the computation of theoretical resistivity and EM sounding curves (Ghosh, 1971b; Das and Ghosh, 1974; Kumar and Das, personal communication, 1975; Johansen, 1975; Nyman and Landisman, 1977; Koefoed et al., 1972). Expressions for various quantities such as apparent resistivity, EM fields or mutual impedance ratio etc., in d.c. and EM prospecting involve Bessel functions of orders 0 and 1. A suitable change of variables in these Hankel transforms results in convolution integrals whereby filters are developed to evaluate these integrals. Computing these Bessel function filters (J, and J,), not associated with any multiplying factor as in resistivity and EM filters referred earlier, at identical abscissae, Anderson (1979) has shown how to compute the related Hankel transforms in which the discrete values of the initial kernels are stored for repeated use in the evaluation of other integrals. Avoidance of re-evaluation of kernels thus saves substantial computer time. The availability of these Bessel function filters (Jo and J,) will lead to the computation of any Hankel transform. In this paper, however, we show a new aspect of the available resistivity and EM filters - their versatility. By simple manipulation of the kernels we can use EM filters to compute resistivity curves or resistivity filters to compute EM sounding curves. The advantages of this versatility will be discussed later. BASIC EQUATIONS
FOR DESIGNING
THE FILTERS
The basis of designing the resistivity and the EM filters is the linear relationship between the kernels and the output functions. The output functions are the apparent resistivity (pa) and the mutual coupling ratios (Z/Z,) in the case of the resistivity and the EM sounding respectively. The expression for Schlumberger apparent resistivity, pas (s) is given by: Pas(S)
=
s*_I-hT(X) J,(hs)dh 0
where s is half the current electrode separation and the resistivity transform T(h) is given by T(h) = T,. The recurrence relation for T, is: Ti-1 =
Ti + pi-1 tanh(Xdi-1) 1 +Ti tanh (Adi- lpi-1
(pi and di being the resistivity and the thickness of the ith layer). Applying Hankel transformation (Fourier-Bessel integral (see e.g., Watson, 1962, section 14.3, eq. 3) to eq. 1 we have (Ghosh, 1971a):
299
T(h)
=
s
PUS J, (As) ds/s
(2)
0
With the change of variables given by: In(s) ln(l/h)
X=
Y=
(3)
eq. 2 assumes the following form: T(Y)
= s
-m
adX)
J, ‘&d-b+X)1
(4)
IdX
Eq. 4 is the convolution integral through which Ghosh (1971a, p. 195, eq. 9) established the linearity between T(y) and pas(x). Here, T(y) = resistivity transform, the output function; ,o~s(x) = Schlumberger apparent resistivity, the input function; and J,{exp[-(y -x)1} = the filter function. Expressions for the mutual coupling ratios (Z/Z,) for various coil configurations are (Ward, 1967; Koefoed et al., 1972): (a) for horizontal coplanar loops: (2/2,)1
= 1 - s3 7
A’R(h) J,(hs)dh
(5)
0
(b) for perpendicular loops: (Z/Zo)11
=
-s3 j=- h2R(h)J,(hs)dX
(6)
0 (c)
for vertical coplanar loops:
(Z/Z,)m = 1 - s2 7
XR(h)J,(hs)dh
(7)
0
and (d) for vertical coaxial loops: (2/2,)1,
i
1-Z
j&(h)J,(Xs)dh 1
-s
0
7
A2R(h)Jo(Xs)dh
(8)
0
where R (A) is the kernel function given by, R (A) = Ro,n(h), in which the first suffix represents the consideration of the impedance in air and the second suffix is the number of subsurface layers. The recurrence relation for Ro,n(h) is:
V(i-l),i + Ri,rz(X) hVW%)l R(i-l)‘n(h)
= 1 + V(i_l),i{Ri,n(h)[eXp(-2diVi)]
}
300
with R,,n(X)
= 0 and;
vj = (X2 + ,;)l’* VT
=
iWpoUi
vi, k = (vi
- vk)/(vi
+ vk)
and di being the conductivity and thickness of the ith layer). With the change of variables given in eq. 3 and rearranging, the above equations assume the following forms: (oi
7
{e-2Y[R(y)]
(Z/Z,)1
= 1 - ezx
(Z/Z&I
= -ezx
(Z/Z,)m
= 1 - J {R(y)} -cc
(Z/Z,,),,
= l-
-cc
7
--Lo
1
{R(Y)){-
[Jo(e(“-Y))]
{-e”-YIJl(e”-Y)]
{R(y)e-2Y}
}dy
{-e ‘(“-Y)[J1(ex-Y)]}dy
e3(x-Y)[J,,(ex-Y)]
}dy]
}dy
(9)
(19)
{-e2(x--Y)[ J, (eX-Y)] }dy
7 {R(y)} c -m
cc s
}{-e-(x--Y)
(11)
-
(12)
-cc
The above equations are convolution integrals. The term within the first set of brackets in the integrands above represents the input function, and that in the second set, the filter function, while various integrals are output functions, The filter functions in eq. 9 and eq. 10 were developed by Koefoed et al. (1972). VERSATILITY
OF THE FILTERING
OPERATIONS
In this section we shall demonstrate the versatile nature of various filtering operations involved in the computation of resistivity and EM sounding curves. It will be shown that by some simple mathematical manipulations we can use the electromagnetic filters to compute the resistivity sounding curves and the resistivity filters to compute the EM sounding curves. Use of EM filters in resistivity
sounding
Schlumberger sounding Once the linearity between T(y) and pas(x) was established through eq. 4,
301
reversing the concept of the input and the output functions, Ghosh (1971b) designed the filter for computing Schlumberger apparent resistivity curves. Combining eq. 1 and eq. 3 we write the Schlumberger apparent resistivity as: PUS
=
ezx_m I- e-YT(y)J,
(eX-Y) (-e-Y)dy
(13)
Rearranging:
J
pas(x) = ex
-m
e-YT(y){-eX-Y[J,(eX-Y)]
}dy
(14)
Comparing this convolution integral with eq. 10 we find that the term within .brackets in the integrand is the same filter function as that used in eq. 10 for the perpendicular loop system. This would mean that convolution of exp(-y)T(y) with this filter set (Koefoed et al., 1972, p.418, table 2) and multiplying the result by exp(x) could yield the Schlumberger apparent resistivity. Let us compare the filter function used in eq. 14 with the Schlumberger resistivity filter (Ghosh, 1971b). For this we rewrite eq. 13 in the following form : pas(x)
=
7
-m
5!‘(y){-e’(“-Y)[Jl(eX-y)]}dy
(15)
When T(y) is the input to the filter, the filter function is given by {-exp[ 2(x-y)] [J, exp(z-y )] } which has been derived by many workers (see introduction). Instead of this resistivity filter, if we use the EM filter as in equation (14), the input to the former [i.e., T(y)] changes to exp(x)exp(--y)?‘(y). This implies that the action of the resistivity filter can be derived from the EM filter by introducing appropriate changes in the input function.
Two-electrode sounding The expression for the two-electrode apparent resistivity, pat(s) is given by Kumar (1974): Pat(s) =
s
s
T(h) J,(hs)dh
(16)
0
where s = AM, the separatibn between the two active electrodes. Introducing the change of variables given by eq. 3, eq. 16 assumes the form: Pat(x) =
ex s -cc
W)Jo(e”-Y)(-e-Y)dy
302
or: pat(x)
=
s
~(y){-eX-YJo(eX-Y)}dy
(17)
-m
Comparing this convolution integral with eq. 9 we find that the term within the brackets in the integrand is the same filter function as that used in eq. 11 for the horizontal coplanar loop system. Convolution of T(y) with this filter set (Koefoed et al., 1972, p.415, table 1) would yield the apparent resistivities for the two-electrode system. In this case the action of the required resistivity filter can be derived from the EM filter without changing the input T(y). In other words, the resistivity filter for the two electrode system is the same as the EM filter for the coplanar horizontal loop system. Wenner sounding The expression for the Wenner apparent resistivity, paw(s) is given by: paw(s) = 2s j: T(h)[ J,(h)
- Jo
Ws)l dh
(18)
0
where s = a, the inter-electrode spacing. Let us write eq. 18 in the following form: Paw(s) = 2&Z&) - &(2s)
(19)
where:
P&(S)= s j”
T(h)Jo(hs)dX
(20)
= 2s 7 T(X)J, (X2s)dh
(21)
0
and: ,032s)
0
Comparing eq. 20 and eq. 21 with eq. 16, we can write eq. 19 as: Paw(s) = 2&t(s) - Pat@) Js’+Q
(22)
Eq. 22 suggests that the Wenner apparent resistivity paw(s) at s could be obtained from the two electrode apparent resistivity values at s and 2s computed by using the EM filter as mentioned in the preceding section. Using the above mentioned techniques, apparent resistivity curves for various electrode configurations have been computed. These are shown in Fig. 1 by the solid lines. Models for which these curves are computed are given in the figure. Blank circles on these curves are the apparent resistivity values from the tables given by Orellana and Mooney (1966) for the Schlumberger
303
-
USING
EM FILTERS
TABULATED
VALUES
TWO-ELECTRODE 121
Fig, 1. Apparent
resistivity
curves
for various
electrode
configurations
computed
by EM
filters.
and the Wenner system and by Mooney and Wetzel(1956) for the two-electrode system. The largest discrepancy between the computed and the tabulated values is found to be one point in the third decimal place. Use of resistivity
filters in EM sounding
Perpendicular loop system Let us consider eq. 10 for the perpendicular loop system which can be rearranged into the following form: (z/Z,),
= -ex
7
e-YR(y){-e2(X-Y)
Jl(ex-Y)}dy
(23)
_m
Comparing the above convolution integral with eq. 15 we find that the term within the brackets in the integrand is the same filter function as that used in eq. 16 for computing the Schlumberger apparent resistivities. The input to this filter is exp(x)exp(-y)R (y) instead of exp(2x)exp(--2y)R(y) as used in eq. 10. Introducing this modification in the input, (Z/Z,)11 can be computed with the help of the Schlumberger filter. Vertical coplanar loop system Comparing the convolution integral for this system given by eq. 11 with eq. 15 we find that the term within the brackets in the integrand of eq. 11 is the Schlumberger resistivity filter and its convolution with the EM kernel, R(y) would yield the output, the integral itself.
304
Vertical
coaxial loop system
Let us consider the dipole--dipole apparent resistivity, POD (Das and Ghosh, 1973, p.397, eq. A(6)): PaD(S) = s2 (1-b)
hT(h)J,(hs)dX-bs
j-
s h’T(h)J,,(As)dh 0
0
given by
1
where b is a constant depending upon the type of dipole-dipole configuration used. For radial dipole-dipole configuration b = 0.5 and the apparent resistivity for this system is: s2 par(s) = 2
[
m s hT(h)J,(hs)dh
-s s
h2T(h)Jo(hs)dX
0
0
(25)
3
For computing par(s) filters are available (Das and Ghosh, 1974, p. 772, table 1; Nyman and Landisman, 1977, p. 1043, table 2), which convolve with T(h) to yield par(s). Comparing eq. 25 with eq. 8 we observe that: (Z/Z,),,
= l-
par(s) I T(A)&(A)
Thus we can use a radial dipole--dipole filter which convolves with R (A) to yield (Z/Z,), . Using the above techniques we have computed the absolute values of the mutual coupling ratios (represented by the solid line) for various coil configurations for the models shown in Fig.?. The blank circles on the curves are
COIL
SEPARATION
= 1000
?,?c
MODELS
L+
‘.?
,3l
-
P
‘31 58
,/
t
\ \-
/
iz/z,!
\
PEiPENDlCULAR (I 1
32 3.3
L\
“ERTKPC COAXIAL
CC-
I
lo-5
m
T’ -4 IO
Fig. 2. Mutual
I 10
0
“SING
RESlSTlVlT\
“SlNG
EM FILTERS
L_
-3
coupling
IO
-2 ratios
FILTERS
I 16’
to”
I [SEC
for various
1
coil configurations
computed
by resistivity
filters.
305
the values computed by using the EM filters (Koefoed et al., 1972; Verma, 1973). The agreement between the solid lines and the corresponding blank circles is good enough except for the small differences at the lower periods in the case of the perpendicular loop system. These differences are due to the use of Ghosh’s filter. An extensive study of the frequency content of the resistivity kernels led Ghosh to choose a large sampling interval of (In lo)/3 which he considered to be the optimum from the computational speed and accuracy point of view. However, this large sampling resulted in a smaller frequency domain which can not accomodate the frequency content of all the electromagnetic sounding curves and thus deviations of the computed values with Ghosh’s filter from the actual values may occur in some cases having high frequency contents. In those cases, resistivity filters with smaller intervals (Das and Kumar; 1975; Johansen, 1975; Nyman and Landisman, 1977) could have been used. However, if the slight inaccuracy involved in the computation on the high frequency side is accepted, Ghosh’s filter could be used advantageously by a field party to compute an EM curve on programmable pocket calculators such as the HP-67 or TI-59. The results obtained in the preceding two sub-sections are summarized in Table I. The table clearly illustrates how various available resistivity and EMfilters can be used to compute responses for other systems. The need to redefine the input function for this purpose is also indicated in the last column. RELATIONSHIPS
Vertical
AMONG
VARIOUS
COIL CONFIGURATIONS
coaxial and vertical coplanar
loop systems
We have seen in the above analysis that the Schlumberger and the radial dipole-dipole filters could be used to convolve the EM kernel R (h) to compute the mutual coupling ratio for the vertical coplanar and the vertical coaxial loop systems, respectively. This leads us to study the relationship between these two systems analogous to the following relation that exists between the Schlumberger and the dipole-dipole apparent resistivities (Al’pin, 1950; Das and Ghosh, 1973): &D(s)
=
hS(S)
6Pas@) - beti6s
Differentiating eq. 7 with respect to s we have:
;
(Z/Zo),, = s
J hR(h)J,(hs)dX
+ s2
s
h’R(X)J,(hs)dX
(27)
0
0
Multiplying eq. 27 by (s/2) and adding it to eq. 7 we have:
(Z/Z,),,
+ ; ; (Z/Zo),, = l-
;
hR(h)J,(Xs)dX
--s j- h’R(h)J,(hs)dh 0
(28)
1
I
Jo + J,
Jl
Jo
Bessel function
Pas(S) =
-s j: hZT(h)JO(hs)dh] 0
0
$ [I hT( s.)J, fhs)dh
Vertical coplanar loops -J, (exp(=-Y) Irexp2(x-_y )_ Resistivity filter Radial dipole-dipole array
EM filters Perpendicular loops -J, (exp(n-y)}exp(x-y)
Resistivity filter Schlumberger array 4, {exp(x-y)]exp2(x-y)
EM filter Horizontal coplanar loops -JO {expk-y))exp(x-y)
Resistivity filter Wenner array --.[J,{exp(x-y)Tt--J,{2exp(x-y))l2exp(x-y)
Filters for different configurations
Summary of results
TABLE
array
Schlumberger
--s
h2R(h)J,(hs)dh
s0
Vertical coaxial loops
array
Schlumberger
Perpendicular loops Vertical coplanar loops
Two-electrode array Wenner array (at two different spacings)
Configuration using the same filter
1
no
n0
yes
no
yes
no no
Need to redefine input function
301
The right hand side of eq. 28 is equal to (Z/Z,),,. (-WG.I),,
Thus:
+ 2s -C 6s (W-O),
= (%%)I,
(29)
Eq. 29 is similar to eq. 26 revealing that: (1) when (Z/Z,)Irr approaches a constant value, (Z/Z,)Iv approaches (Z/Z,JID; (2) (Z/Z,),v passes through the minima and maxima of (Z/Z,)In as the gradient of the latter is zero at these points; and (3) for the rising part of (Z/Z,)III the values of (Z/Z,)Iv are greater than (Z/Z,,),, while in the descending part (Z/Z,)Iv is less than W-W,. COIL
0.0l-q ICY5 Fig. 3. Mutual earth model.
SEPARATION
1
/
1o-4 coupling
I
= 1000m
I,
I
1
10-3[SEC] 10* ratios
for vertical
I
I
/
lo-’ coaxial
103 and vertical
coplanar
loops
over the same
These statements are illustrated in Fig. 3. Through eq. 29 the mutual coupling ratio of the vertical coaxial loop system could be derived from that of the vertical coplanar loop system. This requires the computation of the gradient of the latter system at a point where this transformation is desired. Al’pin (1958) has given the formulae for computing the gradient of the Schlumberger curve in evaluating eq. 26. These formulae can be directly used to obtain the gradient of (Z/Z,), which would yield (Z/Z,)Iv when substituted in eq. 29. Horizontal
coplanar and vertical coplanar loop systems
Multiplying eq. 27 by s we obtain: hR(h)J1 0
(hs)dh + s3
s
0
h2R (X) J,(hs)dh
(39)
308
From eq. 5, 7 and 30 we have: s ;
(Z/Z,),
= l-
P-/w,,
+
1-
~mb),
or: (31) Thus (Z/Z& eq. 31.
can also be computed from (Z/Z,),,
and its gradient through
CONCLUSIONS
We have shown the versatility of resistivity and EM filters in computing the response of a layered earth. The concept could be used not only for the electrode and coil configurations mentioned in this paper but also for computing the components of the EM fields due to a number of other sourcereceiver configurations (Wait, 1972; Ward, 1967). In fact, computation of the components of EM fields due to a number of source-receiver configurations requires the evaluation of Hankel transform integrals which could be obtained by redefining appropriately the input functions and using the resistivity and the EM filters individually or jointly, as required. The versatility shown here should be useful for field geophysicists in computing approximate EM sounding curves on a programmable pocket calculator. However, for precise computations we do not recommend the use of a single filter. “For each purpose the most desirable filter must be sought, and in this search the possibilities implied in the definition of the input function may well be utilized. Not, however, without due consideration of the effect of the input function upon the sampling interval and the termination point of the filter” (Koefoed, personal communication, 1978). ACKNOWLEDGEMENT
The study was supported by UNDP under project No. IND/74/012.
309
NOTATION List of symbols and suffixes Symbol
Explanation
s
electrode spacing in Schlumberger, Wenner, two-electrode and dipoledipole configurations i.e. AB/2, AB/3, AM (two active electrodes) and OQ (centres between two dipoles) respectively. In EM sounding, it is the separation between the two coils (source and receiver) located at the surface of a horizontally stratified earth. variable of integration resistivity transform, an associated kernel EM kernel apparent resistivity Bessel function of order n ln (s) ln (l/h) mutual coupling ratio constant
A
T(A) R(h) Pa
Jfl X Y
(Z/Z, ) b
Suffixes S W t D ; II III IV
-
Schlumberger configuration Wenner configuration two-electrode configuration dipole-dipole configuration radial dipole-dipole configuration horizontal coplanar loop system perpendicular loop system vertical coplanar loop system vertical coaxial loop system
REFERENCES Al’pin, L.M., 1950. The theory of dipole sounding. Geostoptekhizdat, 1950. (Translation in: Dipole methods for measuring earth conductivity. Plenum Press, New York, N.Y., pp. l-60, 1966.) Al’pin, L.M., 1958. Transformation of sounding curves. Prikl. Geofiz., 19: 23-46. (Translation in: Dipole methods for measuring earth conductivity. Plenum Press, New York, N.Y., pp. 61-78, 1966.) Anderson, W.L., 1975. Improved digital filters for evaluating Fourier and Hankel transform integral. U.S. Geol. Surv., Rept., USGS-GD-75-012: 223 pp. Anderson, W.L., 1976. An optimal method for evaluating a class of convolution integrals with related kernels. U.S. Geol. Surv. Rept., USGS-GD-76-003. Das, U.C. and Ghosh, D.P., 1973. A study on the direct interpretation of dipole sounding resistivity measurements over layered earth. Geophys. Prospect., 21: 379-400. Das, U.C. and Ghosh, D.P., 1974. The determination of filter coefficients for the computation of standard curves for dipole resistivity sounding over layered earth by linear digital filtering. Geophys. Prospect., 22: 765-780. Das, U.C., Ghosh, D.P. and Biewinga, D.T., 1974. Transformation of dipole resistivity sounding measurements over layered earth by linear digital filtering. Geophys. Prospect., 22: 416-489.
310 Dean, W.C., 1958. Frequency analysis for gravity and magnetic interpretation. Geophysics, 23: 97-127. Ghosh, D.P., 1970. The application of linear filter theory to the direct interpretation of geoelectric resistivity measurements. Thesis, Technological University, Delft. Ghosh, D.P., 1971a. The application of linear filter theory to the direct interpretation of geoelectric resistivity sounding measurements. Geophys. Prospect., 19: 192-217. Ghosh, D.P., 1971b. 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 sounding. Geophys. Prospect., 23: 449-458. Johansen; H.K. and Srensen, K., 1977. Fast Hankel transforms. E.A.E.G., Zagreb. Koefoed, O., Ghosh, D.P. and Polman, G.J., 1972. Computation of type curves for electromagnetic depth sounding with a horizontal transmitting coil by means of a digital linear filter. Geophys. Prospect., 20: 406-420. Kumar, R., 1974. Direct interpretation of two-electrode resistivity soundings. Geophys. Prospect., 22: 224-237. Mooney, M.M. and Wetzel, W.W., 1956. The potentials about point electrode and apparent resistivity curves for a two, three, and four-layered earth. University of Minnesota Press, Minneapolis, Minn. Nyman, D.C. and Landisman, M., 1977. VES dipole-dipole filter coefficients. Geophysics, 42: 1037-1044. Orellana, E. and Mooney, H.M., 1966. Master tables and curves for vertical electric sounding over layered structures. Intercientia, Madrid. Robinson, E.A. and Treital, S., 1973. The Robinson and Trital Reader. Seismograph Service Corp., Tulsa, Okl., 3rd ed. Swarts, C.A., 1954. Some geometrical properties of residual maps. Geophysics, 15: 46-70. Verma, R.K., 1973. A feasibility study of electromagnetic depth sounding methods. Thesis, Technological University, Delft. Wait, J.R., 1970. Electromagnetic waves in stratified media. Macmillan, New York, N.Y., 2nd ed. Ward, S.H., 1967. Electromagnetic theory for geophysical methods. Mining Geophys., vol. 2, SEG, Tulsa, U.S.A. Watson, G.N., 1966. A treatise on the theory of Bessel functions. Cambridge University Press, Camb., 2nd ed.