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:
The choice of sign in the exponent determines the direction of extrapolation. The corresponding 2-D equations are obtained by setting ky = 0.
For a constant value of kz, equation () is the analytic solution to the ordinary differential equation,
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 , kx, and ky.
For depth variable velocity v(z), kz is considered constant for small depth intervals () where the velocity is approximately constant. Therefore, the 2-D version of equation () can be written as
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 W1, W2, W3, where
Referring to Figure , upward continuation of the wavefield can be implemented in the following sequence of steps:
F* W1 F AP,F represents the Fourier transform of the horizontal space variable and F* represents the inverse Fourier transform of the horizontal space variable.
F* W2 F [F* W1 F AP + BP].
F* W3 F [[F* W2 F [F* W1 F AP + BP]] + CP].
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 , the adjoint datuming operator is
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).