THE MPE‘ECMEOK @F SPECTRRL EETEODS E‘s“? ”2'33 ANALYSIS m3 ER’TEBPEERTEGR 0? WWW DAM: 3: CRETECAL STEIN T313513 F83 "EEE DEGREE GE" PEER. fiECfifiGfifi SERIES {3&me WEEK? 3mm KEME 1973 This is to certify that the thesis entitled The Application of Spectra] Methods in The Analysis and Interpretation of Gravity Data: A Critical Study presented by Robert David Regan has been accepted towards fulfillment of the requirements for Ph D Geology degree in Date Feb. 1, 1973 0-7639 ABSTRACT THE APPLICATION OF SPECTRAL METHODS IN THE ANALYSIS AND INTERPRETATION OF GRAVITY DATA: A CRITICAL STUDY By Robert David Regan Gravity field data can be properly considered as probabalistic data or one outcome of a stochastic process. An essential assumption in the spectral analysis of such data is stationarity. Although gravity data are shown to be nonstationary the objectives of spectral analysis of these data obviate the requirement of stationarity. The transformation of involved gravity anomaly expressions such as those of the two-dimensional and three-dimensional prisms are easily calculated by means of a general transformation formula. The mathe- matical simplicity of the theoretical transforms of these and other ideal gravity anomalies (deterministic data) has distinct interpretational advantages over the more complex spatial domain equations. Despite these advantages, the finite length of practical gravity data precludes the direct application of the frequency domain equations. Calculation of theoretical transforms of segments of ideal anomalies is not feasible and a profile length of approximately six times the depth of the body is required for the finite transform to approximate the theoretical transform. An alternative spectral interpretation method, that is not as restricted by data length, is based on the convolution of the transform of the two-dimensional prism gravity anomaly with the rectangular and Bartlett spectral windows. .\‘ I 1’. 0 ‘ . v'r - I \ JT.’ I f! ‘9 o - ,- u . I i v I n , t . p . . I ., .0. . x I ' I III : l I I , ui— \- -}.. . a l. ‘ l‘ i i _. I . ,. l "I‘ ,. V! V I l I . v . 1 .‘ . ..~ \ - l ~'.l" ' 'I I O Q -. ..W ,. I. I . ‘—,.- ) -, -\‘1‘ . -J .- '4. .i I .l . . . r ' \‘—' .I .0.‘ - ‘. . ~ .7! I .‘II 4"- '. ‘- L "i; "h’lll ' mm! C ".,.“. .. so ' i .I1 . ‘ l l i- I '.l‘ I :.'.:. .-';’ I THE APPLICATION OF SPECTRAL METHODS IN THE ANALYSIS AND INTERPRETATION OF GRAVITY DATA: A CRITICAL STUDY By Robert David Regan A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Geology 1973 DEDICATION This dissertation is dedicated to my family, as small repayment for their many sacrifices. ii ACKNOWLEDGEMENTS I wish to thank Dr. R. J. Hinze for his invaluable direction, interest, and encouragement during the course of this study. Also the support of the U. S. Geological Survey Center of Astrogeology and in particular that of Dr. A. H. Chidester is gratefully acknowledged. TABLE OF CONTENTS LIST OF TABLES .............................................. LIST OF FIGURES ............................................. INTRODUCTION ................................................ Chapter I. SPECTRAL ANALYSIS .................................... Deterministic Data ................................ Introduction ................................... Applied spectral analysis ...................... Finite length ............................... Discrete function ........................... Digital spectral analysis ................... Probabalistic Data ................................ Introduction ................................... Time series, Stochastic process ............. Stochastic process functions ................ Stationarity ................................ Ergodicity .................................. Transformation of random data .................. Applied spectral analysis ...................... Fourier transform ........................... Niener-Khintchine theorem ................... Spectral estimates .......................... II. GRAVITY FIELD DATA .................................. Gravity Data as a Time Series .................... Probabalistic nature of gravity data .......... Probability density function ............... Stationarity of Gravity Field Data ............... Theoretical considerations .................... Examination of gravity profiles ............... Test of self-stationarity ..................... Rationale .................................. Test cases ................................. iv Page vii viii l 4 I7 Chapter 11. Profile results ................................. Model profiles ............................... Field data ................................... Effect of Nonstationarity on Gravity Analysis ...... Calculations .................................... Results ......................................... III. THEORETICAL GRAVITY ANOMALIES ......................... Transformation of Theoretical Gravity Anomalies .... Review of previous work ......................... General transformation .......................... Two-dimensional prism ........................... Theoretical transform ........................ Direct transformation ..................... General transformation .................... Method of interpretation ..................... Three-dimensional prism ......................... Theoretical transform ........................ Vertical line element ........................... Effect of Finite Data Length on Spectral Analysis Interpretation ..................................... Calculation of theoretical anomalies ............ Spectral windows ............................. Two-dimensional cylinder ..................... Rectangular window ........................ Higher order spectral windows ............. Two-dimensional prism ........................ General considerations ....................... General transformation formula ............ Direct transformation ..................... Length of data required to approximate theoretical transform ........................... Two-dimensional cylinder ..................... Slope ..................................... Intercept ................................. Profile length ............................ Two-dimensional prism ........................ Null positions ............................ Depth ..................................... Profile length ............................ Utility of rectangular window ................ Alternative Spectral Interpretation Method for Two- Dimensional Cylinder ............................... 45 93 Chapter Page IV. CONCLUSIONS ......................................... 99 Spectral Analysis of Gravity Field Data .......... 99 Spectral Analysis of Ideal Gravity Anomalies ..... 100 BIBLIOGRAPHY ................................................ 103 APPENDIX A: RUN AND TREND TEST ............................. 106 APPENDIX B: COMPUTER PROGRAM STEST ......................... 108 vi LIST OF TABLES Table Page 1. STEST - Test Case Results ............................... 29 2. Parameters of Model Profiles ............................ 35 3. STEST Resulta for Model Profiles ........................ 37 4. Gravity Attraction of Three Idealized Bodies: Spatial and Frequency Domain Equations .................. 48 vii I I . i ' . ' . - \ . . I . o l o ' . . ' I Figure mmwmm4>wm d—a—I N-‘O l3. 14. 15. 16. LIST OF FIGURES Two Test Series for Computer Program STEST ............ Model Gravity Profiles A and B ........................ Model Gravity Profiles C and D ........................ Gravity Field Profiles ................................ Direct and Indirect Spectral Estimates: Profile 2A ... Direct and Indirect Spectral Estimates: Profile 30 ... Direct and Indirect Spectral Estimates: Profile 48 ... Two-dimensional Prism ................................. An Example of the Method of Multiple Decay Spectra .... Vertical Line Element ................................. Four Spectral Windows ................................. Two-Dimensional Cylinder: Correlation Coefficient Between Calculated and Theoretical Transforms as a Function of k/D and Spectral Windows .................. Two-Dimensional Cylinder: Normalized Depth as a Function of k/D and Spectral Window ................... Two-Dimensional Cylinder: Tn(o)/G(o) as a Function of k/D ................................................ Two-Dimensional Cylinder: Percent Error in Width as a Function of 2k/b and Spectral Window .................. Two-Dimensional Prism: Correlation Coefficient Between Calculated and Theoretical Transforms as a Function of k/d and Spectral Window ................... viii Page 28 33 34 38 41 42 43 52 59 62 67 78 80 82 87 Figure 17. T8. 19. 20. 21. Two-Dimensional Prism: Correlation Coefficient Between Calculated (Rectangular Window) and Theoretical Transforms as a Function of k/d and d/D .......................................... Two-Dimensional Prism: Correlation Coefficient Between Calculated (Rectangular Window) and Theoretical Transforms as a Function of k/D .......... Two-Dimensional Prism: Percentage of Maximum Anomaly at x = 3b as a Function of M = d/b and N = D/b .............................................. Two-Dimensional Cylinder: T2(o)/T1(o) as a Function of k/D ...................................... Two-Dimensional Prism: Percent Error in Depth as a Function of k/D ...................................... ix Page 88 90 92 95 98 INTRODUCTION The use of spectral techniques in the analysis of gravity data has greatly increased in the last few years. Transformations of theoretical anomalies and raw field data have indicated distinct advantages in the use of the frequency domain. The mathematical simplicity of the frequency domain equations for theoretical anomalies provides interpretational advantages over the more complex spatial domain expressions. Determination of frequency content as an aid in the design of filters, and ease of continuation calculations are among the advantages gained in the transformation of gravity field data. However, the transformation techniques employed in these types of analyses are based on theoretical assumptions that impose limitations on practical data analysis. In general, two distinct limitations arise in practical spectral analysis. First, in the case of deterministic data (e.g., isolated gravity anomalies) the finite length of data necessitates the use of a data window to meet the Fourier transform requirement of infinite data length. Second, transformation of probabalistic data (e.g., gravity field data) is based on the Wiener-Khintchine theorem that requires the data be stationary, i.e., the statistical properties are time (space) invariant. As is the case with deterministic data, a data window (lag window) is usually employed in this transformation to satisfy the requirement of infinite data length. l 2 Some attempts have been made to obviate these requirements, such as the pseudo—power spectra suggested by Fraser, gt a1, (1966). However, the usual methods are, in the case of deterministic data, to assume infinite data length and, in the case of probabalistic data, to assume stationarity. The excellent papers by Odegard and Berg (1965) and Agarwal (1968) illustrate the use of these assumptions in gravity data analysis. The Fourier transformation of theoretical gravity anomalies (deterministic data) over idealized bodies by Odegard and Berg is based on an assumption of infinite data length. Agarwal considers gravity field data (probabalistic data) as stationary and transforms the autocorrelation functions to analyse the data in the frequency domain. However, the conclusions reached in these studies are dependent upon the validity (or practicality) of these assumptions. The main objective of this study is to examine the effect of these assumptions in the analysis of theoretical gravity anomalies and gravity field data. The principal assumption in the frequency analysis of gravity field data is stationarity. The highly correlated nature of gravity data would imply that these data are not stationary and the transforma- tion is suspect. To study the assumption, gravity field data are shown to be a time series, then tested for stationarity by a computer test of stationarity developed from several previously suggested methods. The effect-of nonstationarity on the power spectral calculations is studied by comparing direct and indirect spectral estimates. An important consideration in the transformation of theoretical anomalies is the requirement of infinite data length. The practical application of spectral techniques is dependent upon the transforma- tion obtained from data of finite length. The effect of finite data length is studied by considering the equations for the convolution of the transformed gravity anomalies with various spectral windows, and by determining the length of data vs data window required to approximate the ideal transform. Another concern in the application of spectral techniques to the interpretation of gravity data is the generality of the Fourier trans- form technique in the analysis of theoretical anomalies. The idealized bodies considered by Odegard and Berg (1965) had fairly uncomplicated mathematical structure and were easily transformed. To study the application of the transform technique, a general method for trans- formation of theoretical gravity anomalies, modified from Battacharya (1966, 1967), is presented and transformations of gravity anomalies over several idealized bodies, not previously studied, are calculated and methods of interpretation determined. CHAPTER I SPECTRAL ANALYSIS The analysis of data in the frequency domain, i.e., where frequency is the independent parameter, is termed spectral analysis. Although several frequency functions, such as the phase spectrum and power spectrum, are commonly employed in spectral methods, the funda- mental calculation is the initial transformation of the data into the ‘frequency domain. The limitations imposed by this transformation persist in any further spectral calculations. The nature of the data determines the method of initial trans- formation and subsequent frequency analysis techniques. In determin- istic data the transformation of the original data is calculated whereas with probabalistic data the frequency domain is usually accessed through the covariance functions. Limitations on practical spectral analysis are imposed by the theoretical assumptions necessary for the calculation of these transforms. In both cases the Fourier transform requires continuous data of infinite length and the probabalistic data are further constrained to be stationary. Because of the dual nature of the gravity data considered in this study, i.e., the deterministic nature of theoretical anomalies and the probabalistic nature of field data, it will be necessary to briefly discuss the basis, and subsequent limitations in the applied case, for the transformation to the frequency domain of both types of data. Deterministic Data Introduction Deterministic data is defined as that which is capable of being represented by a strict mathematical expression. Because of this mathe- matical expression, the future values of these data are exactly predictable, i.e., the data can be extrapolated and interpolated with zero error. An example of deterministic data would be the gravity attraction of a perfectly smooth, non-rotating, uniformly dense sphere. The gravitational attraction at any point in the space surrounding the sphere can be precisely calculated, and a measurement would only confirm the previously calculated value. The transformation of deterministic data to the frequency domain is accomplished by the Fourier transform. This transform assumes con- tinuous data of infinite length and the transform exists if the Dirichlet conditions are satisfied, i.e., the function is piecewise continuous and absolutely integrable. These conditions are usually satisfied if the data are physically realizable, and are satisfied in the case of individual theoretical gravity anomalies. Because the value of the transform is generally complex, the results are usually represented as some real function of the transform. Such functions as the amplitude or power spectrum and the phase spectrum are commonly used. The relevant definitions and formulae are: g(x): continuous function in the space domain, extending from «n to +m and square integrable G(w): Fourier transform of g(x),= %%_ f g(x)e'xw dx G(w) = R(w) + iI(w) 6 where w is the angular frequency (an), R(m) is the real part of the transform, I(w) is the imaginary part of the transform, Amplitude Spectrum = (R(w)2 + 1(e)2)f/2 Power Spectrum = G(w) . G(w)* = |R(w)2 + I(w)2|, Phase Spectrum = tan']-%gfi( , and * represents the complex conjugate. Some simplicity or ease of calculation is realized when g(x) is a real even function for then the imaginary part of the transform is zero. Applied Spectral Analysis In the case of most theoretical calculations the transforms exist and are readily determined through the Fourier integral and its many properties. However, the discrete finite characteristics of applied data analysis impose limitations on the transformation that must be considered. These limitations are best detailed by considering each one separately as it affects integral transforms and then detailing the more relevant aspects of digital spectral analysis. Finite length_ If only a finite segment, gf(x) = g(x) [-k,k], of the continuous function g(x) is available for transformation, the requirement of infinite length is usually met by multiplying the segment by a data window. Basically a data window can be defined as: o x<-k h(x) = f(x) -k5x5k o x>k 7 with a maximum value of l at x = o, and k w ng(x)f(x)dx= (g(x)h(x)dx. Multiplication of two functions in the spatial domain is equivalent to the convolution of their Fourier transforms in the frequency domain. The transform T(w) of the finite segment of g(x) is so related to the theoretical transform G(w). T(w)=%—1-T- ( gf(x)f(x)e'wxdx =;—— ( g(x)h(x)e'mxdx -k .00 T(w) = G(w) * H(w). As the length of the segment tends to infinity, gf(x) approaches g(x) and T(w) approaches G(w). At any particular segment length T(w) is considered an approximation of G(w) and the effect of the data window can be significant. Discrete function_ The effect of discrete data on the Fourier transform calculation can be illustrated by considering a sampled version of the theoretical function g(x). The sampling function can be defined as a Dirac comb 00 1(X) = :E:: 6(x-nA) n=-w 8 wichI a fixed sample interval, and the sampled theoretical function is 91(X) = 9M - i(X) with transform G.(w) 1 G(w) * Hm) (I) ; G(w-T) ;— Z 6(T- D—i-F-N =-oo ":4” which indicates that discrete sampling of an infinite function produces a band limited spectrum, i.e., gi(x) has a transform that is uniquely determinable only for lwl 5 %—(f 5 E3) (Nyquist frequency). The digital sampling of a space-limited function, i.e., a signal of finite length 2k, further constrains the calculations of the trans- fbrm to frequency intervals of Am = E-(Middleton, 1960, p. 210). Digital spectral analysis The amount of calculation involved in the spectral analysis of physical data usually requires the use of digital computer techniques. In this case the digital Fourier transform, or the discrete represen- tation of the Fourier integral, is evaluated. The application of these methods requires that the signal is digital and of finite duration. The signal can be considered as: gn(x), n = o, l, ......... N, i.e., N + 1 samples with constant sample interval Ax: length of signal is NAx. Then the digital transform is fiznmn] Gm(‘*’)= {if}: 9,, (X)e N“ where m = -N/2, ..... -l, 0, l, ......... N/Z. The calculation is limited by the discrete, finite nature of the signal in that the spectrum is determined only at multiples of 2n/NAX for le 5 n/Ax and represents the digital representation of G(w) * “(w)- Further limitations imposed by numerical analysis arise in the evaluation of the Fourier transform. The direct evaluation of the digital transform, equivalent to the rectangular rule or Simpson's method of numerical integration, produces inaccurate spectral estimates. Hanning [1962] shows that this result occurs in numerical evaluation of integrals of the type b b f f(x) coswx dx or { f(x) simx dx, a a and outlines Filon's method to be used in this case. Another, and more popular, method of evaluating the digital Fourier Transform is the Fast Fourier Transform suggested by Cooley and Tukey (1965). The original Fast Fourier Transform is restricted to N = 2n samples, routinely applied by adding zeros to the data set, and provides the accuracy of Filon's method with an increase in speed of computation. Recently the Fast Fourier Transform has come to mean a class of Fast Fourier Transforms, some operating on any number of data points [e.g., Cutler, 1970]. The main problem with a Fast Fourier IO Transform acting on any number of data points is that while maintaining speed of computation, in general, the accuracy can suffer. Indeed these transforms degenerate to the rectangular rule when the number of samples is prime. Probabalistic Data Introduction In marked contrast to deterministic data, probabalistic (or random) data are not representable by a strict mathematical expression and are not exactly predictable. Every observation is uniformly singular and is the result of a probability law, that is any observation will represent only one of many possible results. A series of observa- tions is called a sample function (or time series) and the collection of all possible sample functions is termed a random (or stochastic) process. A geophysical example of a random process would be the geomagnetic field (external and internal) with individual magnetograms representing sample functions of the process. The methods and objectives of analysis of probabalistic data differ from those employed with deterministic data. Random data analysis does not consider the sample function by itself but instead considers it as one of the set of all possible sample functions, and the goal is to determine the underlying process. However, it is realized that the process can never be sampled with absolute certainty and the functions calculated from the observations are considered statistical estimates of the associated process functions. The various theories of time series analysis are basically concerned with improving these estimates, or providing bounds on the estimates so that the true process function may be determined. 11 Before spectral analysis of random data can be considered it is essential to discuss some of the concepts and terms pertinent to random data. Time series, Stochastic process A strict definition of a time series, or random function, is that it is one realization (outcome) of a stochastic (random) process. A more common definition is that it is a set of probabalistic observa- tions arranged sequentially. The label "time series“ is perhaps a misnomer, and is directly applicable only when the observations are made chronologically. However, for the independent variable, time, there may be substituted any other parameter, e.g., distance. The basic idea of the statistical theory of time series analysis, or random data analysis, is to regard the time series as a set of observations made on a family of random variables, i.e., for each t in T, x (t) is an observed value of a random variable. The set of observations {x (t); t c T} is called a time series. A time series is more properly identified as {x(t,s); teT, 535} where S is a probability space. Thus, the collection of all random functions, the random (or stochastic) process can be considered as being composed of a family of random variables, and viewed as a function of two variables t and s. The sample space associated with this stochastic process is doubly infinite and the set of time functions which can be defined on this space is called an ensemble. 12 Stochastic_process functions The basic functions that describe the Stochastic process are: a) mean of the_process ux(t)g E{X(t)} = L ;’ x(t)f(x,t) dx where L j! is the Lebesgue Integral, f(x,t) is the Probability density function, and E is the Expectation operator. In the general case the means are different at different times and must be calculated for each t. b) autocorrelation function of the prggess Rx(t].t2)g E{x(t1). x(t2)*} where x(t2)* is the complex conjugate of x(tz). c) autocovariance function g:_the_process * Cx(t],t2)c=l E{[x(t]) - ux(t])], [x(tz) - ux(t2)] } In the general case these functions must be calculated for each t1, t2 combination. Stationarity In general the properties of a stochastic process will be time dependent. A simplifying assumption which is often made is that the process has reached some form of steady state in the sense that the 13 statistical properties of the process are independent of absolute time. Stationarity can be pictured as the absence of any time-varying change in the ensemble of member functions as a whole. The stationarity of an individual time series, rather than the entire process, is termed self-stationarity. A sufficient degree of stationarity for most time series analysis is wide-sense stationarity. A process has wide-sense stationarity if its mean function is a constant and its autocorrelation depends only on the time difference, t1 - t2 = T, and not on the absolute value of the respective times, i.e., E{x(t)} = ux(t) = constant and E{x(t+T), x(t)*} = Rx(r). Wide-sense stationarity is also termed stationarity of the second order, i.e., the series is stationary through its second order statistical moments. Stationarity of order n implies that all statis- tical moments less than n depend only on the time differences. Ergodicity Ergodicity relates to the problem of determining the statistics of the stochastic process from the statistics of one time series. A stochastic process for which the statistics are thus determinable is said to be ergodic, and the single time series is representative of the ensemble. The ensemble moments can then be equated to the time moments. 14 For example the time average equals the ensemble average of an ergodic process k 14x”) = k { mm = E{x(t)} Nl—I If a time series is ergodic we need only to measure the time averages which are available rather than the postulated ensemble averages. Because ergodicity is a subclass of stationarity, a time series must be shown to be stationary before the question of ergodicity can be considered. Transformation of Random Data In the use of spectral analysis with random data the Fourier transform of the autocovariance function is routinely calculated rather than that of the observed data. The reason for this approach is that the principal objective of random data analysis is to obtain information on the process and by transformation of the autocovariance function the process power spectrum can be obtained. While Fourier transformation of the observed data is possible, what is obtained is an equivalent sample function in the frequency domain. Although the power spectrum of this sample function can be calculated, i.e., x(w)x(m)*, it is not possible to obtain the process power spectrum, nor any other process functions, in this manner [Middleton, 1960; Thomas, 1969]. In the case of stationary (wide-sense), ergodic random processes, the Wiener-Khintchine Theorem [Middleton, 1960] relates the auto- covariance function of the process to the process power spectrum through the Fourier transform. This is the basis for spectral analysis of random data. 15 Applied Spectral Analysis The limitations in the case of random data can be detailed by considering the two theorems used in spectral analysis and then considering the applied case of spectral estimation. Fourier transform The limitations on the Fourier transform discussed in the section on deterministic data apply here. The theoretical assumptions of the Fourier transform are independent of the nature of the data. Wiener-Khintchine theorem The assumptions involved in the Wiener-Khintchine theorem are statistical process functions and unlike those of the Fourier transform the error introduced by approximations to these assumptions cannot be fully detailed theoretically. If the process is not wide-sense stationary the process power spectrum cannot be calculated. In this case it is not the introduction of error into the calculations that must be of concern but the total inapplicability of the Wiener- Khintchine theorem. The fact that the process functions can never be determined with probability one and that there is no universally applied test for the property of stationarity has led to the employment of several altern- ative approaches in the routine spectral analysis of random data. 1) The process is assumed stationary, 2) Any nonstationary aspects are assumed to be removed by detrending the data, 3) The samples are segmented and each section considered as piecewise stationary. 16 In the event of known nonstationary data the process power spectrum cannot be calculated and alternatively a generalized or instantaneous spectrum is calculated [Bendat and Piersol, 1966, chapter 9]. Spectral estimates The field of spectral estimation, or mechanics of spectral analysis, is properly considered as the statistical equivalent of digital spectral analysis of deterministic data. The methods of spectral estimation are primarily concerned with the estimation of the process power spectrum from the estimate of the autocovariance function of a discrete, finite sample from a stationary, ergodic random process. In spectral estimation the limitations of practical data analysis with the Fourier transform are coupled with limitations of determining the autocovariance function from a finite length sample. Studies of the effect of lag windows (data windows) on the spectrum are the basis of several excellent texts, e.g., Jenkins and Watts (1968), and the results are equally applicable to deterministic data. The relation between statistical estimates in each domain is detailed by confidence limits in spectral calculations. Chapter II GRAVITY FIELD DATA The use of the frequency domain has been shown to haVe distinct advantages in the filtering, continuation, and interpretation of gravity data (e.g., Dean, 1958; Byerly, 1965; Mesko, 1965; and, Kanasewich and Agarwal, 1970). In the application of these types of analyses, field data are considered to be a sample of a stationary random process (e.g., in gravity data; Goldstein, 1962; Agarwal, 1968; and analagously for magnetic data; Horton gt_al,, 1966; Spector and Grant, 1969) and the covariance function transformed. The central assumption in this approach is that of stationarity. The purpose of this chapter is to detail the validity and effect of this assumption. Specifically, the probabalistic nature of gravity field data is presented and its stationarity investigated from theor- etical and empirical considerations. Thus, the effect of nonstationarity in the transform to the frequency domain is studied. Most studies on spectral analysis of probabalistic data are concerned with one-dimen- sional (vector) discrete time series. To be consistent with this work, and fer the sake of mathematical simplicity, one-dimensional gravity field data (i.e., profiles) are considered in this study. Robinson (1967) indicates that the results would be equally applicable to two- dimensional data. 17 [I!.I";1I,| If .~’ 1',”I ’I‘r If"1 11.1 'I‘ I‘ In 1111 " I"! i I (‘4'! .1211. "~‘I 'i I. ‘ 1'» ‘1'! I 1'. 1:71.11 {'11.11’H 1‘ 1. 1', ,[u '1 1 .I'.I‘H~ '5' 151“" '11:! It): 1 '1' 'L'I I1. . Iv.'tfzr~‘ :31 .'1;(1‘1',’n111.'11 .‘»,I’E‘1‘HH:'I 1"111“ I'm. ”Heliuv -1.’T1 \fl’i‘lvum In ”11111.11 .1: I gun‘j'l 'IH’ 11”,”, ,' \,.'uiv1.’~m ",1![\'1i I, lv~swi.‘1wnilw~ ,[I'I [1.1I',(!H'((I :in) I "13’" II 1: '1 .w 1‘ 1 I 1. . "I «I ‘ .. I ll {1‘ 3 r :v'vul/Ir. 1 ’ ‘I ..1 . ,1 I. "r I [DH-‘1’; 1 H K *‘I .I:.‘..”l .' “w, Ellis." . '1 ' *r '1 1.~ 2': t‘w‘. , '1‘ , '1'.1. ‘1': ' [If/II '1 ’1 1) HYI" f(‘Ili'I/I 1;? 111.! 113-.1111. \ml ‘_'I.'Hm\Y-.I1.'I. 31:1 '1" ' 1:1: Haw; .» 1 ‘2 INT. (' ('1 , 1.1.5.15 Tun-1111. iii} *1 11‘1I1'1"I..'«'11 -1 WWI-1.1M 91:51 1:1 ".":"“1|"! mi? 1 .J in .‘ ‘i , ._ ‘ ’ . l): i(:"'I‘¢ ’ I1 1 1:13. . I Inif: 1.'211:'...7<3‘Y(; .‘e'r M131) II-n'r] -'il‘*1.'1."""ir‘i!ti‘ 11,.i'1i'3'w 211:) 11.111 '1“:qu 0111711 01) iii I13“..fl,)'t"3{ 11') 18 Gravity Data as a Time Series To consider a gravity profile as probabalistic data is to equate it to one realization of a stochastic process, or a time series. Toward this end, the probabalistic nature of the data is argued and a probability density function postulated. Probabalistic Nature of Gravity Data The accurate classification of data as either deterministic or probabalistic is a major problem. It may be argued that purely random or purely deterministic data do not exist, and that all data having some elements of randomness or determinism can be placed anywhere in the spectrum between these two ideals. Actual classification is usually determined by considering the attributes of the available data and comparing these to the qualities of deterministic and random data. Robinson (1967) states that observed data can be classified as deterministic provided it has no features of randomness, and classified as random provided it has any features of randomness. The principal intention here is to show that a probabalistic interpretation is not inconsistent with the structure of gravity field data. If gravity field data, or the results of an exploration type survey, are perfectly obtained and reduCed, i.e., no errors in readings, location, or reduction, then the results are representative of sub- surface mass distribution. Certainly, the location of geologic bodies of varying mass is not deterministic, i.e., there is not a general mathematical expression relating their locations, size, and shapes, 19 and from observations at one location these parameters are not predictable. Ideal gravity field data then must also be considered nondeterministic. The relative sparsity of gravity stations and inaccuracies that can arise in obtaining the reduced values can only add to the character- ization of actual gravity field data as probabalistic. The fact that, a priori, the value that will be obtained from a gravity field reading is not known and that the data cannot be extrapolated with zero error is further indication that a probabalistic approach would be applicable in gravity data analysis. With gravity field data the underlying process would be the world wide, or geologically possible, distribution of subsurface bodies of varying mass and dimensions. A suite of measurements would be the representation of a sample drawn from this population according to some probability law. Probability density function Although the exact nature of the probability law is not critical to spectral calculations an assumption of a probability density function, a mathematical representation of the probability law, should be made. A probability density function describes the probability that a random variable can assume certain values. Spector (1968) assigned a uniform probability model to the parameters of geological bodies causing aeromagnetic anomalies. Naidu (1967) in considering the potential field and statistical properties of mass density distribution over random media assumed a normal (Gaussian) distribution for the density distribution. Because each of ._ ——-- .__ .. .4...— these cann< COT“ dens assu of t stuc imp' gra' obs The Sow 8C1 ls Ce DY 20 these is considered likely, and an exact probability density function cannot be known, it is a moot point as to which is to be assumed correct. There are, however, two reasons for utilizing a normal probability density function. The basis of correlation theory rests on the assumption of a normal density function, and the lower order moments of the normal distribution completely describe the process. In this study a normal (or Gaussian) probability density distribution is implicitly assumed for the parameter of geologic bodies causing gravity anomalies. Stationarity of Gravity Field Data In time series analysis of gravity field data, or any body of observational data, the fundamental assumption is that of stationarity. The stationarity of gravity field data is investigated by considering some of the properties of the data and by testing some model and actual field profiles. Theoretical Considerations The primary reason for considering gravity data to be nonstationary is that the underlying probability function is space dependent. In certain geographical areas the gross characteristics of a gravity profile will be known a priori because of previous knowledge of the geology and the gravitational response to that geology. This is not to say that the profile will be known exactly or even to a good approximation, but, a certain range of anomalies can be anticipated. 21 Other considerations also suggest that gravity data can be considered nonstationary. This result is based on the observation that gravity data possess two characteristics of nonstationary data. These are highly correlated data and significant power in the low frequency end of the spectrum (Jenkins and Watts, 1968). The highly correlated nature of gravity data cannot be demon- strated mathematically by the calculation of a general autocorrelation function. It is rather an intuitive consideration based on the potential nature of the data, where it is known that successive observations are not independent of each other and are observed to be highly correlated. Indeed this fact has led some workers (Kaula 1969) to consider Markov theory in the analysis of gravity data. Substantial power in the low frequency range of the spectrum Dw. where D is some constant is indicated by a term of the form e' related to the depth of the body, observed in the transform of gravity anomalies (Odegard and Berg, 1965; this study, chapter III, and in the analysis of magnetic anomalies by Battacharya, 1966). This term assures power in the low frequency end of the spectrum and is in agreement with the suspected high correlation. This can be seen through the Fourier transform relationships of convolution in one domain being equivalent to multiplication in the other. Thus, the -Dw component of the autocorrelation function related to the e term -Dwe-Dw. would be the transform of e This transform would be 2 2 2 4D|(4D +4n r ) (Erdelyi, 1954). 22 Two objections may arise from the assumptions made in equating a . gravity profile to a time series. Because a normal probability density distribution has been assumed does this automatically make the process stationary? If only a time (space) independent probability density function had been considered, then the answer would be in the affirm- ative. However, the choice of the normal probability distribution does not imply a time (space) independence. It can be time (space) dependent, that is it can have a variable mean value function and it can have a space dependent variance. In reality this would correspond to different parameters associated with various geologic provinces which appears reasonable. Another question is whether the capability to reoccupy a point and repeat a reading implies stationarity of the profile. Because distance has been made the independent parameter (equivalent to time in the usual theory) the answer is no. Going back in space would be equivalent to going back in time. Reoccupying a point in the gravity profile would be equivalent to returning in time to the time a random number was generated in a time series. In both instances going back in either the time or space domain results in repeatability. Examination of gravity profiles The theoretical concept of stationarity as outlined in the first chapter would indicate that it is almost impossible to determine stationarity because of limited knowledge of the process. However, in the transformation from theoretical to empirical concepts, theory must be tempered with reasonable practical considerations. In most practical applications the single observed time series is the only information 23 available on the parent process. Hence, ergodicity must be hypothesized and time domain statistics utilized. Also all analyses are performed on the sample series and thus the stationarity of this series is of prime interest. Bendat and Piersol (1966) have termed the concept of stationarity of this one time series as self-stationarity. Thus we speak of the series being stationary rather than the process being stationary. However, the concept of self-stationarity is not restrictive. A necessary condition to extend this self—stationarity to stationarity is that the process be ergodic, i.e., that the series is representative of the process. Ergodicity is impossible to prove (except in special instances) when the entire process is not known. Bendat and Piersol (1966, p. 12) state that “in actual practice, random data representing stationary physical phenomena are generally ergodic". If the assumption of ergodicity is justified then self-stationarity becomes equivalent to stationarity, since the single time series is representative of the ensemble. The concept of self-nonstationarity is no less restrictive, although ergodicity cannot apply in this case, most observed non- stationary processes are also self-nonstationary (Bendat and Piersol, 1966, p. 13). Test of Self-Stationarity To date no rapid, accurate, routine method has been available to test the stationarity of a time series. Such a test is clearly required so that the stationarity or nonstationarity of a series can be readily determined and the validity of the application of time series 24 analysis be established. The method outlined in this section is designed to fulfill this requirement. The test for self-stationarity is a computer programmed test, based on the methods proposed by Bendat and Piersol (1966) and Bryan (1967). The basis of these methods is that in a stationary series certain statistical pr0perties of the time series are considered invariant with time. The tests are for second order (wide-sense) self-stationarity. Bendat and Piersol (1966) suggest that the series be divided into M equal time intervals (either contiguous or non-contiguous) and that the mean and variance of these intervals be calculated. The two series thus formed, composed of means and variances, are then tested for underlying trends or variations by the Run test and the Trend test. (Appendix A). If no trends or variations are suggested by the applica- tion of these tests, the original series is assumed to be stationary. Basically, Bryan's test (1967) is quite similar in that he tests the invariance of the means and variances obtained from independ- ent, equal length segments. Rather than using the sample mean and sample variance he constructs two combinations of the data to serve as estimates of the population mean and population variance. Those two estimates are independent, and independently distributed. Using these two variates, m-a linear function of the data and an unbiased estimate of the population mean, and Q-a quadradic function of the data, he develops a test for the hypothesis that the time series is stationary, and two test variables, L1 for the Neyman-Pearson L test, and F for the F-distribution test. 25 Thus, with a combination of the two tests there exist four independent test statistics for stationarity. A rigorous test of stationarity is obtained by determining an acceptance criterion to be that all test statistics must indicate stationarity. In the proposed test both methods are combined and used with some restrictions. The tables for the Run and Trend tests (Bendat and Piersol, 1966, p. 170) were extrapolated to include the range M = l to M = 200. With values of M less than 12 the results proved unreliable. Hence, a low limit cut off at M = 12 is utilized. The acceptance region was extended to include the lower bound. The Bryan test was extended to include the 97.5 percent confidence interval. Rationale The logic followed in the test of stationarity in computer program STEST as given in Appendix B is: First the sampling interval for independent samples is determined. If independent samples cannot be determined, i.e., the autocorrelation function does not go to zero, the test is aborted and the series can be considered nonstationary if a reasonable number of points has been used. Once the sampling interval has been determined, the series is segmented into independent samples of length N. Initially the series is tested with N = 5, then N is increased to 10, and the final test, if there are enough data points, is for N = 15. The minimum test is for M = N samples. If there are not enough data points for this number of samples, the requirement for independent 26 samples is relaxed slightly (i.e., the sample separation interval is steadily decreased to a limiting value equal to the number of lags necessary for the autocorrelation function to go to 0.1). If at this sample interval there are not N samples, the test is aborted. If this happens at N = 5, it may be an indication of non-stationarity. If we have M independent samples of length N, the series is tested for stationarity at three confidence intervals (95%, 97.5%, 99%) in the following manner: a) if M is greater than or equal to 12, the Bendat and Piersol test statistics and the Bryan test statistics are both utilized. b) if M is less than 12 only the Bryan test statistics are utilized. The series is considered stationary if both the Run and Trend tests show no trends or variations for the mean and variance series and if the two test statistics in the Bryan test indicate stationarity. In the case where only one method is utilized (i.e., M < 12) station- arity is tested on the merits of the Bryan test statistics alone. It should be noted that the 95 percent confidence interval is the most restrictive (i.e., the smallest acceptance region) and the other confidence intervals progressively less restrictive. Also results at the largest N used are preferable since more data points are used in each sample to determine the test statistic and assumption of normality in the Bryan test statistics is more closely approximated. 27 Test cases Seven series with known stationarity characteristics were used as test cases. Box and Jenkins (1970, Appendix) detail six series (labeled A through F) that are known to be stationary or nonstationary and used as examples throughout their book. The seventh test series is the sample of a second order autoregressive process as given in Jenkins and Watts (1968). One of the nonstationary samples of Box and Jenkins and the stationary sample of Jenkins and Watts are shown in Figure l. The results of the program STEST are shown in Table 1. In all cases, except series A, the test accurately indicated the stationarity or nonstationarity of the series. The discrepancy in series A may be attributable to the fact that correlated (i.e., not independent) samples were used. It is interesting to note that the series generated from the second order autoregressive process is stationary and the process itself is stationary. Thus, in this case self-stationarity is indicative of stationarity. Profile Results Four model and three actual field profiles were used to test the hypothesis that gravity data is nonstationary. Model profiles The model profiles were constructed using spheres and two- dimensional cylinders of arbitrary depth, radius, and density contrast placed at various locations in a 30 kilofeet profile. Two profiles had regional structures and two did not, and gravity stations were 28 Emmhm Emgmoca sausanu gem mowemm pmoh 03H "F beamed mmcvummm ogzgmgmaeme mmmooga Fmowsmgo Amcwxcmc use xomv u mmwgwm Am cow omp 00— cm 0 d q q 4 _ _ _ H A j A A q _ q q _ H 4 q a . . % o 0.... :00... 0 cos. c 1 FN + Nrpx m.o . FAx u “x mmmuoea m< gouge ucoomm we mmzpm> omF A< cm 7 O_m_- « 1 d a A a a _ _ O_O_- « a _ A _ _ _ q _ fl _ H . q a a q _ — O O 0 lg 0. o o o o o o o o I O 0 O. 0 o O 000 o o O o 00 O 00 000000 00 00 O o. L “NF-'0!— w_ om um um mm ml 29 TABLE 1 STEST - TEST CASE RESULTS Series A - Nonstationary (Box and Jenkins, 1970) STEST Results - Correlated samples used - A) N = 5 95% confidence interval : nonstationary 97.5% confidence interval : stationary 99% confidence interval : stationary B N = 10 Not enough data points for independent or correlated samples. Series B - Nonstationary (Box and Jenkins, 1970) STEST Results There are not enough data points for independent or correlated samples. Since this occurred for a sample of length 5 and the length of the input series is 369, this may be indicative of nonstationarity. Series C - Nonstationary (Box and Jenkins, 1970) STEST Results A) N = 5 95% confidence interval : nonstationary 97.5% confidence interval : nonstationary 99% confidence interval : nonstationary B N = 10 Not enough data points for independent or correlated samples. 30 TABLE 1 (cont'd) Series D - Nonstationary (Box and Jenkins, 1919) STESTyResult§_ A) N = 5 For N = 5 there are not enough data points for independent or correlated samples. Since this occurred for a sample of length 5 and the length of the input series is 310 this may be indicative of nonstationarity. Series E - Stationary_(Box and Jenkins, 1970) STEST_Results A) N = 5 95% confidence interval : nonstationary 97.5% confidence interval : stationary 99% confidence interval : stationary B N = 10 Not enough data points for independent or correlated samples. Series F - Stationary (Box and Jenkins, 1970) STEST Results A) N = 5 95% confidence interval : stationary 97.5% confidence interval : stationary 99% confidence interval : stationary B) N = 10 Not enough data points for independent or correlated samples. _§gcond Order Autoregressive Process - Stationary_(Jenkins and Watts, 1968) STEST Results A) N = 5 95% confidence interval : stationary 97.5% confidence interval : stationary 99% confidence interval : stationary 31 TABLE 1 (cont'd) B) N = 10 95% confidence interval 97.5% confidence interval 99% confidence interval C) N = 15 95% confidence interval 97.5% confidence interval 99% confidence interval stationary stationary stationary stationary stationary stationary 32 assumed at 100 foot intervals. The gravity profiles are shown in Figures 2 and 3, and the parameters used are detailed in Table 2. All four profiles tested nonstationary and the computer output detailing this is presented in Table 3. Field data Three field profiles were selected from available gravity data obtained at Meteor Center, Arizona; Zuni Salt Lake, New Mexico; and S. P. Crater Quadrangle, Arizona (Figure 4). The different character of the subsurface structure of these areas, reflected in the gravity profiles, offers a good field sample. In all three cases the reduced field values were used. All three profiles tested nonstationary before and after removal of any trends. Effect of Nonstationarity on Gravity Analysis The effect of nonstationarity in the analysis of gravity data was studied by comparing the power spectra calculated by both the direct and indirect methods. The direct method involves the calculation of the product of the Fourier transform of the data with its complex conjugate (usually employed with deterministic data). The indirect method is the standard method for calculating the power spectrum of random data, i.e., the Fourier transformation of the autocovariance function. For finite length samples of stationary random processes, the mathematical equivalence between the two methods can be demonstrated (Robinson, 1967; Jenkins and Watts, 1968; Box and Jenkins, 1970). 33 m use < mmpmeoca apw>mcw quoz "N weaned 8 cm _ . BEN . .89. s Doom: 1 m mrwwoca 88m peed < m—vmocm comer o~.o mN.o o¢ om on on om slebw slebw 34 a was 0 mmewoca xpw>mco quoz "m weaned oooom oooom . seed oooo_ a _ _ _ _ 4 _ _ — fl _ _ _ _ 1 o mrwwogm 1 oooom oooom I u wremogm or m_ mm mm slefiw slefiw 35 TABLE 2 PARAMETERS OF MODEL PROFILES All units in Kilofeet--density contrast 0.5 gm/c.c. Profile A (Figure 2) a) Two-dimensional cylinder centered at x = 5 depth = 9.25 radius = 9.195 b) Sphere centered at x = 9 depth = 2.0 radius = 1.05 c) Two-dimensional cylinder centered at x = 10 depth = 0.85 radius = 0.662 d) Sphere centered at x = 15 depth = 10.0 radius = 8.6 e) Regional trend X(I) = I/10 I = 1,300 Profile 8 (Figure 2) a) Two-dimensional cylinder centered at x = 5 depth = 0.25 radius = 0.1 b) Two-dimensional cylinder centered at x = 9.4 depth = 0.15 radius = 0.06 Profile C (Figure 3) a) Two-dimensional cylinder centered at x = 5 depth = 0.25 radius = 0.1 36 TABLE 2 (cont'd) b) Sphere centered at x = 9 depth = 0.25 radius = 1.05 c) Two-dimensional cylinder centered at x 10 depth = 0.85 radius = 0.662 d) Sphere centered at x = 15 depth = 1.0 radius = 0.6 e) Two-dimensional cylinder centered at x 18.3 depth = 0.25 radius = 0.195 f) Sphere centered at x = 26.4 depth = 1.0 radius = 0.6 Profile 0 (Figure 3) a) Two-dimensional cylinder centered at x 7.5 depth = 0.05 radius = 0.041 b) Two-dimensional cylinder centered at x 10.5 depth = 0.75 radius = 0.5 c) Sphere centered at x 15 depth = 10 radius = 8 depth = 2 l 0 6 d) Sphere centered at x = 25 . 2 radius = 1.05 37 TABLE 3 STEST RESULTS FOR MODEL PROFILES Profile A (Figure 2) There are not enough points for independent or correlated samples. Because this occurred for a sample length of 5, and the length of the input series is 300, this may be indicative of nonstationarity. B) N = 10 C N = 15 Profile B (Figure2) A N = 5 95% confidence interval : nonstationary 97.5% confidence interval : nonstationary 99% confidence interval : nonstationary B) N =_y1 Profile C (Figure 3) A) N Not enough points for independent or correlated samples 5 95% confidence interval : nonstationary 97.5% confidence interval : nonstationary 99% confidence interval : nonstationary 10 Profile 0 (Figure 3) Not enough points for independent or correlated samples There are not enough points for independent or correlated samples. Because this occurred for a sample length of 5, and the length of the input series is 300, this may be indicative of nonstationarity. mgals 38 -185- -186- A) Zuni Salt Lake -l87—- -188_- -189 Q ..2900 Feet -155 — — -160- B) SP Crater A -165 - m ~170- '3 U) E -175 — -180 F Feet -185 LIIlLllllllllllLllllllLll -170-— -171'- C) Meteor Crater I -l72-— ~ U) E» -173" ‘ 9 I . .2900 J -174-— Feet _'I75lllLLllllllilIlllJllJllllLAIIIIIIJJI Figure 4: Gravity Field Profiles 39 This approach was selected because the direct method permits the calculation of the power spectrum without the assumption of stationarity. And, it was anticipated that the two spectral estimates of gravity data obtained from these methods would differ due to non- stationarity. The effect of nonstationarity could then be determined by a comparison of the two results, and its effect on other frequency calculations detailed. Calculations Two model profiles (Figure 2A, Figure 30) and one field profile (Profile B, Figure 4) were selected, assumed to be stationary and their spectral estimates calculated by both methods. In a calculation of this type with digital data of finite length, smoothed spectral estimates must be calculated. This involves inter— polation of the direct spectrum and use of both data and spectral windows (Brumbach, 1968). The calculations were done in a computer program using the rectangular rule form of the digital Fourier transform. This form was chosen to assure mathematical equivalence, particularly with the interpolation of the direct spectrum, and to obtain spectral estimates at the same frequency interval in both methods. Several combinations of data and spectral windows (as defined in Chapter III) were utilized to study any possible effect on the results. Although the direct method is not usually used in the analysis of pure random data because of the inconsistency of the estimates (Jenkins and Watts, 1968) this problem was not observed with the gravity data. Consistent spectral estimates (i.e., a decrease in the variance of the spectral estimate as the number of data points increases) were obtained even with as few as 50 field points. 40 Results In all test cases the spectra calculated by both the direct and indirect methods are virtually identical. There was no effect of the various combinations of windows on the degree of similarity between the spectral estimates. The results are shown in Figures 5, 6, and 7. The difference between the two methods at the higher frequencies is not considered significant and is perhaps related to the use of the rectangular rule in the numerical integration. An explanation can be offered, pgsteriori, that indicates why these spectra do agree. In the analysis of a single profile (single time series), under the assumption of stationarity, the autocorrelation function calculated is the time averaged autocorrelation function which is a function of T only, whether or not the process is stationary (Middleton, 1960, p. 49-50, p. 125). The Fourier transform of this function is mathematically equivalent to the direct method calculations. In terms of time series analysis this means a sample power spectrum (i.e., the power spectrum of one realization of a stochastic process) has been calculated. In spectral analysis of random data the sample power spectra calculated by the indirect method converge to the process power spectrum only if the process is wide-sense stationary and ergodic (Wiener- Khintchine theorem). Because the aim of time series analysis is usually to investigate the underlying process, the process power spectrum is of prine:importance. But, the aim of time series analysis of gravity data is the calculation of the sample power spectrum and not the process power spectrum. In other words the mass density configuration that 41 10000 r rlTlT j 1 1000 llllll 1 T TITITI T 1 1 10 Cycles/Data Interval Figure 5: Direct and Indirect Spectral Estimates: Profile 2A 42 10000 [Tilt] 1000 10 1 1 l J 0.4 0.5 Cycles/Data Interval Figure 6: Direct and Indirect Spectral Estimates: Profile 20 10000 1000 10 1111'. 1 i ITIFIT 43 Figure 7: Direct and Cycles/Data Interval Indirect Spectral Estimates: Profile 48 ELECT to f-\ (“x "w 44 produced the individual gravity profile is under investigation and not the world wide mass density distribution. In this situation stationarity is not a necessary assumption. CHAPTER III THEORETICAL GRAVITY ANOMALIES A simple geometric body is often postulated to be the causative mass in the analysis of an isolated gravity anomaly. The parameters of this body are then estimated from the structure of the anomaly curve by equating it to the theoretical gravity attraction of the body. However, the mathematical complexity of the equations for the gravity attraction of even simple bodies makes parameter determination difficult. Recently, Odegard and Berg (1965) indicated that the Fourier transform of the theoretical gravity anomaly of several idealized bodies offered a simplified interpretation method. This was also shown to be true in the case of the magnetic field over idealized bodies (Bhattacharayya, 1966). However, the applicability of this technique to actual gravity data has never been determined. One major limitation in the application of this technique is the finite length of observed gravity data. Although field data can be accurately collected and reduced, and an assumption of a simple body as a causative mass justified, because the data are of finite length only a portion of the theoretical gravity anomaly can be considered. The entire effect of the postulated body is not observed due to the finite data length, and in some cases due to the influence of other geologic structures. To have practical usefulness a more general transform method based on theoretical anomalies of finite length must be studied. 45 46 In this chapter the Fourier transform method, as it applies to simple idealized bodies, is detailed. First, the frequency domain method is extended to include bodies that have not been previously studied. Then the interpretation of the transforms of theoretical gravity anomalies of finite length is examined. And finally the length of data required for the finite transform to approximate the ideal transform is determined. As in the discussion of gravity field data (Chapter II) the discussion is primarily concerned with one- dimensional data (i.e., profiles). Transformation of Theoretical Gravity Anomalies In the case of theoretical gravity anomalies the exact mathe- matical expression of the anomaly is known and the data is deterministic. The frequency domain is thus accessed by the Fourier transform of the original data and power spectra calculated by the direct method. There are several acceptable forms of Fourier transform pairs (e.g., Lee, 1960; Robinson, 1967). In this study all calculations and equations are based on the transform pair: N—ul =:| F(w) ( f(x)e"°‘*”‘dx f(x) F(w)eiwxdw , II ‘*‘\x 47 and for the two-dimensional case: on Hum-(21;: f [f(xwle’fiuxwfldxdy 11' -00 (X) 00 I / F(u,v)ei(ux+vy)dudv f (x.y) where 2nfx U V 2nfy . Review of Previous Work Odegard and Berg (1965) calculated directly the Fourier transform of theoretical gravity anomalies due to a sphere, fault, and a two- dimensional cylinder. The simplicity of the frequency domain expression is immediately apparent from a comparison of the equations in both domains as given in Table 4. The transform of the gravity anomaly due to a fault is not from Odegard and Berg because their calculation was based on the assumption that the anomaly equation is an odd function. However, it can be shown that the function is neither even nor odd. In the case of the two-dimensional cylinder, the transform as plotted on a log scale is a straight line with slope equal to the depth and intercept equal to the Gaussian mass. Interpretation of the other frequency domain expressions, although more complicated than that of the two-dimensional cylinder, are more managable than with the spatial domain formulas. 48 TABLE 4 GRAVITY ATTRACTION OF THREE IDEALIZED BODIES: SPATIAL AND FREQUENCY DOMAIN EQUATIONS SPATIAL DOMAIN FREQUENCY DOMAIN A)* Two-Dimensional Cylinder 2 DC 2 ‘DC/w/ 9(XI = ZWYRC o—-—-- G(w) = ”Y RC 96 0 +x2 c B)* Sphere _ 4 3 Ds 4 3 9(X) - gnv RS 0 (D 2+x2)3/2 G(w) - ngs owk1(st) C) Fault 2+ d2 g(x) = 2ret3g- 1n[ l + 1;- G(w) = 2Wp{nT6(w) x2+(d2 +T2 ) + d tan“ (g) — (am tan“ [d—é—TJ} + i(-—l- e'(D"T)‘° - J; e'Dwn (I) 00 where Dc = Depth to center of cylinder 0 = Depth to top of fault Rc = Radius of cylinder T = Throw DS = Depth to center of sphere y = Gravity constant R = Radius of sphere p = Density S * from Odegard and Berg (1965) 49 General Transformation The complexity of the mathematics involved in the transformation of theoretical gravity anomalies can be greatly reduced by use of a general transformation. The equation for a theoretical gravity anomaly over an idealized body is developed by solving a volume integral based on the summation of the attraction of all point masses in the body. These volume integrals are the partial derivative of the gravitational potential of the body. Bhattacharyya (1967) has shown that the potential expression can be directly transformed and a general transform for the volume integral obtained. The transformation of the theoretical gravity anomaly is then calculated by carrying out the integration in the frequency domain. When combined with another suggestion of Bhattacharyya (1966), one-dimensional transforms can be obtained over two- or three-dimensional bodies. If we consider the coordinate system where P(x,y,z) is a field point, and q(x1,y],21) is a body point. 50 Then the gravitational potential is (Grant and West, 1965), ] , U 3 ’ ) = 1 d d d o The two-dimensional Fourier transform of U(x,y,z) is co ( . ) = 19;— 1}. .jj:{’ 2 1 2 2 1 d d d U u v 4.. [(x-xp +(y-y11 + in(w) or 91(X) <=> 1wG(w) where f](x) = giggl— Then 2 2 l _ (x+b/2) - (x+b[2) + l D 22+(x+b/2) ( - 2 1n (2 2) g X) yp {02+(x+b/2)2 d2+(x+b/2)2 2 d 2+( x+b/2)2 _ (x-bm2 (x- b/2) _ 1.1n(b +(x- b/2) ) 02+(x-b/2) d +(x- 11/2)2 2 d +(x- b/2) 1 1 1 1 + - - + Dz+(x+b/2)2 02+(x-b/2)2 d2+(x+b/2)2 d2+(x—b/2)2 55 Now the expression is split up and the superposition principle of Fourier transforms employed f11x) + f21x) + f31x) = F1 e'iwcF(w) f(x+c) <=> eiwcF(m) F (w) = l {eiwb/Z [E_e-dw _ E_e-Dw] _ e-imb/Z [E_e-dw _ E_e-Dw]} l d D d D iwg_ -im§_ 1 2 n -dw n -Dw 2 n -dw £_ -Dw F2(w)-§-T?{e (‘58 -a)'e ]-e 5e -me ]} lw_l_)_ -l'u>p_ F300) _ %{e 2 g e"D(.U_ '3' e dw] e 2 [ID]; e"DU)_ Elie‘dMJ} and 1w G(w) = 2Y0 (F](w) + F2(w) + F3(w) ) but, F](w) + F3(w) = 0 so To) G(w) = :TYFB [(eiwb/Z _ e-lwb/Z) (% e-dw _ ge-Dw) ] 1w G(w) = $3 [2 i sin (”—2- (i:- e'dw - (lie-D”) ] so») = 1119 (sin 93-) (e'dw - e'D“) . (5) (A) General transformation Using the formula for the general transformation the mathe— matical ease afforded by this equation is indicated and the equivalence between the two methods demonstrated. Putting the limits of the two- dimensional prism into the integral formula for a profile along the x axis (Equation 2); 0° 0° 2 2 D b/2 -\/u+v z1 -i(ux]+vy]) G(u,o) = g1 I f f e e dxldzldyldv ... -.. -b/2 integrating over x1 , 2 2 D - \lu +v z1 -ivy1 I e e dzldy1dv d G(u,o) = %%E sin 9%- f j: 57 integrating over 2], 212. ub m m l -d u2+v2 -D u2+v2 -ivy1 G(u,o) = 21m $1" —-2— { JW [e - e ]e dy1dv integrating over y], 2 m 2 2 2 G(u,o) = g§§§3-sin 2%. I, 1 [e'd u +v - e'D u +V ] 6(v)dv Vu +v and finally integrating over v. G(u,o) = 232$ (sin 9%) [e’du - e‘Du] U which is the same result achieved in direct transformation of the gravity anomaly. Method of interpretation The frequency domain expression for a profile at right angles to a two-dimensional horizontal prism, i.e., G(w) = 3%2.51n Qg-[e-dw - e'Dw] U) is a less complicated expression than its spatial domain equivalent and theoretically can provide unique solutions to the body parameters. Zn The transformation will have null positions at multiples w = E—-. From these null positions the width of the body can be determined and the term gig-sin gg-divided out of equation 5, leaving 0) 2 walwb = p[e-dw _ e-Dw] . (5a) —-;Lsin-§- (L) 58 If the transform is examined at the high frequency end the curve will be composed entirely of the term efd”, because D is always greater than d and consequently e"Dw damps out more rapidly with frequency. This method is the basis for the multiple decay spectra technique discussed in Kaplan (l956) and illustrated in Odegard and Berg (l965) for the case of a fault. Basically, this technique involves fitting a straight line to the high frequency end of the function given by equation 5a. This line, representing pe'dw, is then extrapolated to the origin. Equation 5a is then subtracted from this line and the result equals e'Dw. Thus, the two depths are determined. Knowing these depths and the width permits determination of the density from the value of the spectrum at the origin where G(o) = 4nypg-[D-d]. An illustration of the various stages of decomposition of the spectra by this method of interpretation is shown in Figure 9. Thus from fairly simple mathematical manipulations of the theoretical transform the body parameters of the two-dimensional horizontal prism can be uniquely determined. Three-Dimensional Prism The general transformation formula permits the calculation of the transformation of the gravity anomaly due to a three-dimensional prism. This method is the only method available to obtain this expression because no general expression for the spatial domain is readily available in the literature. 59 ITlrj 0.1 \ l l l l l J> / A l l l J 0.1 0.2 0.3 0.4 0.5 Cycles/Km A: ———§$91——- = p[e'dw - e'Dw] (equation 5a) C: pe-dw _ {p[e‘dw _ e-DwJ} = pe-Dw Figure 9: An Example of the Method of Multiple Decay Spectra 60 Theoretical transf9:m_ The twoedimensional frequency transform provides a formula that is suitable for interpretation. However, the one-dimensional transform, although undoubtedly simpler in form than the spatial equivalent, is not readily interpretable. Substituting the body parameters into the general transformation formula, for a flat surface (2 = o) D b/Z a/2 1oz] -i(ux1+vy1) G(u,v) = Efi- { I. I, e e dy1dx1dz1 (6) d -b/2 -a/2 where a = width in the y direction. Integrating with respect to y]. D b/z -wZ 'X 1U = 1.9— l ' E 1 1 G(u,v) n (v) Sln 2 f I, e e dx1dz1 d -b/2 similarly for x], - z = gxg_ 1_ . va __ ub w l G(u,v) fl (v) Sln— (l ) sin 2 J’ e dz1 d and integrating over 21, we obtain G(u,v) = $19 (—1 sin 2—3 (L —) s n— ”U —) [e 'dw- e‘D‘i’J (7) 61 inhich is quite similar to the one-dimensional transform over a two- dimensional prism and can be interpreted in an analagous fashion. To calculate a one-dimensional transform for a profile along the x axis °° ( 2 2v; 2 21/; G(u,o) = 219- 1- sin 99- {e'd‘ u +V l e'D (u +v)) l .' va 'rr 2 —s1n TdV 2 2 '1; V (u +v1 2 2 v e( _ er . ub -d (u +v )‘ u,o) - nu Sln §—- { e (19 sin %9- dv '°° (u2+v2)"" v 00 D ( 2 + 2 )1]; - u v . - e 1 . va ‘[ -———————-— (V) Sln §—-dv } 2 Va ) -oo 2 (u +v The integrands are even functions, so; 00 2 ZIIL = 419 . ‘ub_ I -d (u +v ) G(u,o) nu Sln 2 { e (19 sin %9-dv 0 2.2V; V (u+v) 0° 2 2']; -D (u +v ) ' e (1) sin 29-dV} 0 2 2V; V 2 WW) It is not possible to obtain a closed solution for these integrals. However, Battacharayya (1966, p. l20, Equation 9A) gives a series solution for this equation. Using this solution the transform is . . 2 _ 4Y9 . up a. - 2mm 1 a u m , ‘5‘“) ' 'fil’ 5‘" ‘2‘ { 2 E (,2m+1):(-4B )- Km(“B-:1)- m=o ‘l a :E . a u m m=o where 2 2 I/ B1 = (d +a /‘*) a W II 2 2 . 2 (D+a/u)/l m order modified Bessel function of the second kind. 5‘ Although appearing complex it does seem to be simpler in form than the spatial domain equivalent. A hint of the complexity of the spatial domain expression is afforded by considering the restricted expression offered by Nagy (1966). Vertical Line Element Contrary to the previous bodies the simple expression of a vertical linear mass element is not readily improved by transformation and no advantage in interpretation is noted. If we have a vertical line element along the z axis, as shown in Figure 10, Figure l0. Vertical Line Element 63 Then, the gravity attraction along a profile in the x axis is l g(x) = Yo{ } (d1 +x {a ' (d22+x2ft The transformation is readily obtained by calculating the transform of this expression or using the general formula. In the direct Fourier transform of this expression we obtain, using F( 1 ) = Ko(aw) X2+a2 :1|—- 9(a)) = 9,91 [Ko(d1w) - Ko(dzwn (8) where Ko(x) is the zero-order. modified Bessel function of the second kind. In this case the parameters are not readily interpreted and only some approximation methods of G(w) offer any hope of interpretation. Effect of Finite Data Length on Spectral Analysis Interpretation To have practical application the interpretive value of the frequency domain method must be investigated for the case of finite length data. This can be accomplished by considering the transform of a segment of the theoretical gravity anomaly. Ideally interpretation methods similar to those for the transforms of the entire anomaly. as discussed in the previous section, should be developed. However, a 64 study of the theoretical transforms of the segments of the two! dimensional cylinder and prism anomalies suggests that solutions do not exist in closed analytical form, and therefore no simplified interpretation methods are possible. A practical alternative is afforded by determining the length of data required for the empirical transform of the segment to approximate the theoretical transform of the entire anomaly. For these investigations the anomalies due to the two-dimensional horizontal prism and cylinder are considered in detail. These particular bodies were chosen because of the simplicity of their theoretical transforms and associated interpretation methods. Calculation of Theoretical Anomalies To calculate the Fourier transform of a segment of an ideal anomaly either the transform of the signal composed of the theoretical anomaly multiplied by a data window is calculated or the theoretical transform is convolved with a spectral window (Chapter I). For ease in calculation, and to utilize the theoretical transforms already calculated, the latter approach is employed in this section. Spectral windows As discussed in Chapter I the data window plays an important role in practical spectral analysis. In general the selection of data windows based on some optimum criteria, valid for all spectral calcu- lations, is not possible (Jenkins and Watts, 1968). The best data window to be employed must be determined for each particular application. 65 When choosing a data window both the space and frequency domain characteristics must be considered. A good example of the tw0esided nature of a data window is the rectangular window: 0 x<-k h(x) = l -k$x‘k o x>k Although this window has ideal time domain characteristics, that is, it does not distort the data to any degree, the large side lobes in the Fourier transform of this window (spectral window); H(w) = Slgkwk =l|7¢ caused by the sharp truncation of the time function, generally produce poor spectral results. Four spectral windows, rectangular, Bartlett, Tukey, and Parzen, will be considered to study their possible effect on the results obtained in the calculation of the finite transformations. Three of these windows are related to each other by (Jenkins and Watts, l968). single" () ~]_ n 9 Hn(w) n .95 n where n = l = rectangular 2 = Bartlett 3 = Parzen 66 The effect of increasing n is to decrease the side lobes. However, then the spectral window becomes flatter and wider because the first zero occurs at w =,%_-. The Tukey window is described by sings ”W = k wk n _—‘]I€'2— —n- (14%) ) . For convenience in identifying it in this study, and from the relation- ship of its spectral characteristics, we will identify it by n = 3. The four spectral windows are shown in Figure ll. Two-dimensionalgylinder_ The finite transform for the two-dimensional cylinder is calculated by convolving the various spectral windows with the theoretical transform from Table 4. Rectangular window In the case of a rectangular data window _E SW km H](w) _ n kw em) ..9 e‘Dlwl where D = Dc the depth to the center of the cylinder 8 = 2nyRc2p y = gravity constant p = density R = radius of cylinder , 67 ”(w) 2K —- Rectangular ___.____.___. Tukey Bartlett .1\K _______ - Parzen ’1' Figure ll: Four Spectral Windows 68 and the transform of the finite length segment is T1(w) = [g e-Dlal _Ii Sln'k(w3‘0.) d0. n k w-a oo ___ _g i e-DIOLI Sln k(m-a) dOL 2n w-a multiplying by l = e'Dw+Dw T](w) = g—- e'Dw {I e-Dlal+Dw sin k(w-g1_ d " (w-a) a -00 let u = -w+a _ -D -D +u +D T](w) - g;- e w r. e (w I w sin k(-u2 du —m u case #l) w+u>o case #2) m+u"'w U<"(Ll Iw+u| = “+w Iu+w| = -u-w 69 con‘bining on . ‘ "w T](w) = ‘3"; e-Dw r e'Dw sin ku du + g}- e'Dw / e'D(-'u'2w)sin ku du -0) U -00, u = g— e'Dw / e'D‘” sin ku du + 3— e'D‘” I e-Dw sin ku du 2n --———- 2n -—~———- -w u o u '0) + %— eDw ‘ eDu Sln ku du N u 00 00 w in third integral let u = -u, invert limits, and I = I - f w o o o (I) T1(w)-%r— {eDusinku du+§—T?eD‘”{eDusmku du -w U 0 u 00 (A) + 9; e0“) / e'Du sin ku du - % ED“) / fi-Du sin k" d" o u o u o T1(w) = % tan-1 k [e-Dw+e+Dw] + % e'D‘” I e"Du sin ku du D u '0) (L) - £23— eD‘” e'Du sin ku du . (10) n u This equation is the transform of a finite segment of the gravity anomaly. The problem now is to put it in a more tractable form. 70 The main problem is the two exponential integrals. In general a closed analytical solution of these integrals is not possible. Even with an upper limit of the Nyquist frequency, a general solutiOn is not possible. These types of integrals have plagued electrical engineers in antenna theory problems and a volume of numerical solutions was published by Harvard University Computation Laboratory (1949). There are, however, two frequencies where solutions do exist 1) w: At w = o, the expression for T](w) becomes, (X) 71(0) =53 { e'Dlo‘lsingsg da -(X) mowfi ( e'Dam do. TT 0 0. w ems—e . That is at w = o the solution is the Laplace transform of the B spectral window multiplied by E-, or -l T](o) = %- tan (ll) ops- 71 This is also derived by substituting u,= 0 into the general expression. The result agrees with theory in the fact that as the data window become infinite in length the value of the transform of the segment equals the theoretical transform, that is, Lim T](_o) =g— = 6(0) . Using integrals 5.l.35 and 5.l.36 in Abramowitz and Stegun (l964) a solution at w = 1 can be presented. e-D D T](l) =9 e‘D + L IE](-D+ik) - g7 e 151mm (12) 2n where IE](x+iy) is the imaginary part of the exponential integral for complex argument as tabulated in N. B. 5. Applied math series No. 5l (1958). The value at this frequency also agrees with theory in that as the data window becomes infinite in length the transform of the segment of the anomaly equals the theoretical transform Lim 71(1) =-§— e'D = G(w) _ w-l k—mo 72 Higher order spectral windows From the relationship of the spectral windows, problems would be anticipated in the calculation of the finite transform using any higher order windows. The general transform solution would be on (sink(w-a) )n da _ Bk -D . .. Tn(w) - — I e la] “Way 411 -m n The solution would still contain the troublesome exponential integrals. To demonstrate this we can consider the finite transform using the Bartlett window (n = 2). In this case; -1 2 T2(w)=§—k-[%tan §-}oiog(i+fi,—)][e9w+e0w] 2 2 0 sin kw_ w sin Kg_ + i e'D‘” e'Dw 2 do) - i em [ e'D‘” 2 do) Trk 2 Wk 2 'w w 0 w The fact that at w = o the solution is the Laplace transform of the spectral window provides the solution at w = o for the Bartlett window but not for the Parzen window (Erdelyi, 1954). For the Bartlett window; -1 2 T2oo Two-dimensional_prism With the difficulty encountered in an attempt to calculate a general transformation for a finite segment of the two-dimensional cylinder, a similar attempt for the two-dimensional prism is probably futile. The mathematical structure of the two-dimensional prism. although similar to that of the two-dimensional cylinder, is more complex. An idea of the complexity can be determined by considering the convolution integral for the two-dimensional prism and the lowest order spectral window - the rectangular window G(w) = g)9-sin g2. (e-dw - e-Dw) w (0.)) = .219 Sin at) (l) (e'd|a|_ e'Dlal) Sin kéw-gl dd w- The same exponential integrals are encountered and furthermore no simple solution at w = o is realized. Thus, as is the case with the two-dimensional cylinder, no solution in closed analytical form for the transform of a finite segment of the ideal gravity anomaly is possible. 74 General considerations Vfi‘vr v‘ To further investigate the problem of determining a general solution two alternative approaches are considered. General transformation formula The general formula for the two-dimensional Fourier transform can be expanded to include finite length data. If H(u,v) is the transform of the two-dimensional data window that is applied to the data, then the general transform for finite length data observed on the plane (2 = o) is; T(U.V) = G(u,v) * H(U.v) -Z\(a2+82 'l'(onX1+B.Y1) T(u,v) = %%- e e H((v-B), (u-a))dx]dy1dzldad6. For the case of a profile along the x axis, the data window in the y direction can be assumed to be a constant equal to one with a spectral window equal to a delta function thus, 2 2 ‘2] U0. +8 ”-i(GX'I'l'BY‘I) T(u,o) = %%- e e H(u-a) dx1dy1d21dad8 . This equation then represents the one-dimensional general transform for a finite length segment of a theoretical gravity anomaly. Substituting the equation for the spectral windows, yields; 2 2 -21 \/a +8 -1' (ax-I+By]) sin k w-a n Tn(u,o) = (2)—TY; If e e n dx1dy1d21dod8. k(m-a) n 75 The major value to a general expression such as this is that usually by selecting the proper order of integration some mathematical ease is attained. However, in this case, the most tractable results occurred when the order of integration was x], y], z], B and then a and the result is exactly equivalent to the previous equations for G(w) * Hn(m). Direct transformation A second alternative is to compute the direct transformation of the signal composed of the gravity anomaly multiplied by a data window. To investigate this, the simplest case, that of the two-dimensional cylinder multiplied by the rectangular data window was considered. In this case, k D T (w) = Q- I, coswxdx . l n Dz+x2 0 Although appearing simple enough all attempts at obtaining a solution in closed form proved futile. Length of Data Required to Approximate Theoretical Transform An alternative approach to the problem of finite length data is to determine the length of data required for the transform of a segment of the anomaly to closely approximate the theoretical transform. The term "close approximation" being defined by high correlation between the two spectra and/or low resultant errors in interpretation. In their discussion on the use of the theoretical transform in gravity data analysis, Odegard and Berg (1965) indicated the excellent 76 (agreement between the actual transform of a theoretical anomaly of finite length and the theoretical transform of a two-dimensional “ cylinder. However, the gravity profile transformed was calculated to a distance where the gravity attraction was 0.0l maximum value. From a consideration of the formula for the gravityattraction of a two- dimensional cylinder, this implies a profile that has a length equal to 20 times the depth of the cylinder. In their specific example of a 3 km deep cylinder, this would demand a 60 km profile free from disturbing anomalies, an unlikely situation. Spector (1968) considered finite length in his study of spectral analysis of magnetic data by examining the individual factors in a general spectral equation and comparing a numerical solution of the convolution integral, using a Hanning window, with the ideal Spectrum. However, by analyzing individual spectral factors it is not clear how they effect one another and in this approach the effect of the numerical integration approximation should be considered. The method used in this section is to calculate the transforms of segments of theoretical gravity anomalies and compare them with the associated theoretical transforms. In this manner the results are directly comparable to transformations of ideal field anomalies. Also, no a_prigri_choice of data window is made and the four data windows, previously discussed, are studied to determine the utility of each window. 77 The two transforms are compared by considering the method of interpretation used with the theoretical transform. Statistical I measures of the similarity between the two transforms are calculated and errors in interpretation determined. Twoédimensionalgyljnder_ In comparing the finite transform to the theoretical transform the factor k/D, where k is half profile length (half—width of data window) and D is the depth to the center of the cylinder, is critical. The approximation of the finite transform to the theoretical transform is independent of the actual depth and dependent only on the ratio k/D. This is evident from an inspection of the finite length transform calculated in the previous section and was verified empirically. The two characteristics of the spectra that are used in interpretation are the slope and d. c. value (intercept). When plotted on semi-log paper the spectrum is a straight line with slope equal to the depth and intercept equal to one half the Gaussian mass. Slope Transformations of various length segments of the calculated gravity anomaly due to a two-dimensional cylinder, multiplied by the four data windows, were calculated and compared to the associated theoretical transforms. The linear correlation coefficient between the calculated and theoretical transforms, for the four data windows, as a function of k/D is shown in Figure 12. This figure demonstrates the relative merits of the four data windows and reveals for the rectangular data window a 78 sauce: Fmgpumam ucm o\x mo cowpuczm m we mseoemcmgh _muwpmgomgk use cape—aupeu :mmzpmm pcmePmmmou compmpmegou ”emuchaQ chowmcmewouozp ”NF mgzmwe o\x 523 llllll \\ .\ 36.3. \ .“ 332$ llil.) Lm—zmcmpomm \\\\.\ .\ N.o v.0 m6 1ua;o;4;aog uoiaelauaoo 79 very high correlation (0.972) at k/D = 3. All data windows produce a correlation coefficient greater than D.BD at k/D Z 5. The surprising fact that the rectangular data window provided the best estimates was further substantiated by considering other statistical measures such as standard error of fit, deviation factor, and regression line parameters. To determine effect on interpretation, depths were calculated from all calculated transforms by determining the slope of the natural logarithm of the transform (the point w = 0 not considered). The calculation for a normalized depth (i.e., calculated depth/true depth) as a function of k/D is shown in Figure 13. This graph again shows the relative merits of the data windows and indicates that less than 10% error in depth is obtained with the rectangular window at k/D 3 3. Less than 10% is obtained with the Tukey window at k/D 2 4. The Bartlett and Parzen windows produced poor results for this particular application. Intercept The effect of profile length on the intercept can be determined by using the equations for T1(o) and T2(o) (Equations ll,l3) determined in the previous section. = E. -l!s T](o) fl tan D (ll) T2(o) = étan 1';— g—fl {2— ln (1 + (k/D)2) (13) 80 zoucwz Fmguuwam ucm o\x mo coeguczu m we spawn uwNwfimELoz ”cmucw—xo chowmcmswoquH ”mp mgzmwu o\x n o m e m N P _ _ s s _ _ _ o .1N.o :mNcoa llllll amxzk \\ ppm_pemm ....:|||u..:u \\ cm_:mcmpumm \ L¢.o \ \ _\Ak. \ . video PQZIlPNJON 81 OY‘ .. 2 T2(o) = g-. (g-tan I g-- %—%—ln (1 + (k/D) ) and if [(g-tan'] k/D)] in the case of the rectangular window and [fitan’1 %-- %—E-ln (l + (k/D)2)] for the Bartlett window are plotted, then as these functions approach unity the finite transforms approach the theoretical transform. This approach is equivalent to dividing the calculated value at (n: o by the theoretical value and this was used to verify the determination of the equations for T1(o) and T2(o). Because no analytical expression could be obtained for T3(o) or T4(o), this method of dividing the calculated value by the theoretical value was employed to graph the effect of profile length for these windows. These functions are plotted in Figure 14. It is interesting to note that once again the rectangular window appears to be the most efficient. To achieve an error of less than 10% with the rectangular window a sample with a k/D ratio > 6 is required. Although this is twice as long as the sample required for the depth determination there is an alternative approach available. In the case of the rectangular window the slope can be used in the T1(o) formula to calculate the Gaussian mass. If the sample profile has a k/D ratio 3 3, then the 10% error in determining the depth will result in less than a l0% error in determining the Gaussian mass by T1(0)n tan"1 k/D 82 / _.—-O-' Rectangular ———o————u-——- Bartlett Tukey Parzen 1.0 F 0.8 — 0.6 F 0 4 Figure 14: Two-Dimensional Cylinder: Tn(o)/G(o) as a Function of k/D 83 At k/D 3 3, the arctangent fUnction is changing quite slowly and an error of 10% in determining D will result in a less than 10% error in detennining B. For example, a 10% error in D will result in approximately a l0% error in k/D (actually 9 - 11%) so for k/D of 3, the range would be 2.7 to 3.3 The resultant range in arctangent function is l.22 to 1.28 or about 3% error. Profile length The requirement that the profile length be equal to six times the depth (an unknown, to be determined) can be expressed in a more practical form. Setting x = 30C in the equation for the gravity attraction due to a two-dimensional cylinder (Table 4, Formula A), shows this to be the distance at which the anomaly has decreased to ten percent of its maximum value. Two-Dimensionaljrism In the interpretation of the two-dimensional prism transform the null positions are used to determine the width of the body and after removal of this factor, the remaining spectrum is interpreted by the method of multiple decay spectra. Null positions To study the effect of finite length on null determination a number of theoretical anomalies over ideal prisms of various widths, depths, and depth extent were transformed and compared to the theoretical transform. In this phase of the analysis only the null positions were considered to be of interest. The Bartlett and Parzen windows proved to be quite ineffective in this study and only the results of the rectangular and Tukey windows are detailed. 84 Initially the width was determined from the first null of the calculated spectra and compared to the true width. The results for the two windows as a fUnction of 2k/b (profile length/width; equal to k/b/2) are shown in Figure l5. This plot indicates that with either window an error of less than 10% in determining the width from the first null position is obtained for a 2k/b > 6. However, it must be considered that the accurate determination of the first null is affected by spectral resolution. Attempts at spectral interpolation by adding zeros to the gravity profile (with and without zero mean) did not improve the determination of this null position. However, the resolution in determining the width is enchanced by considering two (or more) null points. There was no error in determining the width using the Tukey window at 2k/b Z 2 and with the rectangular window 2k/b 2 3. Thus in this application the Tukey window proved to be more effective and a relatively short profile, of length equal to twice the width provides an accurate estimate of the width . Depth In the analysis of the spectra for null positions there was no consideration of the overall_goodness of fit of the calculated spectrum to the true spectrum. However, the interpretive goal is not only to determine the width but to then utilize that to factor out of the spectrum the effect of width and determine the depths. To accomplish this the entire spectrum must be considered. 2k/b 20 l8 l6 l4 l2 10 85 F F l l L l 50 40 30 20 10 Percent Error in b (first null) Figure 15: Percent Error in Width as a Function of 2k/b and Spectral Window 86 From a consideration of the remaining terms in the spectra, i.e., (e'dw-e'uw) it is anticipated that a result in terms of k/D similar to that for the two-dimensional cylinder would be obtained. If the terms could be considered independently, a profile with k/d 3 3 would be required to determine d to a 10% accuracy and a profile with k/D 2 3 required for the determination of D. In making the comparison between the spectral estimate and the true spectrum a linear correlation coefficient for the entire spectrum was calculated. That is, the width term was not factored out of each spectrum and the method of multiple decay spectra was not employed. From a comparison of the correlation coefficient and normalized depth of the two-dimensional cylinder the 10% error figure is comparable to a correlation coefficient of approximately 0.95. The correlation coefficient between spectral estimates calculated from bodies of various widths and depths, and associated theoretical spectra as a function of k/d is shown in Figure 16. A best fit line was drawn through the data for each window. Again, as in the case of the two-dimensional cylinder, the rectangular data window proved to be superior and a correlation coefficient greater than 0.9 is obtained at k/d > 3. The scatter of data points observed on this plot was thought to be due to bodies of varying depth extent. To examine this the correlation coefficients obtained with the rectangular window were plotted as a function of k/d for various d/D ratios (Figure 17). 87 sauce: Pecuumgm use U\x mo cowpuczu m we mseommcech _eowpmcomzh use nope—zopeu :mmzpmm pcmwuwwemou cowpmpmecou "amps; .mcowmcwswo-03h "or mczmwm e\x m m o m e m N P fl _ _ _ _ 4 _ _ .\ «\\ 4 \d \d \ O _\ s\ _\ a, . a . \\\ . \ d x. . \¢ ¢ 4\\ < .\\u . d <\\ o :mNmeflE. . . . x. 1\< o \x e . . . . \u cmpsmccpumm . \.¢\ 0 0 1111? d 00 I I. O H 111141111411 < o o. . N.o ¢.o 0.0 w.o o.~ auatogagaog uoiqelaauog 88 o\u wee u\x co cowpucsm e we msgoemcceh _cuvumgomzh ace Azoccwz ceF=mcmauwav umumpzupmu :mmzpmm pcmqumwmou cowpmpmccoo "Emwcm Pacommcmewo-ozh “NF mgamwu u). u m m N F _ _ _ 4 - o .4 N.o 1 ed p L 0.0 \\0 o ....o ....o /p .0 s0 0» 0 ( 1 m.o 1111“ W m L O._. quaioigyaog uoizelaaaog 89 This shows the dependence on the depth extent of the body. Some additional scatter was still observed that is attributable to various Zk/b ratios but it was not significant. In general all these points were calculated using profiles of 2k/b Z 1. Figure 17 shows that increasing length of profiles are required for bodies of greater depth extent. Basically by considering the correlation coefficient over the entire spectrum we are requiring that both the d and D slopes be fit accurately. To determine only the depth to the top accurately a profile with k/d Z 3 would probably suffice. To determine the bottom depth (and subsequently both depths) a profile with k/D Z 3 would be required. This is demonstrated by Figure 18 showing the correlation coefficient versus k/D. Profile length For determination of all the parameters of the two-dimensional prism a profile of length (2k) six times b/2 or D, whichever is greater, would be required. In the case of the two-dimensional cylinder the length requirement was translated into terms of the maximum anomaly. Although no such general figure for the two—dimensional prism is obtainable, an idea of the range of values can be obtained by substituting x = Cb, D = Nb, d = Mb, (N>M) where C is some constant into the gravity formula. The expression then becomes 90 Q\v_ ....o £05.0ch n ma mELommcmLk Pmumvmsowch Ucm AzoucP3 empamcmpummv umpm_:u—mu :mmzpmm pcmwuwwwmou cowuepmccou ”Ewes; Pocowmcmewo-ozp ”mp weaned Q\x w n o m ¢ m N F A _ _ d _ _ d O 1 N6 H. 1 so i 0.0 .. 1 m6 quaioi;gaog uoiseladdog 91 2 2 G(Cb) =y2pb[(c+i/2)'niJGLiigiifillfl- - -M2+c +C+1/4 . 2 2 N +C -C+1[fl_ — (C+1/2) ln -M2+C2-C+l/4 + N (Tan'1 1Efil£§l_ - Tan-1 (C’1/21) N - M (tanu1 FC+1 2 - tan'1 C'] 2 )]. 'Ihe profile length at which less than ten percent error in all parameters can be expected was determined to be six times max (b/2,D). If C is set equal to 3, then M xj for ioux JJam >426 N. e; um .omm: mac memu» zpom «a .m.m mH mudcz»~uazuur<»m mum oupmme zur» mu mmfiamm audazcm mt» .muzmozmumozu mmaamm< mute .o<4 mute rm omeumrz~ e< umsuzcm 2mm» m. mmsuem m1» .ouzmsamemo mm cum» 0» exam 0» uuc mce mom scammmuuz macs mu munzsz mrp oz< omh4<2< QZ< FZWZUEDmIQOUUIII>aIWUHm>LQ0m0III>aL wLP mo kmwh Jam 00.3 a“ mhmwk 03h wwwD zP~¢<20Hhh~a<2u~h~k«z anb..x«.m~..mu mwnaum bauzu m:# to Ihuzuau MI» 02¢ m Ikosz no wJazh~uh~u<2u~r4nz~..\..4:— mnxhliimu<4..e~..za oum~ Oh oz\.a~2\24m<.u.uu.m .auawbzu muzwoum4ou bzuuawn..x~.ucnm.xmvp..xm.rotzx.tazo244.m4xx.u()x 4440 Jfixuxzx uuzzz .xxuhmuhuu.xx.hmnpz n.unxx Nm uo «Nomwo.u m431 mu oUZ..XNonomu..u mu cab 00 00w nnw hww Oww own NM Jma m0 oOZ..XN.M~.-N N\Z.oxmonuo.fl 4.ox0..999m£.43x.azxo£zz.szvm<>x JJu.rzu~.ueézx.zazo£ZZ.XBva<>x JJ¢U .m14.+nuzzz fiwzzz «opzuuupzuu 4.~uz.unfl oc uo (Oxzumud onpzuu .\\\..*maannaehmwh 40mmw~o uz< hmmthH wUZmuanOU PtmuuuaooXA.«ammoXeru<£u~h.zxx.;¢z.442.x.m<>x wzthOume CNN ¢¢N 114 c..on..mm..hn..om..mn..cn..nn..~m...m..on..mm..sm..o~..m~..om..nw.. ONNOOONooOuoofluoofiuoovuoomuooMnOONuoouu000uoomuOF.CO.ofioonoomuuoofit «a..n—...«a...o«~..mo~..mo~..~ou..oo~..mou..eou..no~..~ou..~ou..uo. ~..om..mm..sm..mo..cm..no..mo..no..oo..mn..mm..sm..cm..mm..¢m..mm.e ONwooOw0005.095oobfiooOFoomfioochoonhoONFoouho00h.0000030.OOOoOWUOOQQ OOOHUOONU concoonvOooOMOQ mflooomoommoo¢WOONOo Owao ouflo 00m. 00¢. 000.003.. .me..cc..nc..~e..ec..cc..mn..mm..sn..mn..em..nn..mn..nm..ow..e~..s. NOOONoomNoonNOONN couWooONooOuoOQuooOnoon“. Gown. 0"". ONdo CO“. 00. QWOOOQ ..u..n\~x ~Jum ..uvnx.amoaouvxwiUZQJ<)~30m ..uVNx.anooonuxwiUZ¢Jdauzuu «guy—x..—.00u.xx.wuzmq<>~num aoomamxx.Aoowvmxx..oomvuxx.aoomvnx..oow.~x..oom.«x zunmzwzno ao.oouvxxx..o.00u~xx zuuwzuzuo .«vx zu~mzmxnu anvhmubu zuumezuo . «hmwhu.ZDuobduo..omo..mmo..ooc..onm..oom..mmm..m«m..mm. c..moc..~cc..mwc..moe..mmn..osn..omm..oem..own..mom..mcm..o»~..mm~. ..ocw..mw...o.~..mm~..~mu..-~..~c...oc...ocn..omu..om~..o~"..oo~.. .Nu.oww ..WFOOQO ..00..°w..mc..mn..dn.CNN. .NN. .0“. Om". ON“. ON. ON..m. 0' e..n..w..u..~..«..mmum..om-..oauw..omom..ooom..oco~..mom~..osww... onmu..oom...ome_..o~e~..ooog..oeon..ooou..amm...m—mu..moed..oee»... mac—..wmm~..~em~..o~nu..mmwu..m¢-..om~...om”u..omfin..omuu..omo~... moon..m"om..ooo...oem..oem..ouo..omm..omo..omm..oom..cms..-e..o-. ..moc..mco..o~o..ouo..ooc..orm..omm..c_m.."om..mme..mme..mee..mue.. .uocn.omm..omn..men..own..mon..om~..okm..mm~..oc~..kmw..m-..oo~... emu..ue...ocu..om~..ocu..omu..o~...o-..oo—..oo..om..mo..~c..em..e. cOOdq..mn..°m..mN.. WW..QN ..nd .OO“OOQ..O. 0m. 0‘. On. ON. C“. .d\~xx (F(o \.m«~..~"~..o~«..m6n. ..mo~..so~..oo~ ..mOd OOVOH ..now..~on.. HONOOOO". 0J1”. .0“. 05¢. .0W...nwu.. CH“ ..Ww ..“W .00” ..ow ..wm ..FQOOOQ..mm..¢Q. One. 0“”. 0””. Com. 0“”. Ch”. 0”, F.0mh..¢fi..flh..~h..flh..°b..ou..mo..ho..00. Ono. .«o. .m0. .“0. 000.00“.. .wm..hm 0.0“ ..mm ..¢m..mm ..dm..om.OOQOCQQ..Nc. 00*. 019‘. Cowc. Ono. .fi'. .0, 116 plunpzu< .moaam rh~r<2xum .nom.n.mp~u‘ vac op uo N\24uca uaznrzuu nne o» oo conu.fifiux Oo«¢ZJOntaa n~¢.¢~e.ooo.ooniaqiwa.x.uu nnc Or 00 cocn.1w.x o.«+zaxu44a o~¢.-e.ooo.o.n1.~ina.x.uu ouc up ou.o.o.p4..11.x.u~ axl.ww.xn.wfi.x u41.1zu11 nnc on «+4241u21nzz .oo~nnzz.oom.»o.azzvmu «ozzzuaz ugznpzuu cunu.?zz.x do.» 0.. Do coeu.zzz.x moe.~o¢.noc.zx1.zzz.x.u~ Oouuzaa uncurzuu w42~p2uu oo~+h‘uu»¢u.auowux.po..~.x.un Joana unv uo 4.4zzwu one on «1132»; uooupzu «+4zzlzzznz .oo«n:az.oo—.»u.zsz.u~ sexuozz \oumNN..omo~..ooo~.comm»..oomm..o omom..oo¢~..omn~..oomm..oowm..mm~m..om-..oo-..mmom..moom..ommu... owuu..osmn..onww..oomu..shhu..mosu..omo«..omo~..oooa..omm«..oumu... awe"..occu..ooeu..00nu..o«m«..mo~_..ommu..mo-..omflu..o¢~u..mo-... amou.ou¢0u.onuON.oooo.oono.oooo..onm..oam.ooum.comb..oos..om»..o~ha .ousu.pomu.oono.oooo.oosm.oomm..onm.ooom.oooo.oooc.coaceomue..oon.o omnn.oomm.oomm..omn.ooom..omw.omom.oomw..mnm..o-..oow.covaooouu..o mun..0m«.omn~..ow~..o~«..mo..om..mh.cno..oooomm..cc..on.o~n.omw.oua N.ohu.onfi.oon.om.fio.o¢..m.odomN.o0uhN.commN.oooom.oommN..oomN.oow¢n ¢¢Q now can find QNQ MNG mac ONc uNe one QUQ moo Noc One «MC 117 02m (makmz mJturZOU awe «+a 1..." vFWUhu "~13”ahmmhu~Z~ ¢Om N~ zmom omXUJuuhdba 0h xuuhh~mam nemeeamao..xm~.\\\\.p..mu.mu.m..mu.3..mu.o..mu.m«.~..m~vsz..mn.ma.x zoumzuzua .mu.mu.-..m-.—-..m-.-..Muvgs..mu.ou..mu.o zouwzusno nu no .Ouem mm 2 huMHZ~ 104m 4uIb~3 m4u~h<>uww33 omuJJawh2~ x bumdflm .OuwN Uh O¢quu Uéu> mu~m>Iauw0 000(00 I OPOHOZ >m<2u~hL th mu hmmh Jzmooefi mom I I I n Ubum to mom on 20 Q Qwhm #4 s02 wa< Ms an+aciufievuu~n.s.~_- om ulfiuufi z.uuu 0m 00 x.una cm 00 - x~¢hz~ «wuwz4~ JP4JJUJ¢U sOZ afiaquuau+.29ufiv~NN om UUUU UUU UUUUU UUUU UUU 120 malanUU o—m ...ouczamuonzamuc .u.0323wouzamc ...uu— 03 ac Ooouzamoo coowzswo £w£< xu¥< MJ‘NPZUU 00m aa~v0.0u0044n.~.00 oou 000 Oh uo a\\..0whuum< hwwh >h~u m>~hoafio~VX¢2JmonaJwO hon ¢.un~ bow on coonaamo x.uwa com 00 O I I m MP(JJUJnafi.uvr .1.xx~xxax¥.uvmczuwruznwr mam (.«nxx mom on coougawr 4.~u~ com 00 . xeuuw can on UUUUU UU UUUU 121 «xuezu aaqama» mt» zh~u ucdomae o» 44:06 no car» mm»; mu .upmue.«4 en numoa.m w4m ukauuauuuda rpm: pu muum z. «4 ecu.»m~»<»m emu» ux» mu «pmue zsoxsuzu(40~»ath~ uuzwuuucou h4MUamu..x~.nomm.xmvpu14u~huth~ UUZMUUULOU h‘wuauu..xu.domuexmvr~WDJU£UU£H mthP UMPDuzuuooxno\uh<2mcu New . .Nuo.mUM» an; oooH Uh 00ANuokaououzh~u<4uuhm 4uHmum>4~ xwah .m.<.zvm4~