Matrix formulation of adjoint Kirchhoff datuming

Dimitri Bevc

Author has no known email address

ABSTRACT

I present the matrix formulation of Kirchhoff wave equation datuming and show that upward and downward continuation are adjoint to each other. By expressing the datuming operator explicitly in matrix form, it is possible to examine the nature of the Hessian and to show that Kirchhoff datuming is not idempotent.

The Kirchhoff datuming codes I presented in earlier SEP reports (, ) are adjoint operators which pass the dot-product test; however, the adjoint formulation was arrived at by algorithmic considerations, not by formulating the operators in terms of matrix multiplications and taking their conjugate transpose. In this paper, I fill in this theoretical gap by presenting the matrix formulation of the datuming algorithm. I also show how Kirchhoff datuming can be cast as a multi-step depth extrapolation algorithm analogous to phase shift (, ) and finite difference () formulations presented in previous SEP reports.

I begin by briefly reviewing Kirchhoff datuming. I then present the matrix formulation of the operators, and conclude by incorporating these in a multi-step depth extrapolation formulation.

KIRCHHOFF DATUMING The upward continuation of a scalar wavefield can be expressed in integral form as:

$$\EQNLABEL{finallytext} P(0,z,t)=\frac{1}{\pi}\int_{-\infty}^{\infty}dx \frac{z}{v} \frac{1}{r_x} Q(t-\frac{r_x}{v}),$$ (1)

where P(0,z,t) is one trace of the upward continued wavefield at datum elevation z and lateral position x=0 (). The function $Q(t-\frac{r_x}{v})$ is a time delayed filtered version of the input traces on the original datum. r_x is the distance between input and output trace location and v is the propagation velocity. The purpose of the filtering operation is to compensate for the Hankel tail of the output wavelet.

For transformation of a wavefield from one datum to another, a discrete form of equation finallytext is a summation where each input trace P_i, recorded at location i along the lower datum is filtered and delayed by traveltime t_ij. Mathematically:

$$\EQNLABEL{bhtext} {P_j (t)= \sum_i \cos \theta_{ij} A_{ij} Q_i(t-t_{ij})}$$ (2)

A_ij is an amplitude term incorporating spherical divergence and geometric terms, $\theta_{ij}$ is the angle between the normal to the input surface and the straight line traveltime path connecting input location i and output location j, and t_ij is the traveltime between the input and output locations.

Equivalently, the sum can be performed in the frequency domain where each output trace is given by:

$$\EQNLABEL{bhomega} {P_j (\omega)= \sum_i \cos \theta_{ij} A_{ij} e^{i\omega t_{ij}} Q_i(\omega)}.$$ (3)

By making the far-field approximation of the relevant Green's functions, the filtering operation used to compensate for the Hankel tail can be replaced by performing a half-order differential of the data traces (in 2-D). This is computationally attractive because the half-order differential is applied only once and does not depend on the indices i and j.

The algorithm can be applied to zero-offset data or to shot and receiver gathers. Prestack datuming is done by extrapolating the shots and the geophones separately (). The extension to three dimensions is straightforward and has the same algorithmic form. Equations finallytext through bhomega are for upward continuation. For downward continuation the adjoint process is used.

SINGLE OUTPUT TRACE Before writing down the matrix equations for Kirchhoff extrapolation one depth step at a time, I will present the matrix equations for the conventional formulation of Kirchhoff datuming: that is for extrapolation between two surfaces in one step.

Single output trace For illustration I consider the simple case of three input traces and one output trace. Equation bhomega can be expressed in matrix form as

$$\EQNLABEL{kcon} \left[ p(\omega,x_j,z_j) \right] \;=\; \left... ...z_{i=2}) \\ p(\omega,x_{i=3},z_{i=3}) \\ \end{array}\right] ,$$ (4)

where $p(\omega,x_i,z_i)$ is the Fourier transformed input wavefield recorded at datum z_i and $p(\omega,x_j,z_j)$ is an output trace at location (x_j,z_j). A half-order derivative is applied to the input data by the block matrices

$$H \; = \; \left[ \begin{array} {cccccc} \sqrt{-i\omega_1}& &... ... & & &.& \\ & & & & &\sqrt{-i\omega_n} \\ \end{array}\right] .$$ (5)

