Theory Overview

Primary reflections in a common midpoint gather exhibit a hyperbolic moveout as a function of offset (gray lines on the left of Figure 3). The governing equation of the hyperbolic moveout is:  

$$t_x=\sqrt{t_0^2+\frac{x^2}{V_s^2}},$$ (1)

where t_x corresponds to the arrival time of the reflection at offset x, t₀ corresponds to the arrival time at zero offset and V_s is the NMO-stacking velocity. This velocity is the one that best fits the moveout of the hyperbola and is determined by trial and error from among a series of probable velocities. If correctly chosen, this velocity allows the moveout corrected primary reflections to become horizontal (solid black horizontal lines in left of Figure 3).

Clearly the selection of the stacking velocities must be done to correct for the moveout of the primary reflections and not for the multiples. At a given zero-offset arrival time the velocity of a primary reflection is greater than that of a multiple, which according to Equation (1) implies a smaller moveout. This difference in moveout makes it possible to flatten the primary reflections while leaving the multiples under-corrected with a moveout approximately parabolic Hampson (1986). The Parabolic Radon Transform exploits this difference by summing trace amplitudes along parabolas of different zero-offset time and curvature. Hence, the transform can be considered a mathematical operator that maps parabolas in the t-x domain to small regions of the parabolic moveout (p) and zero-offset time ($\tau$) domain.

 

mul_esq2
Figure 3 Schematic representation of multiple suppression by filtering in the parabolic Radon transform domain. Left: t-x domain, right: $\tau$-p domain.

This is schematically shown in Figure 3 which shows that the horizontal events in t-x domain map to a vertical strip in the $\tau$-p domain at p=0. The multiple reflections, on the other hand, are mapped in the $\tau$-p domain to a region away from the p=0 vertical line. This separation allows for the suppression of the multiple energy by zeroing out the $\tau$-p region to the right of the dashed line in Figure 3. The inverse $\tau$-p transform would then return the primaries to the t-x domain.

In practice the process is applied a little differently: it is the energy of the primaries that is suppressed (energy to the left of the dashed line in Figure 3) and inversely transformed to the t-x domain. The primaries are computed by subtracting the multiples from the original data in this domain. This method was first introduced with the name ``inverse velocity stacking'' Hampson (1986).

The mathematical equivalent of the qualitative description given before for the Parabolic Radon Transform is a set of two equations:  

$$y(p,\tau)=\int_{x_{min}}^{x_{max}} z(x,t=\tau+px^2)dx$$ (2)

 

$$z(x,t)=\rho (t)\ast \int_{p_{min}}^{p_{max}} y(p,\tau=t-px^2)dp$$ (3)

The first equation corresponds to the forward transform (from t-x to $\tau$-p) and the second one to the inverse transform (from $\tau$-p to t-x). z and y represent the trace amplitudes in t-x and $\tau$-p domain respectively. x_min and x_max correspond to the minimum and maximum CMP offset, p_min and p_max to the minimum and maximum parabola curvature used in the transform, and, as usual, the symbol $\ast$ denotes convolution. It is interesting to note that except for the difference in sign and the presence of the $\rho$term, the equations for both transforms are basically the same. The term $\rho$ represents a filter that corrects the high frequency loss incurred in the forward transform Claerbout (1995). In the case of continuous functions these transforms are exact inverses of one another. In seismic data processing we deal with sampled information, however, which means that we need to use the discrete equivalents of equations (2) and ( 3):  

$$y(p_i,\tau)=\sum_{k=0}^{N_x-1} z(x_k,t=\tau+p_ix_k^2)\Delta x$$ (4)

 

$$z(x_k,t)=\rho (t)\ast \sum_{t=0}^{N_p-1} y(p_i,\tau=t-p_ix_k^2)\Delta p$$ (5)

where N_x and N_p are the number of traces and parabolas respectively.

The need to work with discrete equations may give raise to aliasing problems Yilmaz (1987) as well as to some stability problems related to the selection of the number of parabolas used in the transform. In commercial software packages the transform is normally implemented in the f-x domain because of issues related to the amplitude of the inverse transform which is computed via a numerical optimization process. The discussion of these details, which are very important for the successful application of the method, are out of the scope of this paper. See for example Anderson (1993) and Alvarez (1995).