Methodology

Canales' 1984 implementation is in the f-x domain, so he makes the implicit assumption that the data is time-stationary, i.e., that each trace is the convolution of a single time-invariant wavelet with the earth's random reflectivity sequence. Computationally, this approach is very efficient. However, since ground roll is highly dispersive, and thus nonstationary, we instead choose a t-x method utilizing nonstationary PEF's Claerbout (1998a); Clapp et al. (1999); Crawley (1999).

Consider the recorded data to be the linear superposition of coherent signal plus coherent noise:  

$$\bf d = s + n.$$ (1)

Also assume that both the signal and noise are predictable, i.e., made up of of one or more local plane wave segments. The prediction error (residual) of the convolution of the signal and the noise with the corresponding nonstationary PEF's $\bf S$ and $\bf N$ is then uncorrelated. Writing these ideas as convolutional ``fitting goals'', we have

$$\begin{eqnarray} \bf r_n \; \equiv \; Nn &\approx& 0 \nonumber \\ \bf r_s \; \equiv \; \epsilon Ss &\approx& 0 \end{eqnarray}$$ (2)

where $\bf S$ and $\bf N$ represent t-x convolution with the nonstationary signal and noise PEF's, respectively, and the scaling factor $\epsilon$ balances the energies of the residuals. Equation (2) can be rewritten in an equivalent, but more familiar, notation as the minimization of a quadratic objective fuction:

$$Q({\bf n,s}) \equiv {\bf \Vert r_n + \epsilon r_s\Vert}^2 \;=\; \bf{\Vert Nn + \epsilon Ss \Vert}^2,$$

but we will hereafter use the notation of equation (2).

Rewriting the constraint equation (1), $\bf n = d - s$, allows us to eliminate $\bf n$ from equation (2), and suggests a regularized least squares optimization problem for the unknown signal $\bf s$

$$\begin{eqnarray} \bf Ns &\approx& \bf Nd \nonumber \\ \bf \epsilon Ss &\approx& 0, \end{eqnarray}$$ (3)

which is the approach used by Abma 1995. It can be easily shown (see Appendix A) that the predicted signal of equation (3) is the same as the optimal Wiener reconstruction Castleman (1996); Leon-Garcia (1994) for the special case of uncorrelated signal and noise.

We now consider some issues involving the calculation of the nonstationary signal and noise PEF's, $\bf S$ and $\bf N$.

• Signal PEF
  Equation (3) treats the signal PEF $\bf S$ as known, but this assumption appears circular since the unknown signal $\bf s$ is itself required in order to estimate $\bf S$. Spitz (1999) proposes an estimate of the signal PEF in terms of a PEF, $\bf D$, estimated from the data, and a PEF, $\bf N$, estimated from a model of the noise.  

  $$\bold S = \bold D \bold N^{-1}$$ (4)

  $\bold N^{-1}$ is nonstationary deconvolution, now better-defined in more than one dimension since the advent of the helix transform Claerbout (1998b). By conceptualizing equation (4) as a division of inverse spectra, Spitz' conjecture makes sense intuitively, that the unknown signal spectrum is the data spectrum scaled by the inverse of the noise spectrum. In the t-x domain, we have seen heuristic proof of equation (4), and Sergey Fomel reports that a general proof is trivial in the f-x domain.
• Noise PEF
  The assumption that the noise PEF $\bf N$ is known is also circular, since the noise is implicitly part of the unknown model. In many circumstances it is possible to explicitly model the noise well enough to estimate $\bf N$, like a convolutional model for multiples (i.e., the Delft SRME method of Verschuur et al. (1992)) Unforunately, variation in near surface earth properties makes it nearly impossible to explicitly model ground roll to cover all observed cases from around the world, so it necessary to obtain the noise model directly from the data.

  The noise model should be a first-order estimate of the noise - if it was perfect, the job would already be done. Ideally the noise model should contain the spatial correlation of all noise events, which will realistically differ in amplitude from the actual noise by a small, arbitrary time- and space-variant scale factor. Of course, a small amount of signal will contaminate the noise model, but as the ground roll will have much higher amplitude than any embedded signal, we assume that the L₂ PEF estimation of $\bf N$ will effectively ignore the embedded signal, though this is not always the case.

  If the ground roll and signal are sufficiently well-separated in temporal frequency, application of a lowpass filter to the data may produce a satisfactory noise model. The degree to which the noise and signal are separated is dependent on many variables - geography, depth to target, geology, and source type, to name a few. Similarly, if a hyperbolic Radon transform focuses primaries well, then it may be possible to mute the primaries in $\tau-p$ space, then inverse transform to obtain a noise model. We have tested both approaches, and find that lowpass filtering is the most reliable, although this assertion is data-dependant. Clapp and Brown (1999) applied an L₁ iterative hyperbolic Radon transform to a multiple-infested 2-D seismic line, and found that the primaries were focused far better than for the usual L₂ least squares DRT.

Using Spitz' 1999 choice for the signal PEF, ${\bf S = DN}^{-1}$, rewrite equation (3)

$$\begin{eqnarray} \bf Ns &\approx& \bf Nd \nonumber \\ \epsilon \bold D \bold N^{-1} \bf s &\approx& 0. \end{eqnarray}$$ (5)

Following Fomel et al. 1997, we precondition this iterative problem to improve performance. Make the change of variables

$$\bold s = (\bold D \bold N^{-1})^{-1} \bold p = \bold N \bold D^{-1} \bold p,$$ (6)

and rewrite equation (5):

$$\begin{eqnarray} {\bf N N D}^{-1} \bold p &\approx& \bf Nd \nonumber \\ \bf \epsilon p &\approx& 0. \end{eqnarray}$$ (7)

This unconstrained least squares optimization can be solved using any gradient-based iterative technique, as in Claerbout (1998a).