Methodology

Consider local windows of the seismic image to be the simple superposition of signal and noise:  

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

The frequency domain representation of the Wiener optimal reconstruction filter for uncorrelated signal and noise is Castleman (1996); Leon-Garcia (1994):  

$$\bold H = \frac{\bf P_s}{\bf P_s + P_n}$$ (2)

where $\bf P_s$ and $\bf P_n$ are the power spectra of the unknown signal and noise, respectively. Multiplication of $\bf H$ with the data spectrum gives an optimal (in the least squares sense) estimate of the spectrum of the unknown signal.

Abma (1995) and Claerbout (1998) solved a constrained least squares problem to separate signal from spatially uncorrelated noise:

$$\begin{eqnarray} \bf Nn &\approx& 0 \nonumber \\ \bf \epsilon Ss &\approx& 0 \\ \mbox{subject to} &\leftrightarrow& \bf d = s+n \nonumber \end{eqnarray}$$ (3)

where the operators $\bf N$ and $\bf S$ represent t-x domain convolution with prediction-error filters (PEF's) which decorrelate the unknown noise $\bf n$ and signal $\bf s$, respectively, and the factor $\epsilon$ balances the energies of the residuals. Explicitly minimizing the quadratic objective function suggested by equation (3) leads to the following expression for the predicted signal:

$$\bf s = \bold N^T \bold N \left( \bold N^T \bold N + \epsilon^2 \bold S^T \bold S \right)^{-1} \ \bold d$$ (4)

Since the frequency response of the PEF approximates the inverse spectrum of the data used to estimate it, we see that Abma's approach is similar to Wiener reconstruction.

If the noise is assumed a priori to be spatially uncorrelated, as in Abma (1995), the noise decorrelator $\bf N$ is simply the identity. Gaussian noise is in the nullspace of the PEF estimation, so the signal decorrelator $\bf S$ can be estimated reliably from the data, i.e., $\bf S=D$, where $\bf D$ is a data decorrelating filter. Otherwise, if the noise is correlated spatially, an explicit noise model is required to estimate $\bf N$, and an approach like the one used by Spitz (1999) to estimate $\bf S$. Modifying equation (3) to reflect Spitz's choice of $\bf S = DN^{-1}$ and applying the constraint $\bf n = d - s$ gives

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

When solved iteratively, the problem can be preconditioned to improve convergence. Following Fomel et al. (1997), we can 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)

Brown et al. (1999) solved equation (7) iteratively to suppress ground roll with complicated moveout patterns, where $\bf S$ and $\bf N$ are nonstationary t-x-domain PEF's. Clapp and Brown (1999) did the same for multiple reflections.

Unfortunately, the estimation of nonstationary PEF's is computationally costly, and it is often difficult to ensure that the filters are minimum-phase, a necessary requirement for stable deconvolution, as in equation (7). For the application at hand, the final result is not the estimated signal and noise, but simply the noise-to-signal ratio. It follows that the separation need not be perfect - just good enough to distinguish between regions of the data with gross similarity to the facies template from the rest of the data. A properly stacked or migrated seismic image should have no ``crossing dips,'' and so can be conceptualized as a single-valued spatial function of local dip angle. Not surprisingly, we have found that simple three-point ``steering filters'' Clapp et al. (1997), work well for the noise and data decorrelating filters, $\bf N$ and $\bf D$, required to solve equation (7). The only thing needed to set up the steering filters is an estimate of the local dip field of the data and facies template, for which the automatic dip scanning technique of Claerbout (1992) produces satisfactory results.

Assuming that a given 2-D wavefield u(t, x) is planar with unknown local dip p, the operator  

$$\frac{\partial}{\partial x} + p\frac{\partial}{\partial t}$$ (8)

will extinguish it. If $\delta_x$ and $\delta_t$ are finite difference stencils for the continuous partial derivatives above, then equation (8) can be rewritten as a convolution, and hence cast as a univariate optimization for p:  

$$r = \left(\delta_x + p\delta_t\right)*u(t,x) \approx 0.$$ (9)

Differentiating the quadratic functional r^T r with respect to p gives an optimal estimate of the local dip:  

$$p = -\frac{\delta_x u \cdot \delta_t u}{\delta_t u \cdot \delta_t u}$$ (10)