WO2000042448A1

WO2000042448A1 - Method of attenuating noise in three dimensional seismic data using a projection filter

Info

Publication number: WO2000042448A1
Application number: PCT/GB2000/000058
Authority: WO
Inventors: Ali ÖZBEK; Ahmet Kemal ÖZDEMIR; Ralf Ferber; Kamil Zerouk
Original assignee: Schlumberger Holdings Limited; Schlumberger Canada Limited; Services Petroliers Schlumberger
Priority date: 1999-01-14
Filing date: 2000-01-12
Publication date: 2000-07-20
Also published as: EP1145046B8; NO20013501L; EP1145046A1; US6668228B1; NO20013501D0; DE60023109D1; AU773131B2; EP1145046B1; NO332712B1; CA2358512C; AU1992300A; CA2358512A1

Abstract

A method of attenuating noise in three dimensional seismic data using a projection algorithm is disclosed. A frequency-space-space ('f-xy') projection algorithm is used which is a generalization of the f-x projection algorithm. The predictability of the seismic signals in the f-xy domain constitutes the basis of the algorithm. Specifically it is demonstrated that if the seismic events are planar in the t-xy domain, then in the f-xy domain they consist of predictable signals in the xy-plane for each frequency f. A crucial step of the 2-D spectral factorization is achieved through the helical coordinate transformation. In addition to the disclosed general algorithm for arbitrary coherent noise, a specialized algorithm for random noise is disclosed . It has been found that the disclosed projection algorithm is effective even in extreme cases of poor signal to noise ratio. The algorithm is also signal preserving when the predictability assumptions hold.

Description

Method of Attenuating Noise in Three Dimensional Seismic Data Using A

Projection Filter

Field of the Invention:

The present invention relates to the field of seismic data processing. In particular, the present invention relates to methods of attenuating noise in three- dimensional seismic data.

Background of the Invention:

A seismic signal which consists of only linear events has an f-x domain representation which is predictable in x for each frequency/ A generalization of this results is easy to show: the/- y domain representation of a 3-D (time-space-space) seismic signal, which consists of only planar events, is perfectly predictable in the xy- plane for each frequency/ The significance of this result is that, some problems such as reducing noise in a 3-D volume can be reduced to a set of 2-D problems in the xy- plane, which are easier to solve. For each frequency/ the noise is attenuated regardless of the data at other frequencies. For this purpose, what is needed is an algorithm which separates 2-D predictable data from additive noise. F-xy prediction (i.e./- y decoή) is conventionally used for this purpose. For example, see the following references: M. Chase, "Random noise reduction by FXY prediction filtering," EAGE Conf. Exp. Abs., Paris, pp. 164-165, 1992; M. Chase, "Random noise reduction by 3-D spatial prediction filtering," SEG Ann. Mtg. Exp. Abs., New

Orleans, pp. 1152-1153, 1992; and N. Gulunay, V. Sudhakar, C. Gerrard, and D.

Monk, "Prediction filtering for 3-K poststack data," SEG Ann, Mtg. Exp. Abs.,

Washington, ZλC, pp. 1183-1186, 1993.

Unfortunately, the/- y prediction methods suffer from model inconsistency problems. The model inconsistency in the/- y prediction method adversely affects signal preservation and noise attenuation when applied to seismic data. This disadvantage is similar to the model inconsistency problem in the 1-D counterpart to the/-xy prediction method, namely the/-* prediction algorithm. For an example of/-jc prediction, see: Canales, "Random noise reduction," 54^th SEG Ann. Mtg. Exp. Abs., Atlanta, pp. 525-527, 1984; and N. Gulunay, "FXDECON and complex wiener prediction filter," SEG Ann, Mtg. Exp. Abs., Houston, pp. 279-281, 1986.

The f-x projection algorithm is described in for example, the following references: R. Soubaras, "Signal-preserving random noise attenuation by the f-x projection," SEG Ann Mtg. Exp. Abs., Los Angeles, pp. 1576-1579, 1994; R Soubaras, "Deterministic and statistical projection filtering for signal-preserving noise attenuation," EAGE Conf. Exp. Abs., Glasgow, A051, 1995; R. Soubaras, "Prestack random and impulsive noise attenuation by f-x projection filtering," SEG Ann. Mtg. Exp. Abs., Houston, pp. 711-714, 1995; R. Soubaras, "the necessary and sufficient condition for lossless sampling," EAGE Conference and Technical Exhibition, Geneva, 1997; and U.S. Patent No. 5,771,203. However, the f-x projection algorithm has the limitation in that it is only applicable to 2-D data. Summary of the Invention:

