Journal of Applied Geophysics 40 Ž1998. 179–186
Walsh spectra of gravity anomalies over some simple sources R.K. Shaw a , B.N.P. Agarwal a
a,)
, B.K. Nandi
b
Department of Applied Geophysics, Indian School of Mines, Dhanbad 826 004, India b Marine Wing, Geological SurÕey of India, Calcutta, India Received 2 December 1997; accepted 24 June 1998
Abstract Walsh transform of gravity anomalies over a point mass, a horizontal and a vertical line mass have been computed to obtain a cyclic shift invariant differential energy density ŽDED. function. Quantitative relations between DED spectral characteristics with depth to centroidrtop of the source have been established. The effects of profile length, sampling interval, random noise and zero padding have been investigated. Applicability of the proposed method has been evaluated through two field examples. q 1998 Elsevier Science B.V. All rights reserved. Keywords: Walsh transform; Gravity anomaly; DED spectrum; Sequency; Shift invariant
1. Introduction Analysis of an observed potential field anomaly in terms of its spectral components has become popular in recent years. The rationale of this approach lies in the fact that a composite anomaly satisfying certain conditions can be decomposed into constituent components with energy distribution function determined by the source parameters. Very often, a kernel of sinusoidal function is used for spectral decomposition. However, it is well known that any complete set of orthogonal functions, e.g., Walsh functions ŽWalsh, 1923. can also be used. In sharp contrast with sinusoids which have continuous values, the Walsh functions possess bistable states having values either q1 or y1. Sinusoids are characterised by their frequencies whereas Walsh functions are represented by their sequencies. Sequency is defined as half the average number of zero crossing per unit time and thus forms a generalized concept of frequency Ž Harmuth, 1972. for quasi-periodic functions. Fig. 1 shows first eight Walsh functions with corresponding sinusoids to compare their similarity and differences ŽAhmed and Rao, 1975.. Although sinusoids are always ordered by their frequencies, Walsh functions possess different orderings depending upon their mode of generation, e.g., sequency, Hadamard and Paley orderings ŽYuen, 1972. . Each ordering is associated with at least one spectral mode. Some
)
Corresponding author. Fax: q91-326-203042, q91-326-202380
0926-9851r98r$ - see front matter q 1998 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 6 - 9 8 5 1 Ž 9 8 . 0 0 0 2 7 - 5
180
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
Fig. 1. First eight sequency ordered Walsh functions and corresponding sinusoids; cal and sal represent cosine and sine like Walsh functions respectively Žafter Ahmed and Rao, 1975..
common modes ŽAhmed and Rao, 1975; Rao et al., 1973. are: Ž a. Sequency ordered Walsh spectrum Žcyclic shift variant.; Žb. Hadamard ordered Walsh spectrum Žcyclic shift invariant.; Žc. Quadratic or Q-spectrum Žcyclic shift invariant. ; Ž d. Optimum quadratic or J-spectrum Ž cyclic shift invariant. . Out of these, Hadamard ordered Walsh spectrum consists of a subset of Q-spectrum. Computation of J-spectrum is very involved procedure and has not yet been applied in analyzing geophysical data. Some applications of Walsh transform in geophysics are in Ža. signal to noise ratio enhancement ŽGubbins et al., 1971. , Ž b. pattern recognition Ž Bath and Burman, 1972. , Ž c. feature extraction Ž Chen, 1972; Chen and Boucher, 1973. , Ž d. data compression Ž Bois, 1972; Wood, 1974. , Ž e. identification of boundaries from well log data ŽLanning and Johnson, 1983. , Ž f. gravity data interpretation Ž Shaw and Agarwal, 1990. , Žg. resistivity filtering Ž Pal, 1991. and Žh. density mapping Ž Keating, 1992. . Shaw and Agarwal Ž1990. have evolved an interpretation scheme called Sequency Octave Analysis Ž SOA. based on characteristics of sequency ordered Walsh energy spectra of gravity anomalies over a few simple geometrical sources. However, Walsh energy spectrum depends upon the position of origin which requires a priori information on horizontal location of the source centroid. In this work, we have utilized Q-spectrum—termed as Walsh spectrum, to interpret gravity anomalies over some idealised sources, viz., a point mass, a horizontal line mass and a semi-infinite vertical line mass.
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
181
2. Theory A modified fast Walsh Transform ŽMFWT. of an N-periodic Ž N s 2 n . data sequence, x i 4 is obtained from ŽAhmed and Rao, 1975; Beauchamp, 1975. g Ž n . s Ž 1rN . d Ž n . x Ž n . , Ž1. TŽ . T where g n s w g Ž0., g Ž1., g Ž2., . . . , g Ž N y 1.x, x Ž n. s w x Ž0., x Ž1., x Ž2., . . . , x Ž N y 1.x and T stands for transpose of a matrix. The transformation matrix dŽ n. is obtained from the recurrence relation d Ž k q 1. s
dŽ k . 2
k r2
dŽ k .
I Ž k . y2 k r2 I Ž k .
with dŽ0. s 1 and I Ž k . represents a 2 k = 2 k identity matrix. Following Ahmed and Rao Ž1975, p. 136., the Walsh spectrum is defined as w0 Ž0. s g 2 Ž0. , ks2 my1
wm Ž 0 . s
Ý
g2Žk.,
m s 1,2,3, . . . ,n
Ž2.
ks2 my1 m
y1yl my1 Ž m . and wmŽ l . s Ý ks2 g Ž k . g Ž l q k . y Ýly1 q j . with m s 1,2, . . . , n and ks2 my1 js0 g 2 q j y l g Ž 2 my1 l s 1,2, . . . 2 . On the basis of the above formulation, Walsh spectrum can be computed for Ž n q 1. independent spectral groups of different sizes with the largest one having Nr4 q 1 independent spectral points. Each spectral point in a group has contribution from constituent sequencies and has been identified by a number, l called the sequency number. As a word of caution, we note that the sequency number has no physical interpretation in terms of sequency and should not be confused with wave number in Fourier transform. Further, this Walsh spectrum cannot be strictly called an energy spectrum as it may assume negative values Ž Eq. Ž 2.. . The above concept differs from the Fourier energy spectrum in
Fig. 2. Differential energy density ŽDED. of Walsh spectra of gravity anomalies over a horizontal line mass at indicated depths.
182
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
Fig. 3. Variation of l max with depth for different sources.
which a spectral point represents contribution to energy from a particular frequency only. The group of Walsh spectrum having maximum spectral points has been used in the present study.
3. Methodology In absence of any analytical solution for Walsh spectrum, the present problem has been solved numerically. Walsh spectra for gravity field over idealised sources, viz., a point mass, a horizontal line mass and a semi infinite vertical line mass have been computed and a few characteristic features are correlated with depth to centroid or top of the source to evolve an interpretation scheme. The effects of profile length, sampling interval, noisy data and zero padding have been investigated. 3.1. Simulated examples The gravity anomalies over a horizontal line mass with depth 3, 5 and 7 km have been sampled every 1 km over a profile length of 128 km. The Walsh spectra for these data have been computed from Eq. Ž2.. The differential energy density ŽDED. of these spectra, computed as the difference in Walsh energy between two consecutive sequency numbers exhibits a peak at l max ŽFig. 2. . It is evident that an increase in depth shifts the peak towards increasing sequency number. To establish a Table 1 Effect of profile length: the gravity anomaly over a horizontal line mass at a depth of 5 km is sampled every 1 km over profiles of different length Profile lengh Žkm.
l max
Computed depth Žkm.
64 128 256 512
6.20 5.82 5.79 5.79
5.25 4.93 4.90 4.90
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
183
Table 2 Effect of sampling interval: the gravity anomaly over a horizontal line mass at a depth of 8 km has been computed with varying sampling intervals over a profile 128 km long Sampling interval Žkm.
l max
Computed depth Žkm.
0.5 1.0 2.0 4.0
19.05 9.53 4.78 2.39
8.06 8.07 8.10 8.10
quantitative relationship between depth and l max , gravity anomalies over the selected models have been analysed for depth varying from 3 to 10 km at an interval of 0.2 km Ž Fig. 3. . The following relations have evolved: Depthrdata spacing
s Ž 1.090 " 0.003 . l max ,
for point mass,
s Ž 0.847 " 0.015 . l max ,
for horizontal line mass,
2 s Ž 0.615 " 0.001 . l max q Ž y0.0096 " 0.0001 . l max for semi infinite vertical line mass.
3.2. Effects of profile length and sampling interÕal The effect of change in profile length has been studied by analysing gravity anomaly over a horizontal line mass at a depth of 5 km and profile lengths 64, 128, 256 and 512 km. The results are summarised in Table 1. Further, the effect of sampling interval has been investigated from gravity anomalies over a horizontal line mass at a depth of 8 km, fixed profile length 128 km, and selecting the sampling intervals 0.5, 1, 2 and 4 km. The results are given in Table 2. This study, extended to include all the models, suggests that a profile length of 10–12 times the depth provides accurate results. Further, the depth to sampling interval ratio should lie between 2.5–10.
Fig. 4. Ža. Effect of random noise. The numbers represent percentage of noise in terms of peak anomaly used to corrupt the signal and Žb. effect of zero padding to profile ends. Numbers indicate actual profile length.
184
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
Table 3 Effect of random noise: gravity anomaly over a horizontal line mass at a depth of 5 km computed at 1 km sampling interval over a 128 km long profile has been mixed with various percentage of random noise in terms of peak anomaly magnitude Percentage of random noise
l max
Computed depth Žkm.
0 2 5 10
5.82 5.78 5.62 5.48
4.93 4.90 4.76 4.64
3.3. Effects of random noise and zero padding To investigate the effect of noise on the proposed technique, the gravity anomaly over a horizontal line mass at depth of 5 km has been mixed with a Gaussian random noise of amplitude 2, 5 and 10% of the peak anomaly. The results of this analysis Ž Fig. 4a and Table 3. reveal that the present scheme can tolerate high noise level. In order to satisfy the computational requirement of MFWT in form of a positive integer power of 2, the simulated data points have been padded by zeroes. Gravity anomalies over a horizontal line mass at 5.0 km depth have been computed at interval of 1 km over profile lengths 68, 88 and 108 km. These data have been padded by requisite number of zeroes to extend the profile to 128 km. DED spectra Ž Fig. 4b. show that even substantial padding of zeroes have insignificant effect on the source depth estimation by the present method.
4. Field examples 4.1. Cuba chromite anomaly The gravity anomaly Ž Fig. 5a. over a chromite ore deposit Ž Robinson and Coruh, 1988. discretised at an interval of 2.15 m by bicubic interpolation scheme ŽBhattacharyya, 1969. padded by zeroes to
Fig. 5. Ža. Residual gravity anomaly over a chromite ore Žafter Robinson and Coruh, 1988., Žb. corresponding DED spectrum and Žc. Fourier energy spectrum. Frequency number represents the harmonics.
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
185
Fig. 6. Ža. Gravity anomaly over Louga area Žafter Nettleton, 1976., Žb. its DED spectrum and Žc. Fourier energy spectrum. Frequency number represents the harmonics.
obtain 64 points. The DED spectrum Ž Fig. 5b. of this anomaly shows a maximum at l max s 9.22 corresponding to a depth of 21.6 m for a spherical source. This agrees well with the drill hole depth of 21.0 m. For comparison, the Fourier energy spectrum of this anomaly has been shown in Fig. 5c which yields a depth of 19.6 m. The depth determination from Fourier energy spectrum involves computation of spectral decay rate which is influenced by personal judgement. However, DED spectrum is considerably smooth and hence determination of l max is less susceptible to error. 4.2. Louga anomaly A gravity anomaly profile ŽFig. 6a. over Louga area ŽNettleton, 1976. has been discretised at an interval of 1.1 km. Its DED spectrum Ž Fig. 6b. exhibits a peak at l max s 7.80 to yield a depth of 9.35 km for a spherical source which agrees well with the results obtained by various authors Ž Shaw and Agarwal, 1990. . Fig. 6c shows the corresponding Fourier energy spectrum which provides a depth of 7.70 km.
5. Conclusion A cyclic shift invariant spectral mode of Walsh transform has been used to interpret gravity anomaly over idealised sources. This spectrum gets divided into smaller subgroups of different sizes.The group with maximum number of spectral points has been used to evolve an interpretation scheme based on differential energy density. The peak positions of DED spectra are related to depths of idealised sources. Studies have shown that a profile length of 10–12 times the depth of the source is sufficient for reliable depth estimation. Sampling interval to depth ratio can be in the range 2.5–10. A random noise of about 10% of peak anomaly does not affect depth computations considerably. Padding large number of zeroes to the end of profiles has no appreciable effect on the scheme. Analyses of two field examples demonstrate the applicability of the technique.
186
R.K. Shaw et al.r Journal of Applied Geophysics 40 (1998) 179–186
Acknowledgements The authors are thankful to Dr. P. Keating and Prof. P.C. Pal for their constructive suggestions leading to substantial improvement of the manuscript. Thanks are due to the Director and to the Head of the Department of Applied Geophysics, Indian School of Mines for providing necessary facilities. BKN is thankful to the Council of Scientific and Industrial Research, New Delhi for financial support.
References Ahmed, N., Rao, K.R., 1975, Orthogonal Transforms for Digital Signal Processing. Springer-Verlag, 263 pp. Bath, M., Burman, S., 1972. Walsh spectroscopy of Rayleigh waves caused by underground explosions, Proc. Symp. Applications of Walsh Functions, Washington, DC, pp. 48–63. Beauchamp, K.G., 1975. Walsh Functions and Their Applications. Academic Press, 236 pp. Bhattacharyya, B.K., 1969. Bicubic spline as a method for treatment of potential field data. Geophysics 34, 402–423. Bois, P., 1972. Analyse sequentille. Geophysical Prospecting 20, 497–513. Chen, C., 1972. Walsh domain processing of marine seismic data, Proc. Symp. Applications of Walsh Functions, Washington, DC, pp. 64–67. Chen, C.H., Boucher, R.E., 1973. Further results on Walsh domain processing of marine seismic data, Proc. Symp. Applications of Walsh Functions, Washington, DC, pp. 253–256. Gubbins, D., Scollar, I., Wisskirchen, P., 1971. Two dimensional digital filtering with Haar and Walsh transforms. Annales de Geophysique 27, 85–104. Harmuth, H.F., 1972. Transmission of Information by Orthogonal Functions. Springer-Verlag, 325 pp. Keating, P., 1992. Density mapping from gravity data using the Walsh transform. Geophysics 57, 637–642. Lanning, E.N., Johnson, D.M., 1983. Automated identification of rock boundaries: an application of Walsh transform to geophysical well log analysis. Geophysics 48, 197–205. Nettleton, L.L., 1976. Gravity and Magnetics in Oil Prospecting. McGraw-Hill Book, 464 pp. Pal, P.C., 1991. A Walsh sequency filtration method for integrating the resistivity log and sounding data. Geophysics 56, 1259–1266. Rao, K.R., Abdussattar, A.L., Ahmed, N., 1973. On cyclic autocorrelation and Walsh Hadamard Transform. IEEE Trans. Electromagnetic Compatibility, EMC 15, 141–146. Robinson, E.S., Coruh, C., 1988. Basic Exploration Geophysics, Wiley. Shaw, R.K., Agarwal, B.N.P., 1990. The application of Walsh transforms to interpret gravity anomalies due to some simple geometrically shaped causative sources: a feasibility study. Geophysics 55, 843–850. Walsh, J.L., 1923. A closed set of orthogonal functions. American Journal of Mathematics 45, 5–22. Wood, L.C., 1974. Seismic data compression methods. Geophysics 39, 499–525. Yuen, C.K., 1972. Remarks on ordering of the Walsh functions. IEEE Trans. on Computers C 21, 1452.