- Email: [email protected]
Corrections to: “A method for analyzing gravity anomalies due to a geologic contact by fourier transform”
Geoexploration, 19 (1981) 67-71 Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands
67
Short Communication CORRECTIONS TO: “A METHOD FOR ANALYZING GRAVITY ANOMALIES DUE TO A GEOLOGIC CONTACT BY FOURIER TRANSFORM”
S. CHACKO Union Oil Company of California, Los Angeles,
Calif. (U.S.A.)
(Received November 6, 1980; accepted December 10, 1980)
ABSTRACT Cbacko, S., 1981. Corrections to: “ A method for analyzing gravity anomalies due to a geologic contact by Fourier transform. Geoexploration, 19: 67-71.
In a recent paper (Chacko and Bhattacharya, 1980), we presented a method to analyze the measured first horizontal derivative of the gravity profile over a geologic contact. On p. 48 of that paper (hereinafter to be called ‘Paper I’), under the section “Comments on Errors”, our mathematical error analysis was not made for ‘white noise’ as claimed, but rather for some constant perturbation in the gravity measurements; such an error may be too selective and therefore not very meaningful. Certainly it would be a better approximation to consider the noise as ‘white’, and in this paper, a procedure is presented whereby the effects of random noise is reduced when using the method described in Paper 1. Assuming that the error e(3t) existing in our measurement g;(x) is totally ‘white noise’ (i.e., zero-mean, uncorrelated), we can write: &(x)
= g&r)
+ e(x)
(1)
where g, (X ) represents the signal. Taking Fourier transforms: F[gS, @)I
= JT&(x)l + JId~)l
Our interest, however, is in the second horizontal gravity gradient and instead of differentiating in the space domain, we can express equivalently in the frequency domain: F&&)3
=
wmi@)1
(2)
where k is the spatial frequency. Now let:
OOlS-7142/81/000~0000/$
02.50
0 1981
Elsevier Scientific Publishing Company
68
Gk. (LC)= autocorrelation = L&) * L&x) Taking Fourier
of gkx (x)
transforms:
= k2F[gx
(xc)* g, (-x) + e(x) * e(-x)]
so that: P’(k)
= k2P(k)
+ k2N(k)
(3)
where P(k) and N(k) are the power spectra of signal and noise respectively, is the power of the effective measurement, and we have assumed no correlation between noise and signal. At this point we will consider N(k) to be negligible compared with P(k); this is justified if the data has been smoothed or treated in any other way to improve the signal-to-noise ratio (stacking of parallel profiles for example), and is a common assumption in seismic data processing (Robinson, 1967). Observe that an improved amplitude spectrum may now be obtained:
P’(k)
A(k) = m) The procedure should be as follows: (1) multiply the Fourier representation of the measurement g:(x) by (‘jk) to obtain the transform of gxx(x); (2) compute the autocorrelation of this spectrum (its power); the square root of this zero-phase spectrum yields the required A(k). Note that A(k) is the same asin eqs. 10 and 11 of Paper 1: A(k)
= (27rGp sin d)eWk”
ln [A(k)]
= In (2nGp
sin d) - kz
(4) (5)
Assuming that the parameter d (slope of the fault plane) is known, this method of obtaining the parameters p and z should yield estimates that are an improvement over values obtained directly from the Fourier transform as suggested in Paper 1. Next we turn to the question of obtaining a good estimate of d from noisy data. As is the case in seismic processing, there is no simple way of obtaining improved signal information from the phase spectrum. To handle this problem the following least-squares error minimization procedure is suggested: Assume (g,,)j to be the effective measurement at a position Xj; each
69
such measurement, three unknowns: (g&j
+
ej
=
plus a correction,
as a function
of
(6)
forj = 1, . . . . VZ
fj(d,&Z)
The unknowns
may be represented
d,p,z must satisfy the requirement:
m
E = c
ey = minimum
j=l
The phase spectrum
(eq. 9, Paper 1):
Q(h) = -d rt ha
(7)
yields an approximate value for the slope d of the fault; using this we can obtam first order approximations for p and z by the procedure outlined above. The final unknowns may be written as: d=d’+Ad p=p’+Ap z=z’+Az
(9)
where (d,p,z) are the true values, (d’,p’,z’) are our approximations, (Ad,Ap,Az) are the corrections we need to determine. By Taylor’s theorem:
fj(d,p,z) = fj(dpk’) +
2
Ad +
‘$
&I
+
2
and
AZ
(9)
ignoring higher order terms. In order to determine (Ad,Ap,Az), we require to know fj(d’,p’,z’) and the derivatives af’/ad, 3 fj/ap, afj/az. From eq. 4 of Paper 1, we have: fj(d’,p’,Z’)
= [2Gp’ sin d’(z’
cos
d’ - Xj sin d’)]/(Jci2 + z’~),
j = 1, . . . . m (10)
so that fj(d’,p’,z’) is determined for the approximate values (d’,p’,z’). To obtain a fj/ad, 3 fj/ap, a fj/az, we use the following approximations: [fj(d’
+ b, P’J’) -
[fj(d’,p’ afj
ci=-
a2
1 C -28,
[fj(d’Tp:z’
fj(d’ - e,,p:z’)l
+ eb,Z’) - fj(d’,p’
+ b)
- fj(d’,p’,z’
- &&)I
- ec
)I
(11)
where the B’s are perturbations in the unknowns, and the six values of fj with perturbations are obtained using eq. (10) above.
70
Also, let: pj = C&~Ij - fjCd’,P’P’I SO
that Pi, 0’ = 1, . , . WI) tie known. Then:
-ej = ajAd + bjAp + CjAZ - Pj ZUSd: E=
g
+
j=l
2 j=l
fajAd + bjAp + CjAz -pi)”
We now minimize aE
E with respect
(14)
to each of the parameters
= 0 = 2 C i3Ad
(ajAd + bjAp + cjAZ - Pj>aj
a.E -
=0 = 2 C
(UjAd + bjAp f cjAX - Pj)bj
= 0 =
(ajAd + bjAP + cjA~ - Pj),cj
(15)
aAP
-
aAZ
2 g
Ad, Ap, AZ:
This yields the relation: -7
'CGjtJj
23 UjCj
r
zbjllj
CCjaj
Ad
EbjCj
CCjCj
A.2
-
-zUjPj+
_iEcjpj,
Note that the above 3 X 3 matrix is symmetric, and thus can be solved efficiently by recursive methods (Levinsan, 1947). Having obtained Ad, Ahp,Az in this rn~n~r, we can impruve the values of d ‘, p’,z’ and then repeat the minimization calculations with these values. The procedure is iterative, and at the end of each iteration we compute:
71
+ Dipi j=
9
(17)
1, . . . . m
so that E’ = E/m, where E’ represents the deviation still existing in the set (d’,p’,&). The iterative procedure is complete when E’ is less than some predefined tolerance. ACKNOWLEDGEMENT
I am I would helpful ing the
grateful to Union Oil Company for permission to publish this paper. also like to thank Ananthram Swami (Union Oil Company) for discussions on the treatment of noise in seismic data, and for reviewmanuscript.
REFERENCES Chacko, S. and Bhattacharya, B., 1980. A method for analyzing gravity anomalies due to a geologic contact by Fourier transform. Geoexploration, 18: 43-50. Levinson, N., 1947. The Wiener RMS (root-mean-square) error criterion in filter design and production. J. Math. Phys., 25: 261-278. Robinson, E.A., 1967. Predictive decomposition of time series with application to seismic exploration. Geophysics, 32: 418-484.
67
Short Communication CORRECTIONS TO: “A METHOD FOR ANALYZING GRAVITY ANOMALIES DUE TO A GEOLOGIC CONTACT BY FOURIER TRANSFORM”
S. CHACKO Union Oil Company of California, Los Angeles,
Calif. (U.S.A.)
(Received November 6, 1980; accepted December 10, 1980)
ABSTRACT Cbacko, S., 1981. Corrections to: “ A method for analyzing gravity anomalies due to a geologic contact by Fourier transform. Geoexploration, 19: 67-71.
In a recent paper (Chacko and Bhattacharya, 1980), we presented a method to analyze the measured first horizontal derivative of the gravity profile over a geologic contact. On p. 48 of that paper (hereinafter to be called ‘Paper I’), under the section “Comments on Errors”, our mathematical error analysis was not made for ‘white noise’ as claimed, but rather for some constant perturbation in the gravity measurements; such an error may be too selective and therefore not very meaningful. Certainly it would be a better approximation to consider the noise as ‘white’, and in this paper, a procedure is presented whereby the effects of random noise is reduced when using the method described in Paper 1. Assuming that the error e(3t) existing in our measurement g;(x) is totally ‘white noise’ (i.e., zero-mean, uncorrelated), we can write: &(x)
= g&r)
+ e(x)
(1)
where g, (X ) represents the signal. Taking Fourier transforms: F[gS, @)I
= JT&(x)l + JId~)l
Our interest, however, is in the second horizontal gravity gradient and instead of differentiating in the space domain, we can express equivalently in the frequency domain: F&&)3
=
wmi@)1
(2)
where k is the spatial frequency. Now let:
OOlS-7142/81/000~0000/$
02.50
0 1981
Elsevier Scientific Publishing Company
68
Gk. (LC)= autocorrelation = L&) * L&x) Taking Fourier
of gkx (x)
transforms:
= k2F[gx
(xc)* g, (-x) + e(x) * e(-x)]
so that: P’(k)
= k2P(k)
+ k2N(k)
(3)
where P(k) and N(k) are the power spectra of signal and noise respectively, is the power of the effective measurement, and we have assumed no correlation between noise and signal. At this point we will consider N(k) to be negligible compared with P(k); this is justified if the data has been smoothed or treated in any other way to improve the signal-to-noise ratio (stacking of parallel profiles for example), and is a common assumption in seismic data processing (Robinson, 1967). Observe that an improved amplitude spectrum may now be obtained:
P’(k)
A(k) = m) The procedure should be as follows: (1) multiply the Fourier representation of the measurement g:(x) by (‘jk) to obtain the transform of gxx(x); (2) compute the autocorrelation of this spectrum (its power); the square root of this zero-phase spectrum yields the required A(k). Note that A(k) is the same asin eqs. 10 and 11 of Paper 1: A(k)
= (27rGp sin d)eWk”
ln [A(k)]
= In (2nGp
sin d) - kz
(4) (5)
Assuming that the parameter d (slope of the fault plane) is known, this method of obtaining the parameters p and z should yield estimates that are an improvement over values obtained directly from the Fourier transform as suggested in Paper 1. Next we turn to the question of obtaining a good estimate of d from noisy data. As is the case in seismic processing, there is no simple way of obtaining improved signal information from the phase spectrum. To handle this problem the following least-squares error minimization procedure is suggested: Assume (g,,)j to be the effective measurement at a position Xj; each
69
such measurement, three unknowns: (g&j
+
ej
=
plus a correction,
as a function
of
(6)
forj = 1, . . . . VZ
fj(d,&Z)
The unknowns
may be represented
d,p,z must satisfy the requirement:
m
E = c
ey = minimum
j=l
The phase spectrum
(eq. 9, Paper 1):
Q(h) = -d rt ha
(7)
yields an approximate value for the slope d of the fault; using this we can obtam first order approximations for p and z by the procedure outlined above. The final unknowns may be written as: d=d’+Ad p=p’+Ap z=z’+Az
(9)
where (d,p,z) are the true values, (d’,p’,z’) are our approximations, (Ad,Ap,Az) are the corrections we need to determine. By Taylor’s theorem:
fj(d,p,z) = fj(dpk’) +
2
Ad +
‘$
&I
+
2
and
AZ
(9)
ignoring higher order terms. In order to determine (Ad,Ap,Az), we require to know fj(d’,p’,z’) and the derivatives af’/ad, 3 fj/ap, afj/az. From eq. 4 of Paper 1, we have: fj(d’,p’,Z’)
= [2Gp’ sin d’(z’
cos
d’ - Xj sin d’)]/(Jci2 + z’~),
j = 1, . . . . m (10)
so that fj(d’,p’,z’) is determined for the approximate values (d’,p’,z’). To obtain a fj/ad, 3 fj/ap, a fj/az, we use the following approximations: [fj(d’
+ b, P’J’) -
[fj(d’,p’ afj
ci=-
a2
1 C -28,
[fj(d’Tp:z’
fj(d’ - e,,p:z’)l
+ eb,Z’) - fj(d’,p’
+ b)
- fj(d’,p’,z’
- &&)I
- ec
)I
(11)
where the B’s are perturbations in the unknowns, and the six values of fj with perturbations are obtained using eq. (10) above.
70
Also, let: pj = C&~Ij - fjCd’,P’P’I SO
that Pi, 0’ = 1, . , . WI) tie known. Then:
-ej = ajAd + bjAp + CjAZ - Pj ZUSd: E=
g
+
j=l
2 j=l
fajAd + bjAp + CjAz -pi)”
We now minimize aE
E with respect
(14)
to each of the parameters
= 0 = 2 C i3Ad
(ajAd + bjAp + cjAZ - Pj>aj
a.E -
=0 = 2 C
(UjAd + bjAp f cjAX - Pj)bj
= 0 =
(ajAd + bjAP + cjA~ - Pj),cj
(15)
aAP
-
aAZ
2 g
Ad, Ap, AZ:
This yields the relation: -7
'CGjtJj
23 UjCj
r
zbjllj
CCjaj
Ad
EbjCj
CCjCj
A.2
-
-zUjPj+
_iEcjpj,
Note that the above 3 X 3 matrix is symmetric, and thus can be solved efficiently by recursive methods (Levinsan, 1947). Having obtained Ad, Ahp,Az in this rn~n~r, we can impruve the values of d ‘, p’,z’ and then repeat the minimization calculations with these values. The procedure is iterative, and at the end of each iteration we compute:
71
+ Dipi j=
9
(17)
1, . . . . m
so that E’ = E/m, where E’ represents the deviation still existing in the set (d’,p’,&). The iterative procedure is complete when E’ is less than some predefined tolerance. ACKNOWLEDGEMENT
I am I would helpful ing the
grateful to Union Oil Company for permission to publish this paper. also like to thank Ananthram Swami (Union Oil Company) for discussions on the treatment of noise in seismic data, and for reviewmanuscript.
REFERENCES Chacko, S. and Bhattacharya, B., 1980. A method for analyzing gravity anomalies due to a geologic contact by Fourier transform. Geoexploration, 18: 43-50. Levinson, N., 1947. The Wiener RMS (root-mean-square) error criterion in filter design and production. J. Math. Phys., 25: 261-278. Robinson, E.A., 1967. Predictive decomposition of time series with application to seismic exploration. Geophysics, 32: 418-484.