Thus, it is an object of the invention to provide a noise attenuation algorithm that does not suffer from the same problems and limitations of the prior art. In particular, it is an object of the invention to provide an algorithm for attenuating noise in three-dimensional seismic data which does not suffer from model inconsistency problems known in conventional methods.

According to the invention, a method of attenuating noise in three dimensional seismic data is provided. The method includes receiving seismic data representing data gathered in at least two spatial dimensions and a time dimension. The seismic data includes both a noise component and a seismic signal component. The latter of which represents signals originating from at least one seismic disturbance. Values are computed for use in a projection filter, which is used to estimate the noise component of the seismic data. Spectral factorization is then performed in at least two dimensions to obtain additional two-dimensional values for use in the projection filter. The noise component in said received seismic data is estimated by using the projection filter which includes at least some of the additional two-dimensional values. The estimated noise component is then subtracted from the received data to obtain attenuated seismic data having a decreased noise component.

According to a preferred embodiment of the invention, the method also includes causing at least one seismic disturbance, recording raw data from a plurality of sensors distributed in at least two spatial dimensions; and then processing the recorded raw data to form said seismic data.

According to a preferred embodiment of the invention, the method also includes performing a Fourier Transform with respect to time of the seismic data to obtain frequency domain seismic data, selecting a single frequency from said frequency domain seismic data; and repeating for each desired frequency said steps of computing values, performing spectral factorization, estimating the noise component, and subtracting the estimated noise.

According to a preferred embodiment of the invention, the method also includes creating an initial estimate for an initial spectral factor sequence of values to be used in the projection filter, applying an all-pole filter based on the inverse of the square of the initial spectral factor sequence of values to obtain an intermediate sequence of values; computing an autocorrelation of the intermediate sequence of values; and finding coefficients for use in the projection filter by solving normal equations using the autocorrelation of the intermediate sequence of values.

According to a preferred embodiment of the invention, the spectral factorization is performed using a helical coordinate transform on the autocorrelation of the intermediate sequence of values to obtain a one-dimensional sequence, and the factorization is performed on the one-dimensional sequence to obtain a one- dimensional factor, which is mapped into two-dimensions using an inverse of the helical coordinate transform to obtain a two-dimensional factor which represents some of said additional two-dimensional values.

According to a preferred embodiment of the invention, the noise component of the seismic data is primarily random noise, and the projection filter estimates primarily random noise.

According to another preferred embodiment of the invention, the noise component of the seismic data is primarily coherent noise, and the projection filter estimates primarily coherent noise.

Brief Description of the Drawings:

Figure 1 illustrates the allowable support for a 2-D causal sequence;

Figures 2a and 2b illustrate the support of the linear predictor and prediction error filter, respectively, with parameters X and Y;

Figure 3 illustrates an example of a 2-D autocorrelation sequence on the left and its helical coordinate transform on the right;

Figure 4 illustrates an example of a 1-D sequence on the left and its inverse helical coordinate transform on the right; Figure 5 illustrates an example of a helical coordinate transform or the φ-

Radon projection of an autocorrelation sequence;

Figure 6 shows the space domain implementation of the f-xy projection filtering, according to a preferred embodiment of the invention;

Figure 7 shows the computation of the projection filter output, according to a preferred embodiment of the invention;

Figure 8 illustrates an example of one in-line of the noise-free seismic data cube on the left and the noisy seismic data cube on the right;

Figure 9 illustrates the wavenumber spectra of the noisy data on the left, and the corresponding noise estimation filter on the right at the frequency f = 30 Hz;

Figure 10 illustrates the outputs of the/- y prediction algorithm on the left and the/- y projection algorithm on the right;

Figure 11 illustrates the estimation errors of the/- y prediction algorithm on the left and the/- y projection algorithm on the right;

Figure 12 illustrates one in-line of a real data cube; Figure 13 illustrates the output of the/-xy prediction algorithm when applied to the real data set;

