A reversible transform for seismic data processing

Theory

Many conventional seismic imaging steps are special cases of nonstationary shifts, which can be performed simultaneously with the Fourier transform via nonstationary convolution or combination (Margrave, 1998). This section shows how the nonstationary shift operator can be combined with the Fourier transform kernel to create a general data processing transform kernel.

It is important to define a convention for the forward and inverse Fourier transform before proceeding. The shift-direction determined by the sign in a shifting exponential will reverse if the opposite sign convention for the Fourier transform is used. We use the forward Fourier transform convention,

$$F \left( \beta \right) = \int_{-\infty}^{\infty} \, f \left( p \right) \, e^{-i \, \beta \, p} \, dp,$$ (1)

where $F(\beta )$ is the complex Fourier transform (commonly a function of frequency) of the input signal, $f(p)$ (commonly a function of time). With the Fourier transform convention defined in equation 1, a positive sign in the shifting exponential causes a backward shift of the input function. The inverse Fourier transform,

$$g \left( q \right) = \frac{1}{2\pi }\int_{-\infty}^{\infty} \, F \left( \beta \right) \, e^{i \, \beta \, q } \, d\beta ,$$ (2)

has an output function, $g(q)$ , that is strictly different from the original input function. A special property of the Fourier transform is that the input and output axes, $p$ and $q$ , are identical, and $f(p)=g(q)$ , meaning that the input function passes unchanged when put through both the forward and inverse Fourier transforms alone.

In general, any processing step that requires a mapping from an input coordinate that is a function of the output coordinate, to that output coordinate, can be viewed as a nonstationary shift by combination. Examples of such processing steps include normal moveout (NMO), dip-moveout (DMO), and frequency-wavenumber (f-k) migration. NMO is a nonstationary shift in the time-space domain, while f-k migration performs a similar shift in the frequency-wavenumber domain. To remain general though, it is better at this point to think in terms of input and output coordinates rather than input and output times, so that frequency, wavenumber, space, or time shifts can be admitted, depending on the desired domain.

To emphasize that the input and output coordinates do not necessarily correspond to time, we use the above $p-q$ notation. A general nonstationary data processing shift in this notation is,

$$\Delta \left( q \right)=p\left( q \right)-q,$$ (3)

where $p$ and $q$ are the input and output coordinates, respectively. Using these general coordinates, the nonstationary data processing shift applied to an input function, $f(p)$ , is,

$$g \left( q \right) \equiv f \left( q + \Delta \left( q \right) \r... ...\, e^{i \, \beta \, \Delta \left( q \right)} \, e^{i \, \beta \, q} \, d \beta,$$ (4)

where $F(\beta )$ is the Fourier transform of $f(p)$ , $\beta$ is the Fourier dual of $p$ , and the output signal, $g(q)$ is the desired result of processing the input signal $f(p)$ . Unlike the Fourier transform, $p$ and $q$ are not identical, although they represent the same dimension. The convenient form of $\Delta (q)$ allows the exponentials to be reduced by combining the Fourier transform kernel with the nonstationary shifting exponential. The result is a concise form of the forward data processing transform,

$$g \left( q \right) = \frac{1}{2 \, \pi} \int_{-\infty}^{\infty} F \left( \beta \right) \, e^{i \, \beta \, p \left( q \right)} \, d \beta,$$ (5)

which confirms the desired result,

$$g \left( q \right) = f \left( p \right).$$ (6)

The forward transform given in 5 is in the mixed-domain. For example, if the input function is defined over time, 5 takes the frequency spectrum of that input function, and outputs a new time series. This should not be alarming, as the Fourier transform itself is a mixed-domain transform. It may seem questionable how $p \left( q \right)$ is handled in the exponent. However, since the nonstationary shift is formulated here as a nonstationary combination, and the integration is over $\beta$ , nothing new is needed to handle the relation between $p$ and $q$ .

Although nonstationary combination kinematically performs nonstationary shifts correctly, amplitudes are not changed, and therefore energy is not conserved between input and output traces. Figure 1 demonstrates how the area under the input signal changes under a nonstationary combination shift. This nonphysical energy change within a trace is inherent to any processing step that applies or approximates a nonstationary shift. Rayleigh's theorem (also known as Plancherel's theorem) states that the integral of the square of an input function must be the same as the integral of its squared Fourier transform (Karl, 1989). Physically, this means that energy should be conserved between the input and Fourier domains, and clearly, the forward data processing transform violates this. Although nonphysical in most cases, nonstationary combination provides exactly the type of shift desired for the NMO correction and many other seismic data processing steps. The well-known issue of NMO stretch (Barnes, 1992) is a clear example of this type of energy change, where the seismic amplitudes remain unchanged after each trace has undergone the time-variant shift of the NMO correction. Relative to the input signal, this energy change means distortions in both the input and Fourier domains. These distortions are accounted for and handled by the inverse transform developed below, which relies on physically valid nonstationary convolution instead.