If the calculation were performed in the time domain and the far-field approximation were not invoked, each of the H matrices would be replaced by a different matrix representing an i,j-variable convolution operator.

The block matrices Z₁₁, Z₁₂, and Z₁₃ perform time shifts and amplitude scaling of individual data traces. For a given time shift t_ij, the time-shift matrix is given by

$$\EQNLABEL{shift} Z_{ij} \; = \; \left[ \begin{array} {cccccc... ...\\ & & & & &A_{ij}e^{i\omega_n t_{ij}} \\ \end{array}\right] .$$ (6)

The A_ij factors are amplitude scaling terms which embody the effects of spherical divergence and obliquity.

Finally, the summation of the shifted input traces is represented by matrix multiplication with the row vector having block identity matrices I as its elements. The summation matrix and the Z_ij matrix can be combined to simplify equation kcon to:

$$\EQNLABEL{kcon2} \left[ p(\omega,x_j,z_j) \right] \;=\; \lef... ...z_{i=2}) \\ p(\omega,x_{i=3},z_{i=3}) \\ \end{array}\right] .$$ (7)

TIME-DOMAIN DATUMING AND MIGRATION For completeness, I show how the matrices differ when datuming is implemented in the time domain.

If Kirchhoff datuming is implemented in the time domain, the Hankel tail can be compensated for by performing an m-long i,j-dependent convolution. Then the matrix H becomes

$$H_{ij} \; = \; \left[ \begin{array} {ccccccc} a_1& & & & & ... ... . \\ & & & & & & . \\ & & & & & &a_m \\ \end{array}\right] .$$ (8)

The time-shift t_ij is represented by a sparse matrix which is non-zero along one particular diagonal (or band) corresponding to the particular shift. For example, a time-shift matrix corresponding to a shift of two time samples is given by:  

$$Z_{ij} \; = \; \left[ \begin{array} {ccccccccc} 0&0& & & & &... ... & & & &.&.&0&0&0 \\ & & & & &0&1&0&0 \\ \end{array}\right] .$$ (9)

If the data are linearly interpolated, the matrix will be banded (see for example Claerbout ).

Kirchhoff time migration Time migration by hyperbola summation is similar to Kirchhoff datuming except that in time migration, the traces are shifted by a time-variant amount given by $t = \sqrt{t_0^2 + x^2/v^2}$.In order for equation kcon2 to represent time migration, the hyperbolic shifts would make the time-shift matrix look something like this:  

$$Z_{ij} \; = \; \left[ \begin{array} {ccccccccc} &1& & & & &... ...\\ & & & & & & & & \\ & & & & & & & & \\ \end{array}\right] .$$ (10)

MULTIPLE OUTPUT TRACES In this section I generalize equation kcon2 to multiple output traces. For illustrative purposes, I write down the matrix equation for three input traces and three output traces. The wavefield at the upper datum z is

$$\EQNLABEL{kmult} \left[ \begin{array} {c} p(\omega,x_{j=1}... ...z_{i=2}) \\ p(\omega,x_{i=3},z_{i=3}) \\ \end{array}\right] .$$ (11)

The time-shift matrices Z_ij correspond to equation shift with ij corresponding to the time shift between input and output locations i and j.

Equation  kmult represents the upward extrapolation of a wavefield between two arbitrarily shaped surfaces.

Downward continuation The downward continuation operator is derived by taking the complex conjugate transpose of the operator in equation kmult. This yields the adjoint operator:

$$\EQNLABEL{kadj} \left[ \begin{array} {ccc} H^*& & \\ &H^... ...2}^* \\ Z_{13}^* & Z_{23}^* & Z_{33}^* \\ \end{array}\right] .$$ (12)

Equation kadj is the downward continuation operator for the extrapolation of a wavefield between two arbitrarily shaped surfaces.

IDEMPOTENCE AND THE HESSIAN

Downward continuation followed by upward continuation is given by the multiplication of the operators in equations kmult and  kadj:

$$\EQNLABEL{pseudo-unitary} \left[ \begin{array} {ccc} Z_{11... ...2}^* \\ Z_{13}^* & Z_{23}^* & Z_{33}^* \\ \end{array}\right] .$$ (13)