Figure 14 illustrates the output of the f-xy projection algorithm when applied to the real data set, according to a preferred embodiment of the invention;

Figure 15 illustrates the noise estimate of the/- y prediction algorithm, magnified by 2X;

Figure 16 illustrate the noise estimate of the/- y projection algorithm, magnified by 2X, according to a preferred embodiment of the invention;

Figure 17 is a flow chart illustrating a method of attenuating noise in three dimensional seismic data using a projection filter, according to a preferred embodiment of the invention;

Figure 18 is a flow chart illustrating a method of computing values for a noise estimation filter, according to a preferred embodiment of the invention;

Figure 19 is a flow chart illustrating a method of recording raw data from seismic disturbances, according to a preferred embodiment of the invention. Detailed Description of the Invention

1 Review of Some 2-D Signal Processing Concepts

In this section, we review some 2-D signal processing concepts which will be useful in our development of the f-xy projection algorithm. Most of these concepts are extensions of their 1-D counterparts, although the generalizations are not always unique or straightforward.

Causality In 1-D, a sequence h[x] is causal if h[x] = 0 for x < 0. The extension of this definition to 2-D is not unique. Here we follow a common convention by calling a 2-D sequence causal, if its support is constrained into the shaded region in Fig. 1.

Stability A 2-D filter is stable in the bounded-input, bounded-output sense if its impulse response is absolutely summable.

Minimum— phase We define a 2-D minimum-phase (min-^> ) filter to be a 2-D, causal, stable filter which has a causal, stable inverse.

In the next section we review the prediction error filters and the methods used in their computation.

1.1 2-D Prediction Error Filters (PEF)

Consider the following linear prediction problem:

