Filter design

A particularly efficient implementation of offset continuation results from a log-stretch transform of the time coordinate Bolondi et al. (1982), followed by a Fourier transform of the stretched time axis. After these transforms, equation () from Chapter  takes the form  

$$h \, \left( {\partial^2 \tilde{P} \over \partial y^2} - {\... ... i\,\Omega \, {\partial \tilde{P} \over {\partial h}} = 0 \;,$$ (125)

where $\Omega$ is the corresponding frequency, h is the half-offset, y is the midpoint, and $\tilde{P} (y,h,\Omega)$ is the transformed data Fomel (1995b, 2000b). As in other F-X methods, equation () can be applied independently and in parallel on different frequency slices.

Analogously to the case of the plane-wave-destructor filters discussed in the previous section, we can construct an effective offset-continuation finite-difference filter by studying first the problem of wave extrapolation between neighboring offsets. In the frequency-wavenumber domain, the extrapolation operator is defined in accordance with equation (), as follows:  

$$\widehat{\widehat{P}}(h_2) = \widehat{\widehat{P}}(h_1)\, Z_{\lambda}(kh_2)/Z_{\lambda}(kh_1)\;,$$ (126)

where $\lambda = (1 + i \Omega)/2$, and $Z_\lambda$ is the special function defined in equation (). The wavenumber k corresponds to the midpoint y in the original data domain. In the high-frequency asymptotics, operator () takes the form  

$$\widehat{\widehat{P}}(h_2) = \widehat{\widehat{P}}(h_1)\, F(... ...ga\,\psi\left(2 k h_2/\Omega - 2 k h_1/\Omega\right)\right]}\;,$$ (127)

where functions F and $\psi$ are defined in equations () and ().

Returning to the original domain, I approximate the continuation operator with a finite-difference filter of the form  

$$\hat{P}_{h+1}(Z_y) = \hat{P}_{h} (Z_y) \frac{G_1(Z_y)}{G_2(Z_y)}\;,$$ (128)

which is somewhat analogous to (). The coefficients of the filters G₁(Z_y) and G₂(Z_y) are found by fitting the Taylor series coefficients of the filter response around the zero wavenumber. In the simplest case of 3-point filters, this procedure uses four Taylor series coefficients and leads to the following expressions:

$$\begin{eqnarray} G_1(Z_y) & = & 1 - \frac{1 - c_1(\Omega) h_2^2 + c_2(\Omega) h... ...) h_1^2 + c_2(\Omega) h_2^2}{12}\, \left(Z_y + Z_y^{-1}\right)\;,\end{eqnarray}$$ (129) (130)

where

$$c_1(\Omega) = \frac{3\,(\Omega^2 + 9 - 4 i\,\Omega)}{\Omega^2\,(3+i\,\Omega)}$$

and

$$c_2(\Omega) = \frac{3\,(\Omega^2 - 27 - 8 i\,\Omega)}{\Omega^2\,(3+i\,\Omega)}\;.$$

Figure  compares the phase characteristic of the finite-difference extrapolators () with the phase characteristics of the exact operator () and the asymptotic operator (). The match between different phases is poor for very low frequencies (left plot in Figure ) but sufficiently accurate for frequencies in the typical bandwidth of seismic data (right plot in Figure ).

Figure  compares impulse responses of the inverse DMO operator constructed by the asymptotic $\Omega-k$ operator with those constructed by finite-difference offset continuation. Neglecting subtle phase inaccuracies at large dips, the two images look similar, which indicates the high accuracy of the proposed finite-difference scheme.

When applied on the offset-midpoint plane of an individual frequency slice, the one-dimensional implicit filter () transforms to a two-dimensional explicit filter with the 2-D Z-transform  

$$G(Z_y,Z_h) = G_1(Z_y) - Z_h G_2(Z_y)\;,$$ (131)

analogous to filter () for the case of local plane-wave destruction. Convolution with filter () is the regularization operator that I propose for regularizing prestack seismic data.

 

arg
Figure 30 Phase of the implicit offset-continuation operators in comparison with the exact solution. The offset increment is assumed to be equal to the midpoint spacing. The left plot corresponds to $\Omega=1$, the right plot to $\Omega=10$.

 

off-imp
Figure 31 Inverse DMO impulse responses computed by the Fourier method (left) and by finite-difference offset continuation (right). The offset is 1 km.