Download Power Spectral Densities and Cross-Spectra of Gravity Gradient Components - Prof. Xiaobing and more Papers Engineering in PDF only on Docsity!
GEOPHYSICS, VOL. 71, NO. 1 (JANUARY-FEBRUARY 2006); P. J11 – J22, 9 FIGS., 1 TABLE.
10.1190/1.
Spectral analysis of gravity gradiometry profiles
James While 1 , Andrew Jackson 1 , Dirk Smit^2 , and Ed Biegert^3
ABSTRACT
The gravity gradient tensor (whose components are the second derivatives of the gravitational potential) is a symmetric tensor that, ignoring the constraint im- posed by Laplace’s equation, contains only six inde- pendent components. When measured on a horizontal plane, these components generate, in the spectral do- main, six power spectral densities (PSDs) and fifteen cross-spectra. The cross-spectra can be split into two groups: a real group and a pure imaginary group. If the source distribution is statistically stationary, 1D spectra can be found from the 2D spectra via the slice theorem. The PSDs form two power-sum rules that link all gradi- ent components. The power-sum rules, in combination with further equalities between the power and cross- spectra, reduce the number of independent spectra to 13, a number reduced to seven if the power spectrum of the potential is assumed isotropic. The power-sum rules, cross-spectral phases, and coherence between compo- nents all provide information on the internal consis- tency of a set of gradiometry measurements. This in- formation can be used to assess the noise, to determine the isotropy, and, for a self-similar source, to calculate the scaling factor and average depth. When applied to a data set collected in the North Sea, the power-sum rules reveal high-frequency noise that is distributed among only three of the gradient components; additionally, the coherences reveal the source to be anisotropic with a nonzero correlation length.
INTRODUCTION
Measurements of the second derivatives of the earth’s gravi- tational potential, or gravity gradiometry, have long been part of geophysical exploration. In the early 20th century, Eotv ¨os torsion balances were widely used to measure the horizontal
Manuscript received by the Editor June 8, 2004; revised manuscript received April 25, 2005; published online January 12, 2006. 1 University of Leeds, School of Earth and the Environment: Earth Sciences, Leeds, LS2 9JT U. K. E-mail: jwhile@earth.leeds.ac.uk; jack- son@earth.leeds.ac.uk. (^2) Shell International Exploration and Production Inc. Volmerlaan 8, P.O. Box 60, Rijswijk, 2280 AB, Netherlands. E-mail: dirk.smit@shell.com. (^3) Shell Oil Company, 3737 Bellaire Blvd., Houston, Texas 77025. E-mail: ed.biegert@shell.com. © c 2006 Society of Exploration Geophysicists. All rights reserved.
gravity gradient before being superseded by gravimeters in the 1930s (Telford et al., 1990). More recently, practical gradiome- ter systems capable of rapidly measuring all components of the gravity gradient tensor have been developed (Jekeli, 1993; Bell et al., 1997). The great interest in gradiometer systems stems principally from the insensitivity of gradiometers to linear ac- celerations. This insensitivity makes them ideal for use in mov- ing environments such as on a ship or aircraft, where the in- evitable accelerations of the platform make vertical gravity measurements difficult and error prone. A review of the char- acteristics of gradiometry can be found in Bell et al. (1997). In total, nine gradients form the elements of the gravity gra- dient tensor Γ:
U xx U xy U xz
U yx U yy U yz
U zx U zy U zz
,^ (1)
where U is the gravitational potential, and the subscripts re- fer to differentiation with respect to the Cartesian axes. In free space the gravitational potential obeys Laplace’s equa- tion and the gradient tensor is traceless. Also, as a result of the commutativity of the differential operators, the gradient tensor is symmetric. Because of these properties, the gradi- ent tensor contains only five independent components. How- ever, measurement noise will result in the apparent violation of Laplace’s equation, so we treat Uzz as a sixth independent component. Further details on the gradient tensor are given in Pedersen and Rasmussen (1990). All current commercial gradiometers are based on a design that uses opposing pairs of accelerometers set on a rotating disk. By combining the rotation of the disk and measurements of the differences in acceleration between the accelerometers, we can measure one gradient component and the difference between two of the other gradient components (Lee, 2001). A minimum of three disks set at an angle to each other (Brze- zowski and Heller, 1988) are required to measure the full gra- dient tensor; this results in a complex system with a minimum of 15 parameters (six accelerations and nine angles). For such
J
J12 While et al.
a complicated system we ask, “How consistent are the mea- sured gradient components with each other, considering that all of the components are bound together by the form of the potential?” This paper presents relationships that exist in the spectral domain and that provide a means of answering this question. The theory section details the formalism that is at the heart of our paper — namely, the transformation of Γ into the spec- tral domain. In the spectral domain the six components of the gradient tensor generate six power spectral densities (PSDs) and 15 cross-spectra, which form the basis of our data analysis. Using the derived form of the PSDs, we show that two rules relate the Uzz component to the other gradient components via simple addition; these rules are akin to the power-sum rule given for the magnetic case by Parker and O’Brien (1997). In addition, we show that six further equalities exist between the different spectra and that all of the cross-spectra are either real or pure imaginary. These relationships are given for both 2D gridded data and, assuming statistical stationarity of the source distribution, 1D line data. If we additionally assume isotropy, the 1D spectra simplify considerably. We demonstrate the spectral relationships using two differ- ent numerical realizations of stochastic models, characterized by the value of a scaling constant β. One model is based on a source distribution that is uncorrelated on all length scales (β = 0), and another model is based on a self-similar source distribution (Maus and Dimri, 1994) where β = 3. The power- sum rules are corroborated in both simulations, and the co- herences of the cross-spectra are shown to agree with theoret- ical values. However, the difference between the two models shows that increasing the scaling constant steepens the decay of the PSDs. More interestingly, as the scaling constant be- comes large, the coherences between the different gradient components smooth out so they are constant with a value of 1 or 1/3 over the entire spectral domain. Last, we show the spectral results of a shipborne gradient survey conducted over a small region of the North Sea. The form of the derived spectra means that the potential cannot be self-similar or isotropic and that either the source distribution cannot be considered statistically stationary or that there must be significant error present in the measured gradient compo- nents. In particular, the power-sum rules indicate that error is not evenly distributed among the components.
THEORY
We stated earlier that the gradient tensor is symmetric, and we derive here the spectral characteristics of its six indepen- dent components. We do not use the traceless property of the gradient tensor because noise in real data will mean that Laplace’s equation will not be exactly obeyed, and it is useful to have an analytic expression for the PSD of the Uzz com- ponent. We work on a flat earth using Cartesian coordinates (x, y, z) where the positive z-direction points downward into the earth. We define the 2D Fourier transform as
−∞
−∞
q(x, y)e−^2 πi(k^ x^ x+k^ y^ y)dk x dk y ,
where kx and ky , the wavenumbers in the x- and y-directions, are the components of the wave vector k.
The upward-continuation formula between a surface at height z = 0 and another surface at height z = h can be written as the convolution (Blakely, 1996)
U (x, y, h) = ν(x, y, h) ∗ U (x, y, 0), (3)
where
ν(x, y, z) =
z
2 π(x^2 + y^2 + z^2 ) 3 /^2
We now use the convolution theorem to write the upward- continuation formula, equation 3, as
This equation is valid for any arbitrary field; however, we concentrate on second-order statistically stationary fields for which U (x, y, 0) is treated as a random variable with a mean and variance independent of x and y. The appropri- ate spectral-domain representation for a statistically station- ary variable is the PSD (SX), defined as the Fourier transform of the autocovariance function of the stationary random function X (Percival and Walden, 1993). We rewrite the Fourier-domain convolution formula, equation 5, as the PSD equation
The Fourier transform of ν (Blakely, 1996) is
F [ν(x, y, h)] = e−^2 π| k |h^ , (7)
where | k | is the magnitude of the wave vector | k | =
k^2 x + k^2 y. By substituting the transform of ν (equation 7) into the power spectral form, equation 6, and differentiating, the six PSDs are found to be
SU xx (x, y, h) = 16 π^4 k^4 x e−^4 π| k |h^ S 0 (k x , k y ), (8)
SU yy (x, y, h) = 16 π^4 k^4 y e−^4 π| k |h^ S 0 (k x , k y ), (9)
SU zz(x, y, h) = 16 π^4 | k |^4 e−^4 π| k |h^ S 0 (k x , k y ), (10)
SU xy (x, y, h) = 16 π^4 k^2 x k^2 y e−^4 π| k |h^ S 0 (k x , k y ), (11)
SU xz(x, y, h) = 16 π^4 k^2 x | k |^2 e−^4 π| k |h^ S 0 (k x , k y ), (12)
SU yz(x, y, h) = 16 π^4 k^2 y | k |^2 e−^4 π| k |h^ S 0 (k x , k y ), (13)
where S 0 (kx , ky ) = SU (x, y, 0) is the PSD of the potential field at h = 0. Because | k |^2 = k^2 x + k^2 y , we see that equations 8– 13 are related to each other by
which are the 2D power-sum rules. Parker and O’Brien (1997) show that for a statistically stationary field, the 1D PSD P X of a single data line
J14 While et al.
From inspection of the cross-spectra and PSD equations, we see that
P [U xy ](k x ) = P×[U xx , U yy ](k x ), (24)
P [U xz](k x ) = −P×[U xx , U zz](k x ), (25)
P [U yz](k x ) = −P×[U yy , U zz](k x ), (26)
P×[U yy , U xz](k x ) = P×[U xy , U yz](k x ), (27)
P×[U xx , U yz](k x ) = P×[U xy , U xz](k x ), (28)
P×[U zz , U xy ](k x ) = −P×[U xz , U yz](k x ). (29)
These relationships combined with the power-sum rules (equation 22) mean that only 13 independent gradient spec- tra exist. Additionally, we see from the cross-spectral equa- tions that the cross-spectra are all either real or pure imagi- nary (see Table 1). If we write the cross-spectra in terms of their amplitude and phase, we see that real cross-spectra can only have phases of 0 or π and pure imaginary cross-spectra only have phases of ±(π/2). These phase results provide pow-
erful checks on the consistency of any gradient signal and are summarized in Table 1.
Isotropy An important special case of this spectral theory is when the 2D PSD is assumed isotropic [i.e., S 0 (kx , ky ) = S 0 (k), where k is the distance, in wavenumber space, from the origin]. When isotropy holds the 1D PSDs (equations 16–21) can, via the sub- stitution k^2 y = k^2 − k^2 x , be written
P [U xx ](k x ) = 32 π^4
k (^) x
k x^4 k(k^2 − k x^2 )−^
1
2 e−^4 πkh^ S 0 (k)dk,
P [U yy ](k x ) = 32 π^4
k (^) x
k(k^2 − k x^2 )
3
2 e−^4 πkh^ S 0 (k)dk, (31)
P [U zz](k x ) = 32 π^4
k (^) x
k^5 (k^2 − k x^2 )−^
1
2 e−^4 πkh^ S 0 (k)dk,
P [U xy ](k x ) = 32 π^4
k (^) x
k x^2 k(k^2 − k x^2 )
1
2 e−^4 πkh^ S 0 (k)dk,
Table 1. Properties of PSDs and cross-spectra of gravity gradients. In column 2, symbols that match identify which spectra are equated by the interrelations, equations 24–29.
Allowed General form Isotropic form Spectra Equalities^1 phases Equation^2 Equation Figure^3
P Uxx 0 16 30 P Uyy 0 17 31 P Uzz Two rules, 0 18 32 2, 4 equation 22 P Uxy 0 19 33 P Uxz ♦ 0 20 34 P Uyz 0 21 35 P×Uxx , Uyy 0 , π A-1 A-16 3a, 5a P×Uxx , Uzz ♦ 0 , π A-2 A-17 3b, 5b P×Uxx , Uxy 0 , π A-4 0
P×Uxx , Uxz
π 2
3 π 2
A-5 A-19 3d, 5d
P×Uxx , Uyz ♣
π 2
3 π 2
A-6 0
P×Uyy , Uzz 0 , π A-3 A-18 3c, 5c P×Uyy , Uxy 0 , π A-7 0
P×Uyy , Uxz
π 2
3 π 2
A-8 A-20 3e, 5e
P×Uyy , Uyz
π 2
3 π 2
A-9 0
P×Uzz , Uxy ♠ 0 , π A-10 0
P×Uzz , Uxz
π 2
3 π 2
A-11 A-21 3f, 5f
P×Uzz , Uyz
π 2
3 π 2
A-12 0
P×Uxy , Uxz ♣
π 2
3 π 2
A-13 0
P×Uxy , Uyz
π 2
3 π 2
A-14 A-22 3g, 5g
P×Uxz , Uyz ♠ 0 , π A-15 0
(^1) Spectra equated by interrelations. (^2) Zero indicates a zero spectrum. (^3) First number denotes coherence-free case, B = 0; second number refers to self-similar case, β = 3.
Spectral Analysis of Gravity Gradient J
P [U xz](k x ) = 32 π^4
k (^) x
k x^2 k^3 (k^2 − k x^2 )−^
1
2 e−^4 πkh^ S 0 (k)dk,
P [U yz](k x ) = 32 π^4
k (^) x
k^3 (k^2 − k x^2 )
(^12)
e−^4 πkh^ S 0 (k)dk.
Making the same substitution also allows the cross-spectra to be written in isotropic form (see Appendix A). Eight of the cross-spectra vanish identically, meaning that their com- ponents are completely uncorrelated with one another. As a result, only seven cross-spectra and six PSDs are sensitive to the form of S 0 (k). However, because of the two power-sum rules (equation 22) and the interrelations (equations 24–27), the number of independent variables reduces to seven, as in Table 1. Two interrelations, equations 28 and 29, are lost be- cause they are zero under the assumption of isotropy. By specifying the form of S 0 (k), we can derive the form of the PSDs and cross-spectra and calculate the coherence from equation 23. For example, Appendix B derives the form of S 0 (k) for a self-similar source distribution, where S 0 (k) is pro- portional to k−(β+3)^. The equations above and in Appendix A provide a com- plete description of the spectra of the gradient tensor, the properties of which are summarized in Table 1. We use these properties to diagnose signals of nonstationary or nongeo- physical origin. Such signals immediately reveal themselves by causing violations in the spectral equalities (equations 22, 24–
- and by causing low coherence and unrealistic phases of the cross-spectra. Additionally, if a form of S 0 (kx , ky ) is cho- sen, then the spectral equations can be used to invert for the average depth h and other relevant parameters, such as the scaling constant β for self-similar sources.
Two-component gradiometers
In addition to full tensor gravity gradiometers, a number of systems (e.g., Lee, 2001) use a single horizontal spinning disk and only measure the Uxy component and the inline compo- nent (Uxx −Uyy ). For completeness we state without proof that the 1D and 2D PSDs for the difference component are
− 2 SU xy , (36)
P [U xx − U yy ](k x ) = P [U xx ](k x ) + P [Uyy ](k x )
− 2 P [U xy ](k x ). (37)
By combining the inline PSD (equations 36 and 37) and the power-sum rules (equation 22), we derive the relationships
P [U zz](k x ) = P [U xx − U yy ](k x ) + 4 P [Uxy ](k x ). (39)
We give these results as a simple way of calculating the PSD of the Uzz gradient component from a two-component gradiome- ter. We note in passing that the power-sum rules (equation 14) and the inline PSD (equation 38) are all invariant under ro-
tations around the z -axis (this follows because Uzz must be invariant under such rotations). In fact, the inline PSD and the second power-sum rule bear remarkable similarity to the differential curvature amplitude R =
(Uxx − Uyy ) 2 + 4 U (^) xy^2
and the amplitude of the horizontal gradient H =
U (^) xz^2 + U (^) yz^2 (Dransfield, 1994). These invariants have been suggested as alternatives to using just the Cartesian components of the gra- dient tensor in gravity gradient inversion schemes (Condi and Talwani, 1999). By using Parseval’s theorem (Parker, 1994), we relate R and H to the PSD of Uzz:
V
S[Uzz](k x , k y )dk x dk y =
V
S[R](k x , k y )dk x dk y
V
S[H ](k x , k y )dk x dk y ,
where V denotes all of space.
SYNTHETIC EXAMPLE
To test our theory, we chose a suitable form for S 0 (kx , ky ) and constructed a numerical model of the potential. The model was then used to produce simulations of a single line gradient survey. A number of papers [e.g., Pilkington and To- doeschuck (1990); Maus and Dimri (1994)] propose modeling the subsurface density distribution as self-similar, with a 3D PSD (Z[−]) described by a power law
where β is a scaling parameter, κ is a constant, and ξ is the vector (kx , ky , kz). We show in Appendix B that, assuming a self-similar source distribution, S 0 (kx , ky ) ∝ κk−(3+β)^. This is the form used in our numerical computations. We present here the results from two numerical models: a model where β = 0, corresponding to a correlation-free source, and another model where β = 3, which is the best- fitting value for the data described in the real data section.
Correlation-free source model, β = 0
By setting β = 0, we follow the lead of Parker and O’Brien (1997) and model the source distribution as being a realiza- tion of a stationary stochastic process that is uncorrelated at all length scales. This choice provides a limiting case for the potential, the other limit being a layer of constant density (a Bouguer slab) which has the easily verified but uninteresting property that all of the gradient components, and thus their spectra, are zero. A correlation-free source is isotropic and is therefore described by the simplified set of spectral equa- tions 30–35 and A-16–A-22. For the computation we set κ = 1 and constructed a 2D grid in the spectral domain with random phase and amplitude given by k−^
(^32) except at k = 0, which was set to zero ampli- tude. We continued the grid upward to a prescribed height and then differentiated in the spectral domain to obtain grids for the six gradient tensor components. We then inverse Fourier transformed the six grids to obtain the gradients in the space domain, from which we extracted and used a single profile to
Spectral Analysis of Gravity Gradient J
tensor compared with the theoretical phases (listed in Table 1) and the theoretical coherences calculated using equation 23 by substituting our theoretical spectra into the coherence equa- tion 23. Most of the theoretical coherences increase with wavenumber, though at different rates. This is analogous to the effect seen by Parker and O’Brien (1997), where increas- ing coherence is predicted between magnetic vector measure- ments above a statistically stationary correlation-free source. However, the coherence of P×Uxx , Uyy (Figure 3a) only increases to 1/3 as opposed to 1. Also, the coherence of P×Uyy , Uzz (Figure 3c) decreases from 1 to 1/3, and P×Uyy , Uxz (Figure 3e) increases from zero to a maxi- mum before decreasing to 1/3 at the highest wavenumbers. Figure 3 shows that the best estimates of the phase and co- herence occur when the theoretical coherence values are high. This agrees with Priestley (1992) and shows how the variance of the cross-spectra, and thus the coherence, increases with decreasing coherence. So while our theory can predict the co- herence and phase, it also forces estimates of these quantities to be poor when the theoretical coherence is low. Nonetheless, Figure 3 displays good agreement with the theory.
Data approximation model, β = 3
Using the method described in the real data section, we cal- culated a best-fitting value of β = 3 for our North Sea data. We then used this value in a second numerical simulation. Figure 4 shows how well the power-sum rules (equation 22) are obeyed for this simulation. The rules hold and are in better agreement than for the β = 0 case (Figure 2), especially at low frequencies. The other main feature of this figure is the faster rate of decay of the signal with wavenumber. This is because k in spectral equations 30–35 is to a greater negative power. We also show in Figure 5 how the coherence and phase vary for the β = 3 case. Although there is visible low- frequency bias in the Figure 5 plots, the plots strongly agree with the theoretical predictions, especially when the coher- ence is high. More importantly, the figure shows that at high frequencies the coherences converge to the same values as the β = 0 case. However, at low frequencies the coherences are much flatter than before. By calculating the theoretical co- herence for a wide range of β values, we find that as β ap- proaches infinity, the coherences become constant over all wavenumber space, in particular: γ 2 = 1 /3 for P×[Uxx , Uyy ], P×[Uyy , Uzz], and P×[Uyy , Uxz] and γ 2 = 1 for P×[Uxx , Uzz], P×[Uxx , Uxz], P×[Uzz , Uxz], and P×[Uxx , Uyy ].
REAL DATA EXAMPLE
To see how the spectra of real data behave, we examined the spectra of a shipborne gradient survey conducted off the east coast of Scotland in the North Sea. A Bell Geospace full tensor gradiometer was used to measure the five indepen- dent components of the gradient tensor; Uzz was calculated from Uxx and Uyy using Laplace’s equation. Figure 6 shows the Uzz component measured by this survey, in both unfiltered (Figure 6a) and, preempting a later part of this paper, filtered (Figure 6b) forms. The survey was conducted over a 10- × 12-km approxi- mately rectangular area, generating some 150 lines of data with an along-line data spacing of approximately 40 m. A northeast-down coordinate system was used for the gradient
instrument, while the survey lines were oriented on a bearing of approximately 55◦. The water depth at lowest astronom- ical tide varied from 79 to 88 m across the survey site. The data had been previously line leveled and subjected to a noise- reduction algorithm. As mentioned in the preceding section, an inversion was conducted on the data to estimate β and h. To obtain these values, we substituted S 0 (kx , ky ) = κk−(3+β)^ into the PSD equa- tion for Uzz (equation 10) and then linearized the resulting equation by taking the logarithm of both sides. We then used a least-squares inversion to fit the radial spectrum of our data to this equation. The best-fitting values, calculated via this method, were found to be β = 3 .0 and h = 124 m. We see in Figure 7 that these values closely reproduce the radially av- eraged spectra. The fitted values are also reasonable under other considerations: h is greater than the water depth and
Figure 4. The power spectral density of the Uzz component of the gradient tensor (thin line), compared with power-sums (a) P Uxz + P Uyz (thick line) and (b) P Uxx + P Uyy + 2 P Uxy (thick line). All spectra are calcu- lated with 20 sinusoidal tapers and β = 3. The PSDs are of a simulated survey in the x-direction set at a height h one grid unit above the numerical model described in the text. As in Figure 2, the power sums have almost identical power, thus demonstrating the high level of agreement required by the power-sum rules (equation 22).
J18 While et al.
β is nearly 2.9, the predicted value for a Spector and Grant ensemble of blocks (Fedi et al., 1997). However, the data are anisotropic and are not well described by a self-similar approx- imation. As given, the formulas for the 1D spectra treat the x- axis as parallel to the survey direction. Therefore, to test the power-sum rules (equation 22) and the interrelations (equa- tions 24–29), we rotated the measurements into an along-track x, across-track y, and down z coordinate system. To perform the rotation, we found the line of best fit for individual data lines and used the bearing of the best-fitting line to rotate in- dividual data points. Cubic spline interpolation was then used to resample the data points to evenly spaced intervals along x, from which the spectra were calculated. Figure 8 shows the result of a test of the power-sum rules on a sample line (estimated with five sinusoidal tapers), taken from the North Sea data set. Though we only show one line,
Figure 5. The coherence (solid line) and phase (dotted line) of the seven self-similar (β = 3) nonzero isotropic cross-spectra. For comparison, the predicted values of the coherence (long dashes) and phase (short dashes) are given. The cross-spectra tend to be either real or pure imaginary, and the coherences agree with theoretical values calculated using equation 23. The higher value of β flattens out the plots, though all of the spec- tra still agree with the predicted values. All spectra were calcu- lated at a height h of one grid unit above the numerical model described in the text; low-frequency bias is visible in the spec- tra near kx h = 10 −^3.
the figure’s results are representative of the rest of the data. Both power-sum rules mostly hold within the plotted 95% confidence intervals and hold well out to a wavenumber of 0.001 m−^1. However, beyond this wavenumber the second power-sum rule (right-hand side of equation 22, Figure 8a) yields a result that is systematically higher than P Uzz,
Figure 6. The Uzz component of the North Sea gradient sur- vey: (a) unfiltered; (b) filtered with a second-order low-pass Butterworth filter with a cutoff wavelength of 0.001 m−^1. The scale is in Eov ¨os units (1 Eu = 10 −^9 s−^2 ). Plot (b) displays the features we consider geophysically reasonable.
Figure 7. A comparison between the true radially averaged spectrum (solid line) and our best-fitting, estimate (dashed line) β = 3. 0 , h = 124 m. The scale factor κ was poorly re- solved by our inversion; consequently, the vertical alignment of the spectra was determined by eye. The scale on the y -axis is in Eov ¨os units (1 Eu = 10 −^9 s−^2 ). The best-fitting spectrum is in good agreement with the radially averaged spectrum.
J20 While et al.
DISCUSSION AND CONCLUSIONS
We have given expressions for the 2D power spectral densi- ties of the components of the gravity gradient tensor and have shown that they are related to each other by two power-sum rules (equation 14). Under the assumption of stationarity, the power-sum rules may also be applied to 1D line data. Ad- ditionally, by giving the forms of the cross-spectra, we have shown that six interrelations involve the cross-spectra (equa- tions 24–29) and that the cross-spectra can be separated into two groups of equations: real and pure imaginary. These spec- tral relationships are similar in nature to relationships that ex- ist, in the spectral domain, for vector magnetic data. Because the spectral relationships are based on only a few assumptions — a flat earth approximation, noise-free data, and (for 1D data) stationarity — they are very powerful. For instance, in answer to the question posed in the introduction — how consistent are measured gradient components, consider- ing they are bound together by the form of the potential — the power-sum rules provide a method of checking the inter- nal consistency of the measured gradient tensor. Any signal for which the power-sum rules are broken is either nonsta- tionary (1D only) or nongeophysical and, as we did for the North Sea data, should be filtered out. If, in addition, isotropy of S[U (x, y)](kx , ky ) is assumed, eight of the cross-spectra, and thus their coherences, become zero; if any of these eight spec- tra are nonzero, then the assumption of isotropy is invalid and the spectra are anisotropic. Certain assumptions about the source distribution, such as the single self-similar scaling expo- nent, assume isotropy and are strictly invalid in the anisotropic case. The North Sea data in this paper were shown to be an ex- ample of an anisotropic source. If S 0 (kx , ky ) is well approximated, then estimates for the av- erage depth h can be made. Possible assumptions are self- similar sources (see Appendix B), correlation-free sources, or ensemble averages. From our North Sea data set we used a self-similar source assumption and calculated an average depth of 124 m, about 40 m below the sea bed. Care should be taken with this number, however, as the North Sea data were shown to be anisotropic. Our calculations show that a self-similar source distribution can be taken into account and that increasing the scaling expo- nent β affects the spectra in predictable ways: it steepens the decay of the PSDs and smooths the coherencies to constants.
ACKNOWLEDGMENTS
This work grew out of discussions with Robert Parker, whose advice and software are gratefully acknowledged. J. While was funded by a Nature Environment Research Council (NERC) studentship. The authors also thank the three reviewers for their helpful comments and suggestions.
APPENDIX A
CROSS-SPECTRA AND SPECIAL
ISOTROPIC FORMS
We present here the 15 1D cross-spectra of the gradient ten- sor. To obtain the equivalent 2D spectra, one can rewrite the
1D spectra without the integration:
P×[U xx , U yy ](k x ) =
−∞
16 π^4 k^2 x k^2 y
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-1)
P×[U xx , U zz](k x ) =
−∞
− 16 π^4 k^2 x | k |^2
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-2)
P×[U yy , U zz](k x ) =
−∞
− 16 π^4 k^2 y | k |^2
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-3)
P×[U xx , U xy ](k x ) =
−∞
16 π^4 k^3 x k y
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-4)
P×[U xx , U xz](k x ) =
−∞
i 16 π^4 k^3 x | k |
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-5)
P×[U xx , U yz](k x ) =
−∞
i 16 π^4 k^2 x k y | k |
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-6)
P×[U yy , U xy ](k x ) =
−∞
16 π^4 k^3 y k x
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-7)
P×[U yy , U xz](k x ) =
−∞
i 16 π^4 k^2 y k x | k |
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-8)
P×[U yy , U yz](k x ) =
−∞
i 16 π^4 k^3 y | k |e−^4 πkh^ S 0 (k x , k y )dk y ,
(A-9)
P×[U zz , U xy ](k x ) =
−∞
− 16 π^4 k y k x | k |^2
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-10)
P×[U zz , U xz](k x ) =
−∞
−i 16 π^4 | k |^3 k x
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-11)
P×[U zz , U yz](k x ) =
−∞
−i 16 π^4 | k |^3 k y
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-12)
P×[U xy , U xz](k x ) =
−∞
i 16 π^4 k^2 x k y | k |
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-13)
Spectral Analysis of Gravity Gradient J
P×[U xy , U yz](k x ) =
−∞
i 16 π^4 k^2 y k x | k |
× e−^4 πkh^ S 0 (k x , k y )dk y , (A-14)
P×[U xz , U yz](k x ) =
∞
16 π^4 k y k x | k |^2
× e−^4 πkh^ S 0 (k x , k y )dk y. (A-15)
By assuming isotropy, the cross-spectra become
P×[U xx , U yy ](k x ) = 32 π^4
k (^) x
k x^2 k
k^2 − k^2 x
2
× e−^4 πkh^ S 0 (k)dk, (A-16)
P×[U xx , U zz](k x ) = 32 π^4
k (^) x
−k^2 x k^3
k^2 − k^2 x
× e−^4 πkh^ S 0 (k)dk, (A-17)
P×[U yy , U zz](k x ) = 32 π^4
k (^) x
−k^3
k^2 − k^2 x
2
× e−^4 πkh^ S 0 (k)dk, (A-18)
P×[U xx , U xz](k x ) = 32 π^4
k (^) x
ik^3 x k^2
k^2 − k^2 x
× e−^4 πkh^ S 0 (k)dk, (A-19)
P×[U yy , U xz](k x ) = 32 π^4
k (^) x
ik x k^2
k^2 − k^2 x
2
× e−^4 πkh^ S 0 (k)dk, (A-20)
P×[U zz , U xz](k x ) = 32 π^4
k (^) x
−ik^4 k x
k^2 − k x^2
× e−^4 πkh^ S 0 (k)dk, (A-21)
P×[U xy , U yz](k x ) = 32 π^4
k (^) x
ik x k^2
k^2 − k^2 x
2
× e−^4 πkh^ S 0 (k)dk. (A-22)
APPENDIX B
POWER SPECTRAL DENSITY OF THE POTENTIAL
FOR A SELF-SIMILAR ISOTROPIC SOURCE
The autocovariance function C( r i , r j ) = E[U ( r i ), U ( r j )], where the potential U is a function of position vector r ) of a statistically stationary magnetic potential, is given in Jack- son (1990). Jackson’s result when rewritten for a statistically stationary gravitational potential is
C( r i , r j ) = G^2
V
V
| s ′^ − r i |
E[ρ( s ′), ρ( s )]
×
| s − r j |
d^3 sd^3 s′, (B-1)
where G is the gravitational constant, r i and r j are the position vectors of two observation points with components (xi , yi , 0) and (xj , yj , 0), ρ is the density, s and s ′^ are position vectors within the source volume V , with components (s 1 , s 2 , s 3 ) and (s′ 1 , s 2 ′, s′ 3 ), and E denotes expectation. Under the assumptions of isotropy and stationarity, the au- tocovariance of the density is purely a function of the separa- tion between the two points, i.e., E[ρ( s ′), ρ( s )] = a(| s ′^ − s |), where a(x) is an arbitrary function of x. Making this sub- stitution and reordering the integrals of the autocovariance, (equation B-1), we obtain
C( r i , r j ) = G^2
v′
| s ′^ − r i |
0
−∞
−∞
a(| s ′^ − s |)
×
(s 1 − xj ) 2 + (s 2 − yj ) 2 + s^23
ds 1 ds 2 ds 3 d^3 s′.
(B-2)
Using Parseval’s theorem (Parker, 1994) and the shift theorem of Fourier transforms, we convert the ds 1 and ds 2 integrals into integrals involving the wavenumbers ks 1 and ks 2 :
C( r i , r j ) = G^2
v′
| s ′^ − r i |
0
−∞
−∞
e−^2 iπ(k^ s^1 s
′ 1 +k^ s 2 s
′ 2 )
×A(k, |s 3 − s′ 3 |)e^2 iπ(k^ s^1 xj^ +k^ s^2 yj^ )
× B(k, s 3 )∗dk s 1 dk s 2 ds 3 d^3 s′, (B-3)
where A(k, |s 3 − s′ 3 |) is the 2D Fourier transform of a(s 1 , s 2 , |s 3 − s 3 ′|), B(k, s 3 ) is the 2D Fourier transform of 1 /
s 12 + s^22 + s 32 , k =
k^2 s 1 + k^2 s 2 , and the asterisk denotes complex conjugation. We now reorder the integrals to obtain
C( r i , r j ) = G^2
−∞
−∞
0
0
× A(k, |s 3 − s 3 ′|)B(k, s 3 )∗e^2 iπ(k^ s^1 xj^ +k^ s^2 yj^ )
×
−∞
−∞
| s ′^ − r i |
e−^2 iπ(k^ s^1 s
′ 1 +k^ s 2 s
′
2 )ds 1 ′ds′ 2
× ds′ 3 ds 3 dk s 1 dk s 2. (B-4)
From inspection, the bracketed integrals in equation B-4 form the 2D Fourier transform of 1/| s ′^ − r i |, which we write as B(k, s 3 ′ ). Using the shift theorem and making the substitutions x = xi − xj and y = yi − yj , we find
C( r i , r j ) = G^2
−∞
−∞
0
0
A(k, |s 3 − s 3 ′|)B(k, s 3 )∗
× B(k, s 3 ′ )ds′ 3 ds 3 e−^2 iπ(k^ s^1 x+k^ s^2 y)dk s 1 dk s 2.
(B-5)
Making the substitutions −kx = ks 1 and −ky = ks 2 , we obtain
C( r i , r j ) =
−∞
−∞
G^2
0
0
A(k, |s 3 − s 3 ′|)B(k, s 3 )∗
× B(k, s′ 3 )ds 3 ′ds 3
e^2 iπ(k^ x^ x+k^ y^ y)dk x dk y.
(B-6)
Written in the form of equation B-6, the autocovariance is a 2D inverse Fourier transform of the bracketed integrals. By definition the 2D power spectral density, if it exists, is the