h[^χι V\ = ~ _* aχγ{m, n]h[x - m, y - n] , (1)

where aχγ[m, n] denotes the coefficients of the linear predictor with support 5 y and parameters X, Y as illustrated Fig. 2. The corresponding PEF

Aχγ(z_x, Z_y) = l + ^aXY[^X,

> (2)

with the support shown in Fig. 2 is causal according to our definition of causality. This linear PEF is also called the forward linear prediction error filter to emphasize the fact that only the past samples (recall our causality definition) of the predicted sample are used in the prediction.

The filter coefficients, are found by solving the normal equations

r_h[k, l] = - ∑ a_XY[m, n}r_h[k - m, l - n) for (k, I) _ S_XY , (3)

(m,n)ζSχγ where r_/t[/c, -] is the autocorrelation function of h[x, y}. If the autocorrelation sequence is computed using the formula

r_h[k, l] = - ∑ h[x + k, y + l}h^*[x, y] , for (k, I) € S_XY , (4) ^χ,y then the resulting method for computation of the PEF coefficients is called the co- variance method or the autocorrelation method, depending the range of summation in (4). In the covariance method, the range of summation is chosen such that only the known (available) data sample are used in the computation. In the autocorrelation method, this range is extended by assuming that unknown (unavailable) data samples are zero. That is, the known segment of the input sequence h[x, y] is extrapolated with zeros.

A third way, to compute the PEF coefficients is the modified covariance method or the forward-backward method. In this method the filter coefficients are still given by the normal equations (3), however the autocorrelation sequence is computed as

r_h[k, l] = - ∑(h[x + k, y + l]h^* [x, y} + h[x, y]h^*[x - k, y - l}) (5) ^χ,y

Again, the range of summation is chosen such that only the known values of the data are used in the summation.

1.2 2-D Spectral Factorization

The 2-D spectral factorization problem can be phrased as follows: Given an autocorrelation sequence r_c[k, I] with a real and nonnegative Fourier transform R_c(k_x, k_y) > 0, find a causal and min-φ sequence c[x, y) such that:

r_c[k, l] = c[x, y] * *c^*[-x, -y] (6)

= ∑ c[x + k, y + l]c^*[x, y] . (7) ^χ,y

We denote the 2-D convolution operator herein as **. The same problem can be equivalently expressed in the Fourier Transform domain as

R_c(k_x, k_y) = \C(k_x, k_y)\² , (8)

or in the ^-transform domain,

R_c(z_x, z_y) = C(z_x, z_y)C^*(l/z_x ^*, l/z;) . (9)

The 2-D spectral factorization problem differs significantly from its 1-D counterpart. For example, in the 1-D case, it is always possible to find the causal, min-φ spectral factor due to the fundamental theorem of algebra. Briefly, the fundamental theorem of algebra states that a polynomial of order N can be always factored as the product of N first order polynomials over the field of complex numbers. Using this theorem, the 1-D spectral factorization problem is easy to solve: First find all the zeros (roots) of the autocorrelation function, then synthesize a sequence with the subset of these zeros which are smaller than 1 in magnitude (i.e., inside the unit circle). An implication of this spectral factorization algorithm is that, a finite extent autocorrelation sequence has a finite extent spectral factor.

The 2-D spectral factorization problem is complicated by the lack of a theorem similar to the "fundamental theorem of algebra" which applies to only 1-D polynomials. Another complication arises because a "finite" extent sequence may have an "infinite" extent spectral factor. To solve this difficult problem, some researchers have proposed to convert this 2-D problem into a 1-D problem, solve this 1-D simplified problem, then map the obtained solution back to 2-D.

According to a preferred embodiment of the invention, we utilize an approach known in the field of statistical signal processing to solve the 2-D spectral factorization problem - see, e.g., J. L. Roux and F. Dubus, "On the relationship between 1-D and 2-D spectral factorization - application to the analysis of the stability of 2-D transfer functions," Signal Processing, vol. 16, pp. 219-232, 1989, and R. R. Galikere, E. I. Plotkin, and M. N. S. Swamy, "Two-dimensional spectral factorization in the radon space," IEEE Stat. Signal and Array Processing 7^th SP Workshop, pp. 15-18, June 1994.

In this method, first the zero-padded autocorrelation sequence r_c[k, I] is mapped to a 1-D sequence using the helical coordinate transformation as shown in Fig. 3 . Basically, this transformation involves concatenating the rows of the zero-padded autocorrelation sequence to obtain a 1-D sequence. Alternatively, we may view the helical coordinate transform of a 2-D sequence as its -Radon projection shown in Fig. 5 , where the angle φ is determined by the amount of zero-padding. Once the helical coordinate transform is obtained, its 1-D spectral factor is computed (Fig. 4) using one of the 1-D spectral factorization algorithms, such as the Wilson's method - see, e.g., G. Wilson, "Factorization of the covariance generating function of a pure moving average process," SIAM J. Numer. Anal., vol. 6, no. 1, pp. 1-7, March 1969. Finally, this 1-D sequence is mapped back to a 2-D sequence by inverse helical coordinate transformation as shown in Fig. 4. The resulting sequence is the same as the true 2-D spectral factor if it is of finite extent. Otherwise, the obtained 2-D sequence is only an approximation to the true spectral factor. In this case the accuracy of the approximation can be improved by increasing the amount of zero padding used in helical coordinate transformation.

The term helical coordinate transformation was first introduced by J. Claerbout in "Multidimensional recursive filters via a helix," Geophysics, vol. 63, pp. 1532— 1541, Sep-Oct 1998, where a number of applications were shown. The use this approach for 2-D spectral factorization as shown in this invention is new in the fields of seismics and geophysics; its use for filtering of 3-D data through f-xy projection is new for any field.

2 Model of the Seismic Traces

Consider the output of an areal array of receivers, which can be represented as a 3-D signal u(t, x, y) in the time x space x space domain. This trace ensemble u(t, x, y) can be modeled as a clean seismic component u_σ(t, x, y) in additive noise e(t, x, y):

u(t, x, y) = u₀(t, x, y) + e(t, x, y) . (10)

Furthermore, it will be assumed that the noise-free trace ensemble u₀(t, x, y) can be represented as the superposition of p arbitrary wavelets V_j(t), that are delayed with arbitrary delays d_j(x, y) at each value of the space variables x and y: p u₀(t, x, y) = V_jit - d_jix^ y)) . (11) Here the delay functions, d₃ (x, y), determine the shape of the seismic events in a 3-D volume. If we take the Fourier transform of u₀(t, x, y) along the time dimension, we find the f-xy domain representation of the clean seismic signal as:

U₀(f, x, y) = ∑ V₃(f)e-'^2*fd>^ , (12)

where V₃(f) denotes the Fourier transform of the wavelet v₃(t). If we assume that events are locally planar, then the delay functions can be represented as linear functions of the space variables x and y:

