Phase-shift datuming

The 3-D Kirchhoff integral (equation ()) can be written as a convolution which, as shown by Schneider (1978), can be Fourier transformed over x, y, and t coordinates so that it becomes complex multiplication in the frequency-wavenumber domain:  

$$P(k_x,k_y,z,\omega)\;=\; P(k_x,k_y,z_0,\omega) \; \cdot \; e^{i k_z \Delta z}.$$ (13)

The term k_z embodies the dispersion relationship, and is given by

$$k_z\;=\; \sqrt{\frac{\omega^2}{v^2} - k_x^2 - k_y^2}.$$

The choice of sign in the exponent determines the direction of extrapolation. The corresponding 2-D equations are obtained by setting k_y = 0.

For a constant value of k_z, equation () is the analytic solution to the ordinary differential equation,

$$\frac{\partial^2 P(k_x,k_y,z,\omega)}{\partial (\Delta z)^2} + k_z^2 P(k_x,k_y,z,\omega) \; = \; 0,$$

which forms the basis of Gazdag's (1978) phase-shift wave-extrapolation formulation. This form of the wave equation is valid for all values of $\omega$, k_x, and k_y.

For depth variable velocity v(z), k_z is considered constant for small depth intervals ($\Delta z$) where the velocity is approximately constant. Therefore, the 2-D version of equation () can be written as  

$$P(k_x,z=z_0 - \Delta z,\omega)={{P(k_x,z=z_0,\omega)}e^{ik_z \Delta z}}.$$ (14)

This form of the equation is used to downward or upward extrapolate the wavefield for small depth intervals.

The phase-shift datuming operation can be written in matrix form by concatenating a sequence of successive operator matrices. For simplicity, I consider upward continuation from a surface with three depth levels, to a flat output datum (Figure ). The algorithm can be generalized for any number of depth levels and any datum geometry.

The operation of phase-shift datuming begins by picking data that corresponds to depth levels on a numerical grid. This picking can be represented by multiplication with matrices that correspond to the shape of the datum. For the case shown in Figure , these matrices are

$$A= \left[ \begin{array} {cccccc} 1& & & & & \\ &1& & & & \\... ... & & &.& & \\ & & & &1& \\ & & & & &1 \\ \end{array}\right].$$

The phase-shift extrapolation operation corresponding to equation () is represented by the matrices W₁, W₂, W₃, where  

$$W_i= \left[ \begin{array} {cccccc} e^{ik_{zi1} \Delta z} &... ...& & \\ & & & & &e^{ik_{zi6} \Delta z} \\ \end{array}\right] .$$ (15)

Each diagonal contains the phase-shifting exponentials necessary to upward extrapolate the wavefield to a depth level. The value of (k_z)_ij is given by the dispersion relation,

$${(k_z)}_{ij}={\sqrt{{{\omega^2} \over v(z_i)^2}-k^2_{xj}}}.$$

Referring to Figure , upward continuation of the wavefield $P = P(x,z_s,\omega)$can be implemented in the following sequence of steps:

1.

The traces from level A are upward continued to level B by

F^* W₁ F AP,

F represents the Fourier transform of the horizontal space variable and F^* represents the inverse Fourier transform of the horizontal space variable.

2.

The traces from level B, plus the traces from the previous step, are upward continued to level C by

F^* W₂ F [F^* W₁ F AP + BP].

3.

Finally, the traces from level C , plus the traces from the previous step, are upward continued to the flat output datum:

F^* W₃ F [[F^* W₂ F [F^* W₁ F AP + BP]] + CP].

So, for each value of $\omega$, the output wavefield $P(x,z_{\rm dat},\omega)$is given by  

$$\displaystyle{ \left[ \begin{array} {c} F^* \\ \end{array... ...] \left[ \begin{array} {c} P(x,z_s,\omega)\end{array}\right]. }$$ (16)

where $P(x,z_s,\omega)$ is a column vector representing the wavefield recorded on the topographic datum.

 

gazwed Figure 5 Schematic representation of phase-shift upward continuation corresponding to implementation of equation (). The data are extrapolated from level A by matrix W1, from level B by matrix W2, and from level C by matrix W3.

The adjoint algorithm can be obtained by reversing the order of matrix multiplication in equation () and transposing and conjugating each matrix. The matrices A, B, C, are diagonal matrices and therefore equal to their transpose. For each value of the frequency $\omega$, the adjoint datuming operator is  

$$\displaystyle{ \left[ \begin{array} {ccc} A&B&C \end{arra... ...}\right] \left[ \begin{array} {c} F \\ \end{array}\right]. }$$ (17)

The matrices W^*_i are the adjoints of the diagonal matrices W_i, which are explicitly represented in equation ().

The matrix formulation presented here is similar to that developed by Popovici (1992) and by Ji and Claerbout (1992), except that their forward operators take data from a flat surface and upward continued to an irregular datum. Their formulations are more closely aligned with Reshef's (1991) algorithm for migration from irregular surfaces. In fact, Popovici's and Ji and Claerbout's forward operators are necessarily the adjoint of Reshef's downward continuation operator.

Regardless of how these methods are formulated, the underlying principle is the same. A computational grid is laid out, and the extrapolation begins either above or below the topography, the data are inserted into the extrapolation as the wavefield passes the grid points which coincide with the topography.

Prestack datuming can be performed as described for the Kirchhoff method by extrapolating shot and receiver gathers separately. However, since the phase-shift operator is applied recursively for each $\Delta z$ interval, prestack datuming can also be performed by definition of shot and receiver wavenumbers and implementation of equation () so that shot gathers and receiver gathers are upward continued separately over each $\Delta z$ interval. This is the adjoint of the double square root (DSR) operator (Yilmaz and Claerbout, 1980).