If the datuming is performed between level datums with constant velocity, the matrix equations can be simplified by making the following observations and substitutions:

$$\begin{array} {c} Z_{1} = Z_{11} = Z_{22} = Z_{33} , \\ Z_{2... ...} = Z_{23} = Z_{32} ,\\ Z_{3} = Z_{13} = Z_{31} .\\ \end{array}$$

Performing the indicated substitutions, yields the simplified matrices:

$$\EQNLABEL{pseudo2} \left[ \begin{array} {ccc} Z_{1} & Z_{2... ..._{2}^* \\ Z_{3}^* & Z_{2}^* & Z_{1}^* \\ \end{array}\right] .$$ (14)

Using the notation $A_i^2\;=\;Z_{i}Z_{i}^*$and carrying out the multiplication we obtain the following expression for the Hessian matrix:

$$\EQNLABEL{hessian} H \;=\; \left[ \begin{array} {ccc} A_1^2... ...Z_1^*+Z_1Z_2^* & A_3^2 + A_2^2 + A_1^2 \\ \end{array}\right] .$$ (15)

This banded matrix is diagonally dominant with real valued terms on the main diagonal and complex terms along the other diagonals. The datuming operator is not idempotent because of the non-zero terms off the main diagonal. The consequence of this is that if the Kirchhoff datuming operator and its adjoint are applied repeatedly, energy will be lost with each application.

STEP-BY-STEP WAVEFIELD EXTRAPOLATION

The operator in equation kmult can be used to extrapolate a wavefield between two horizontal planes. In this case, each application of the extrapolation algorithm is expressed by the matrix operator W_i which is the operator in equation kmult. For three depth levels, the Kirchhoff datuming algorithm can be formulated as

$$\EQNLABEL{steps} \left[ \begin{array} {c} p(\omega,x_j,z_j)... ...ft[ \begin{array} {c} p(\omega,x_i,z_i) \end{array}\right] ,$$ (16)

The algorithm can be generalized for any number of depth levels and any datum geometry. A, B, and C are geometry matrices which correspond to the datum shape. These matrices represent the action of selecting data traces on the new datum. The matrices A, B, and C have non-zero elements only at certain locations along their diagonals where the wavefield is to be extracted. For a more detailed description of these matrices, see Popovici .

For downward continuation, the conjugate-transpose of equation steps is taken to obtain the adjoint process:

$$\EQNLABEL{stepadj} \begin{array} {c} p(\omega,x_i,z_i)\end{... ...\left[ \begin{array} {c} p(\omega,x_j,z_j)\end{array}\right] ,$$ (17)

where $p(\omega,x,z)$ is data recorded at some irregular datum. The matrices W^*_i correspond to the operator in equation kadj.

DISCUSSION AND CONCLUSIONS For real problems, the operator representations of the preceding section result in huge, albeit sparse, matrices. In the actual computer implementation of Kirchhoff datuming these huge matrix operators are not actually stored. The numerical algorithms apply these linear operators indirectly and much more efficiently.

Equations kmult and kadj represent upward and downward continuation of a scalar wavefield in one step. These are the matrix equivalents of the frequency domain datuming implementation I have presented in previous papers  (, ). A corresponding representation of Berryhill's time domain algorithm can be obtained by changing the H matrices to convolutional operators and casting the rest of the matrices in equations kmult and kadj in the time domain (, ).

Kirchhoff datuming is represented as extrapolation of one depth level at a time in equations steps and stepadj. While computationally less efficient than the one-step method of equations kmult and kadj, it allows greater flexibility in specifying the velocity model. With some approximations, the wavefield can be extrapolated through a v(x,y,z) velocity model. The multi-step operator results in steeper Kirchhoff summation trajectories for each step than the one-step method. This may allow steeper events to be handled with less aperture, therefore gaining some computational advantage.

By examining the simplified datuming operator for the special case of flat datums and constant velocity, I observe that the datuming operator is not idempotent and that the Hessian is diagonally banded.

ACKNOWLEDGEMENTS I am grateful to Stew Levin for critically reviewing this paper and to Jun Ji, Dave Nichols and Mihai Popovici for enlightening discussions.

[SEP,GEOTLE,MISC]