d₃(x, y) = a₃ + b₃x + c₃y . (13)

Note that, if the events are not-planar, we process the events through a spatial window so that the seismic events in this window are linear. Substituting (13) into (12) we find that in the f-xy domain the seismic signal can be represented as the sum of p complex exponentials

Uo(f, x, v) = ∑ A₃(f)e-^fi*W^>x+c>ri , (14)

with complex amplitudes A₃

f⁾ = IΛ(/)l^'^(/) • (15)

Consequently, a 3-D seismic signal which consists of only planar events, is perfectly predictable in the xy-plane for each frequency /.

Thus, in the f-xy domain the additive noise model can be written as

U(f, x, y) = U₀(f x, y) -r E(f, x, y) , (16)

where U₀(f, x, y) is predictable in the xy-plane for each frequency /. As we show in the next section, this fact establishes the theoretical basis for the use of f-xy domain noise attenuation algorithms. 3 f-xy Projection Algorithm for Noise Attenuation

Separation of noise from the signal is an important problem in geophysics. Usually it is necessary to separate the desired seismic data, from the irrelevant noise. A major type of background noise found in seismic data is usually additive in nature. Attenuation of this background noise simplifies the interpretation of seismic data and enables the use of further signal processing algorithms.

Attenuation of the background noise can be achieved either in the f-xy domain, or in the t-xy domain. The f-xy domain algorithms have the advantage of separating a three dimensional problem into independent two dimensional problems. According to a preferrd embodiment, we shall be dealing with f-xy domain random noise attenuation. After describing the general framework in the next section, we shall present the f-xy projection algorithm for noise attenuation.

3.1 f-xy Domain Noise Attenuation Algorithms

As discussed above, the additive noise model and the planar seismic signals assumption lead to a predictable signal model in the f-xy domain. This observation motivates the use of the following algorithm for noise attenuation in the f-xy domain:

1. Find the f-xy domain representation of the data by computing Fourier transformation of the noisy data u(t, x, y) along the time dimension:

U(f, x, y) = J u(t, x, y)e-^^t dt . (17)

2. For each frequency /, attenuate the noise component in (16) using the data only at that frequency

U₀(f, x, y) = T[U(f, x, y)) , V/, (18) where T is an unspecified 2-D (i.e., space x space domain) noise attenuation algorithm and U₀(f, x, y) is the f-xy domain representation of the seismic signal estimate.

3. Find the t-xy domain representation of the signal estimate by computing the inverse Fourier transform of U₀(f, x, y) along the time dimension:

u₀(f, x, y) = j U₀(f, x, y ²^ d/ . (19)

The practical implication of the idea embodied in this algorithm is that, instead of solving the 3-D noise attenuation problem, it is sufficient to solve a set of 2-D problems in the xy-plane, which are easier to solve: for each frequency /, the noise is attenuated regardless of the data at other frequencies (second step of the algorithm) . Thus, in the rest of this description we concentrate on the attenuation of the noise at an arbitrary, but fixed frequency /. For notational clarity, at this arbitrary but fixed frequency /, we suppress the frequency dependence of the variables. For instance we denote the f-xy domain representations of the signals u(t, x, y), u₀(t, x, y) and e(/, x, y) as

u[x, y] = U(f, x, y) : Available data at frequency / (20) u₀[x, y] = U₀(f, x, y) : Noise-free seismic component at frequency / (21) e[x, y] — E(f, x, y) : Random noise component at frequency / . (22)

With this convention (16) simplifies as

u[x, y) = u₀[x, y) + e[x, y] , (23)

where u₀[x, y] is a perfectly predictable signal. Hence the second step of a general f-xy domain noise attenuation algorithm can be restated as the separation of the additive noise e[x, y] in (23) from the predictable seismic signal u₀[x, y}. Although the f-xy prediction algorithm is commonly used with this objective in mind, this algorithm suffers from the model inconsistency problem like its 1-D counterpart, the f-x prediction algorithm. In the next section, we derive the f-xy projection algorithm as an alternative and consistent method. According to a preffered embodiment, the predictability assumption is relaxed, so that the resulting algorithm is more applicable to real data.

3.2 The f-xy Projection Algorithm for Noise Attenuation

In this section, we derive the f-xy projection algorithm for noise attenuation. We first derive the general equations for an arbitrary noise model, and then specialize to the random noise case. To develop a consistent algorithm, in the next section we give the details of the assumed data model.

3.2J Data Model

We will assume that u[x, y] consists of a quasi-predictable (rather then perfectly predictable) seismic component u₀[x, y] in additive noise e[x, y]. The noise can in general be random or coherent.

u x, y] = u₀[x, y\ + e[x, y] (24)

By quasi predictability we mean that the seismic component, u₀[x, y], filtered with PEF a[x, y] produces an output, eg[x, y], that is small in magnitude:

a[x, y] * u₀[x, y] = e_xg_l[x, y} ^ 0 , (25) where β_\ is a small number. By modeling the noise in general as coherent, we assume that when the noise is filtered with a PEF b[x, y] it also produces an output that is small in magnitude:

b[x, y] * e[x, y] = e₂#₂[£, y] ~ 0 . (26)

Based on this quasi-predictable model of the seismic signals, in the next section we derive the f-xy projection algorithm for noise attenuation.

3.2.2 Derivation of the Algorithm

Assume that we know the PEF a[x, y] for u₀[x, y] and PEF b[x, y] for e[x, y}. Define the cost function V in the wavenumber domain:

V = \G_x{k_x, k_y)\² + \G₂{k_x, k_y)\² , (27)

where G_\{k_x, k_y) and G₂{ _x, k_y) are the Fourier transforms of the signals g [x, y] and g₂[x, y], respectively. Substituting the definitions of gι[x, y] and g_2.£, y] into (27) we obtain:

