Computers & Geosciences Vol. 13. No. 5. pp. 549-560. 1987 0098-3004 87 $3.00 + 0.00 Printed in Great Bntmn. All rights reserved Copyright ~ 1987 Pergamon Journals Ltd SHORT NOTE A BASIC PROGRAM FOR 2-D SPECTRAL ANALYSIS OF GRAVITY DATA AND SOURCE-DEPTH ESTIMATION K. DIMITRIADI$ I, G . - A TSELENTIS t and K. THANASSOULAS 2 tAthens University, Geophysics Department. Ilisia. Athens 15701 and 2Institute of Geological Sciences. Mesoghion 70 Avenue. Athens. Greece (Received 23 December 1986; accepted 8 ,tfa), 1987) The spectral methods of analysis have been employed The Fourier transform of these data results in a set increasingly in recent years. In these methods, the of real XR and imaginary Xt amplitudes by which the characteristics of the observed anomalies are studied field values given at the grid points (x,y) can be by first transforming the data from the space to the represented by the sum: frequency domain and then analyzing their frequency characteristics (Bath, 1974). g(x,y) = ~ Z, X~ cos [(2n/DX.N)(kx + my)] The amplitude and phase relationships among the various frequencies has been used extensively by + Xff sin [(2n/DX'N)(kx + my) l many workers for the interpretation of gravity data, where DX is the grid interval. particularly in the situation of downward continu- Equation (2) can be written as follows: ation and source-depth estimation (Spector and Grant, 1970;TreiteI, Clement. and Kaui, 1971;Green, g(x,y) = ~Z C~'cos [(2n/DX.N)(kx + my) 1972; Hahn, Kind, and Mishra, 1976; Pal. Khurana, * " and Unnikrishnan, 1979; Negi, Agrawal, and Rao. - P~] (3) 1983; Bose and Sengupta, 1984; Tselentis, Dimitri- adis, and Drakopoulos, 1986). where P is the appropriate phase angle, and Following the approaches of Battacharya (1966) C~ = [(XsX) 2 + (X'f)2] ~/2. (4) and Treitel, Clement, and Kaul (1971), the power spectrum, when amplitude is on a logarithmic scale It is obvious that each C is the amplitude of a versus a linear scale for the frequency, may show partial field wave with wavelength DX'N/(k 2 + m 2) and frequency F = (k2 + m2). frequency intervals where the logarithms of the am- plitudes may be represented by a linear function of In order to calculate the radial spectrum for each frequency, with amplitudes decreasing with increasing data set we start first by calculating the 2-D power frequencies. The slope of the straightline is propor- spectrum: tional to the depth to the top of the body. Thus, if k SP(L J) = [XR(L j)2 + XI(L J)~'] denotes the wavenumber and S(k) the power spec- trum, the depth d to the source can be estimated from where XR(L J) is the real part and XI(L J) the imagin- the relation S(k) = f(k), by employing the formula: ary part of the data set at the point LJ. The radial spectrum is calculated by superposing In S(k) = - 2kd. (I) the 2-D spectrum by a number of omocentric rings It is obvious that the same approach can be with center the point (I,I) (upper left point of the followed for the situation of 2-dimensional data by matrix SP) which is the lowest frequency component computing the radial spectrum of all the particular of the data set (mean val.), and with radial distances waves falling within a certain frequency range as ex- 0.O-0.5:(wavenumber = 0.0 cycl/grid interval) plained next. 0.5-1.5:(wavenumber = I.O/(N.DX) cycl/grid in- terval) METHOD OF ANALYSIS 1.5-2.5:(wavenumber = 2.0/(N. DX)cycl/grid in- terval The gravity field values for a block of N x N (Pal. Khurana, and Unnikrishnan, 1979; equally spaced (gridded) data, are transformed from Hahn, Kind, and Mishra, 1976) the space domain to the frequency domain by means Elements of the matrix with 0.5 < of the 2-dimensional discrete Fast Fourier Transform (I 2 + J')~ : < = !.5 are averaged, and so forth to described by Cooley and Tukey (1965). the Nyquist wavenumber N/2 (see Fig. I). 549