NSArea Figure 1. Element area change caused by a nonstationary combination shift. The values of $f$ at $p_1$ and $p_2$ are mapped to $g$ at $q_1$ and $q_2$ , respectively. Since $\Delta f$ does not equal $\Delta g$ , the areas $A_f$ and $A_g$ are not equal. Therefore, by Rayleigh's Theorem, energy is not conserved. Scaling $g \left ( q \right )$ by $\alpha (q)$ makes the areas the same.

To quantify how the energy has been changed by the forward transform, take the ratio of areas of a shifted and unshifted rectangular element as in Figure 1:

$$\alpha \left( q \right) = \frac{\Delta p \left( \frac{f \left( p_... ...\Delta q \left( \frac{g \left( q_1 \right) + g \left( q_2 \right)}{2} \right)}.$$ (7)

Equation 6 shows that the shifted values of the input function are not changed. This corresponds to the height of the rectangular elements shown in Figure 1 remaining unchanged by the forward transform ( $f(p_i)=g(q_i)$ ), giving,

$$\alpha \left( q \right) = \frac{\Delta p}{\Delta q}.$$ (8)

In the continuous limit, 8 goes to,

$$\alpha \left( q \right) =\frac{d p}{d q}.$$ (9)

The output function, $g \left ( q \right )$ , is of course distorted by different amounts for various $q$ since the applied shift is nonstationary. By scaling $g(q)$ by $\alpha (q)$ and integrating, we claim that the energy change can be accounted for via Rayleigh's theorem in the inverse transform.

Now that there is a mechanism for removing amplitude distortions, the rest of the inverse transform development comes from removing the kinematic shift. The output function of the forward transform, $g \left ( q \right )$ , becomes the input function of the inverse transform. Since the forward shift is applied simultaneously with the inverse Fourier transform, the reverse shift should be applied simultaneously with the forward Fourier transform. That is, the inverse data processing transform should naturally recover the original input spectrum, $F \left( \beta \right)$ , rather than go directly back to the original input signal.

We first apply the Fourier transform to $g \left ( q \right )$ in terms of its own coordinate, $q$ , then reverse the nonstationary shift simultaneously,

$$G\left(\gamma \right) = \int_{-\infty}^{\infty} \, \alpha \left( ... ...e^{-i \, \gamma \, q} \, e^{-i \, \gamma \, \Delta \left( q \right)} \, dq \, ,$$ (10)

where $\gamma$ is the Fourier dual of $q$ , and $G(\gamma )$ is the desired uncorrected spectrum. Two interesting simplifications come from 10. First, we can again exploit the convenient form of $\Delta (q)$ to reduce the exponential giving,

$$G \left( \gamma \right) = \int_{-\infty}^{\infty} \, \alpha \left( q \right) g \left( q \right) \, e^{-i \, \gamma \, p(q)} \, dq \, .$$ (11)

Doing so forces integrating $p \left( q \right)$ in the exponential over $q$ , which may at first seem like an unattractive approach, but we must do so to account for the nonphysical distortion caused by the forward transform. The second simplification helps, as the scaling function, $\alpha (q)=\frac{d p}{d q}$ conveniently acts as a change-of-integration-variable factor, giving,

$$G\left( \gamma \right) = \int_{-\infty}^{\infty} g \left( q \right)e^{-i \, \gamma \, p\left( q \right)}dp .$$ (12)

By Rayleigh's theorem, the energy change can be accounted for by scaling $g \left ( q \right )$ by $\alpha \left( q \right)$ inside an integral over $q$ . This observation, combined with 6, shows that 12 is equivalent to a regular Fourier transform of $f(p)$ . Therefore, we can recover the unprocessed Fourier domain function-the goal of the mixed-domain general inverse data processing transform,

$$G\left( \gamma \right) = F\left( \beta \right).$$ (13)

The original unprocessed input function can then be exactly recovered by a regular inverse Fourier transform:

$$h\left( q \right) = \frac{1}{2 \, \pi}\int_{-\infty}^{\infty}G \left( \gamma \right) e^{i \, \gamma \, q}d\gamma,$$ (14)

where, $h(q)$ is the final unprocessed function of the inverse Fourier transform. Equation 14 is just a regular inverse Fourier transform of $G(\gamma )$ over its own coordinate, $\gamma$ . Therefore, using 13, it is clear that 14 is equivalent to a regular IFT of $F(\beta )$ over $\beta$ , and that the desired result of the inverse transform is confirmed:

$$h\left( q \right)=f\left( q \right).$$ (15)

The kinematics of this general data processing transform are clear under nonstationary filtering theory. We claim that for any data processing step implemented by transform, the kinematics of the step can also be clearly understood, and that $\alpha \left( q \right)$ trivially predicts the amplitude corrections which are otherwise difficult to quantify using conventional filtering theory.

A reversible transform for seismic data processing