V = k_y)E(k_x, k_y)\² . (28)

In this last equation we observe that € and e determine the relative importance of the two terms in the cost function:

• If ei << e₂ then the first term dominates the cost function, therefore the minimization of the cost function should produce A(k_x, k_y)U₀(k_x, k_y) = 0. This produces a perfectly predictable seismic component that is closest to the data. • If 6χ >> β than the second term dominates the cost function. This reflects our belief that residual noise level (i.e., prediction error resulting from application of b[x, y) to e[x, y\) is low. In this case the solution for U₀(k_x, k_y) will be only quasi-predictable, but it will be closer to the observed data.

If we use the additive noise model then we can express the cost function in (28) only as a function of the noise component:

V = . (29)

Minimization of V in (29) over the noise component gives the estimate of E(k_x, k_y):

\A(k_x, k_y)\

E(k_x, k_y) U (K_x, Ky (30)

\A(k_x, k_y) \ⁱ + e \B(k_, k_y)\^< (31)

where N(k_x, k_y) is the noise estimation filter and e is eι/e₂. Note that, only the relative value of β and e₂ is important, because only their ratio appears in (31). N(k_x, k_y) is a projection filter because its spectrum mainly consists of l's and O's apart from the transition zones:

is

Finally the estimate of the seismic component is computed using the additive noise model:

Uo(k_x, k_y) — U (k_x, k_y) — i_ k_x, _y) . (J )

In the rest of this work, we shall assume that f-xy projection algorithm is used in random noise attenuation. In this case the noise component is unpredictable, hence B(k_x, k_y) = 1. Therefore the noise estimation filter is given as:

"<*^■■*-> ⁼ Ϊ ^ ? ^■ <³⁴> and the estimate of the random noise component is obtained as:

E(k_x, k_y) = N(k_x, k_y)U(k_x, k_y) (35)

3.2.3 Space Domain Implementation of f-xy Projection Algorithm

First we define a min-ςό filter C(k_x, k_y) such that:

\C(k_x, k_y)\² = \A(k_x, k_y)\² + e² , (37)

i.e, the autocorrelation of the filter c[x, y] is the prewhitened autocorrelation of the PEF a[x, y\. In -domain this relation can be expressed as:

C(z_x, z_y)C*(l/z_x ^*, l/_Z;) = A(z_x, z_y)A^*(l/z_x ^*, l/z;) + e² . (38)

As explained in above, computing such a C(z_x, z_y) from its autocorrelation is called spectral factorization. With this definition of C(z_x, z_y), the expression for the computation of the random noise component in (36) becomes

E,z z ) = ^A'W* Q (*^■ > Utz z ) (39) ^{h {z}" ^z'^] C-(l/z_{i t} l/z;) C(z„ z_v) ^{U (z}" ^z>> ^{• (39}>

in -domain. This expression suggests that, we can compute e[x, y] by passing the input through a min-φ filter A(z_x, z_y)/C(z_x, z_y) and computing an intermediate output. Then this intermediate output is filtered with A^*(l/z^*, l/z^*)/C^*(l/z*, 1/z*) "backwards" in space which gives the output e[x, y}. These filtering operations are shown in Fig. 6.

The min-φ property of C(z_x, z_y) ensures that first filter is stable when filtering is performed "forwards" in space and the second filter is stable when the filtering is performed "backwards" in space.

3.2.4 Computation of the PEF

So far we have assumed that the PEF were known beforehand. However, this is rarely the case. If the PEF is not known it should be calculated by minimizing the output energy of the projection filter.

In Fourier domain we have the following expression:

However, since \A(k_x, k_y)\² /\C(k_x, k_y)\² is a projection filter its square can be approximated by itself:

— |C-J/c_-, ky) (42) where the definition of g[x, y] is given in Fig. 6. Thus

