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:

(13) |

For a constant value of *k*_{z},
equation () is the analytic solution to the ordinary
differential equation,

For depth variable velocity *v*(*z*), *k*_{z} is considered
constant for small depth intervals ()
where the velocity is approximately constant.
Therefore, the 2-D version of equation () can be written as

(14) |

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

The phase-shift extrapolation operation corresponding to equation () is represented by the matrices(15) |

Referring to Figure , upward continuation of the wavefield can be implemented in the following sequence of steps:

- 1.
- The traces from level
*A*are upward continued to level*B*by*F*^{*}*W*_{1}*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*_{2}*F*[*F*^{*}*W*_{1}*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*_{3}*F*[[*F*^{*}*W*_{2}*F*[*F*^{*}*W*_{1}*F AP*+*BP*]] +*CP*].

(16) |

gazwed
Schematic representation of phase-shift upward continuation corresponding
to implementation of equation (). The data are
extrapolated from level A by matrix Figure 5 W, from level B
by matrix _{1}W, and from level C by matrix _{2}W.
_{3} |

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 , the adjoint datuming operator is

(17) |

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 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 interval. This is the adjoint of the double square root (DSR) operator (Yilmaz and Claerbout, 1980).

2/12/2001