^χ>y ^χ,y so instead of minimizing the energy of e[x, y] we may proceed by minimizing the energy of g[x, y] which is a simpler task.

For this simplified case an iterative algorithm can be given:

4 Synthetic Data Example

In this section, we show the results of a simulation to compare the performance of the f-xy projection algorithm with that of the f-xy prediction algorithm. The noisy input record is obtained by adding white noise to a clean, synthetic data record. The signal to noise ratio is -1.2 dB (peak-to-peak) and -12.9 dB (RMS). The 3-D data cube used in this simulation consists of 60 x 60 traces and 251 time samples. Fig. 8 shows one in-line of the data cube. The zero-dipping events have non-zero dip in the cross-line direction. The wavenumber spectra given in Fig. 9 is obtained with e = 0.3, and PEF order 5, using the noisy signal at / = 30 Hz. As it can be verified, the spectrum of the noise estimation filter is close to 0 at the wavenumbers where the signal component is present. The output signals and estimation errors of the two algorithms given in Fig. 10 and Fig. 11 justify the conclusion that the f-xy projection algorithm has a superior performance than f-xy prediction.

5 Real Data Example

In this section, we compare the results of the application of the f-xy projection algorithm with that of the f-xy prediction algorithm on a real data set. Fig. 12 shows one in-line of the real data cube. The 3-D data set consists of 200 x 21 traces and Algorithm 1 Computation of the PEF.

Inputs:

The samples of the noisy data u[x, y] in a rectangular grid 0 < x < N_x - 1 and 0 < y < N_y - 1.

The initial estimate for c[x,y], denoted as c₀[x,y}. Initialization:

Set the iteration number q = 0. Steps of the Algorithm: repeat

Set q = q + 1.

Compute the output of the stable filter l/C_q(z_x,z_y) (see Fig.7):

Find PEF a_q[x,y] for η[x,y] using one of the methods described above in the section on 2-D prediction error filters.

Find the min-φ filter c_q+ι[x,y] which has the following autocorrelation

Cq+ι[x,y] *c^*+ι[-x,-y] = a_q[x,y] * a^* _q[-x, -y] + e²δ[x,y] .

until {∑_xy |c₉₊₁[x,y] - ,y )/(∑_χy |_C(7+1[x,y]|²) < Threshold Output:

The filters a[x,y] = a_q[x,y] and c[x,y) = c_q[x,y\.

Comments: Since an initial estimate of c[x,y] may not be available, it may be set to c[x, y] = δ[x, y] for the frequency / = 0. For the other frequencies, the initial estimate for c[x, y] may be chosen as the converged c[x] filter of the previous frequency. 1001 time samples. Only the first 100 traces of an in-line are shown. Fig. 13 and Fig. 14 depict the outputs of the f-xy prediction and the f-xy projection algorithms, respectively, f-xy projection has produced a cleaner image, and some weak events like the one at 0.35 seconds is better defined. More revealing are the difference plots Fig. 15 and Fig. 16 , which show the noise estimates by f-xy prediction and f-xy projection respectively. These images have been plotted with a 2x magnification to show more detail. It is clear from these figures that f-xy projection has removed more random noise, whereas f-xy prediction has removed coherent seismic reflection energy.

According to the invention, an f-xy projection algorithm has been provided. The algorithm exploits the predictability of the seismic signals in the f-xy domain. With computer simulations conducted on the synthetic data, it has been shown that the algorithm can be used at very low SNRs and it outperforms the f-xy prediction algorithm. While the preferred embodiments of the invention have been described, the examples and the particular algorithm described are merely illustrative and are not intended to limit the present invention.

Claims

Claims What is claimed is:

1. A method of attenuating noise in three dimensional seismic data comprising the steps of: receiving seismic data representing data gathered in at least two spatial dimensions and a time dimension, said seismic data including a seismic signal component and a noise component, said seismic signal component representing signals originating from at least one seismic disturbance; computing values for use in a projection filter, said projection filter used to estimate the noise component of said seismic data; performing spectral factorization in at least two dimensions to obtain additional two-dimensional values for use in the projection filter; estimating the noise component in said received seismic data by using the projection filter which includes at least some of the additional two- dimensional values; and subtracting the estimated noise component from the received data to obtain attenuated seismic data having a decreased noise component.

2. The method of claim 1 further comprising the steps of: causing the at least one seismic disturbance; recording raw data from a plurality of sensors distributed in at least two spatial dimensions; and processing the recorded raw data to form said seismic data.

3. The method of claim 2 further comprising the steps of: performing a Fourier Transform with respect to time of said seismic data to obtain frequency domain seismic data; selecting a single frequency from said frequency domain seismic data; and repeating for each desired frequency said steps of computing values, performing spectral factorization, estimating the noise component, and subtracting the estimated noise.

4. The method of claim 3 further comprising the steps of: creating an initial estimate for an initial spectral factor sequence of values to be used in said projection filter; applying an all-pole filter based on the inverse of the square of the initial spectral factor sequence of values to obtain an intermediate sequence of values; computing an autocorrelation of the intermediate sequence of values; and finding coefficients for use in the projection filter by solving normal equations using the autocorrelation of the intermediate sequence of values.

5. The method of claim 4 wherein said step of performing spectral factorization comprises: using a helical coordinate transform on the autocorrelation of the intermediate sequence of values to obtain a one-dimensional sequence; performing spectral factorization on the one-dimensional sequence to obtain a one-dimensional factor; and mapping the one-dimensional factor into two-dimensions using an inverse of the helical coordinate transform to obtain a two-dimensional factor which represents some of said additional two-dimensional values.

6. The method of claim 5 wherein the noise component of said seismic data is primarily random noise, and wherein the projection filter estimates primarily random noise.

7. The method of claim 5 wherein the noise component of said seismic data is primarily coherent noise, and wherein the projection filter estimates primarily coherent noise.

8. A method of attenuating noise in three dimensional seismic data comprising the steps of: a. causing at least one seismic disturbance; b. recording raw data from a plurality of sensors distributed in at least two spatial dimensions; c. processing the recorded raw data to produce seismic data, said seismic data capable of being expressed in at least two spatial dimensions and a time dimension, said seismic data including a seismic signal component and a noise component, said seismic signal component representing signals originating from said at least one seismic disturbance; d. performing a Fourier Transform with respect to time of said seismic data to obtain frequency domain seismic data; e. selecting a single frequency from said frequency domain seismic data; f. creating an initial estimate for an initial spectral factor sequence of values to be used in a projection filter; g. applying an all-pole filter based on the inverse of the square of the initial spectral factor to obtain an intermediate sequence of values; h. computing an autocorrelation of the intermediate sequence of values; i. finding coefficients for use in said projection filter by solving normal equations using the autocorrelation of the intermediate sequence of values; j. obtaining a one-dimensional sequence a using a helical coordinate transform on the autocorrelation of the intermediate sequence of values; k. performing spectral factorization on the one-dimensional sequence to obtain a one-dimensional factor;

1. mapping the one-dimensional factor into two-dimensions using an inverse of the helical coordinate transform to obtain a two-dimensional factor which represents additional values for use in said projection filter; m. estimating the noise component in said received seismic data by using the projection filter which includes at least some of the additional values and at least some of the coefficients; n. subtracting the estimated noise component from the received data to obtain attenuated seismic data having a decreased noise component; and o. repeating steps d trough m for each desired frequency.