Discrete matrix formulation of Kirchhoff datuming

For transformation of a wavefield from one datum to another, as depicted in Figure , a discrete form of equation () is  

$$P_j(\omega)\;=\;\sum_i A_{ij}\, \sqrt{-i\omega}\,P_i(\omega) e^{i\omega \tau_{ij}}.$$ (6)

A_ij is an amplitude term incorporating divergence and obliquity, and $\tau_{ij}$ is the traveltime between input location i and output location j. The equivalent time domain equation is  

$$P_j(t)\;=\;\sum_i A_{ij}\, D_{\frac{1}{2}} \ast P_i(t-\tau_{ij}).$$ (7)

If the far-field approximation were not invoked, the $D_{\frac{1}{2}}$ operator would be replaced by an i,j-variable convolution operator. The 3-D equations are algorithmically identical, except that the half-order derivative is replaced by a full derivative, and the amplitude terms

$$A_{ij}\;=\; \frac{\Delta x}{\sqrt{2\pi r_{ij} v}}\, {\rm cos}\,\theta_{ij},$$

are replaced by

$$A_{ij}\;=\; \frac{\Delta x\,\Delta y}{2\pi r_{ij} v}\, {\rm cos}\,\theta_{ij}.$$

$\theta_{ij}$ is the angle between the normal to the input surface and the traveltime path connecting input location i and output location j. The ${\rm cos}\,\theta_{ij}$ factor arises from the identity

$$\frac{\partial r}{\partial n}\;=\; \nabla r \cdot {\bf n}.$$

 

kwed Figure 2 Schematic representation of Kirchhoff upward continuation corresponding to implementation of equation () or equation (). Each input trace Pi is filtered, time shifted according to $\tau_{ij}$, scaled, and summed into output trace Pj.

For the special illustrative case of three input traces, equation () can be expressed in matrix form as

$$\left[ P(x_j,z_j,\omega) \right] \;=\; \left[ \begin{array}... ...,\omega) \\ P(x_{i=3},z_{i=3},\omega) \\ \end{array}\right] ,$$

where $P(x_i,z_i,\omega)$ is the Fourier-transformed input wavefield recorded at datum z_i and $P(x_j,z_j,\omega)$ 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] .$$

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

$$Z_{ij} \; = \; \left[ \begin{array} {cccccc} A_{ij}e^{i\omeg... ... & & & & &A_{ij}e^{i\omega_n \tau_{ij}} \\ \end{array}\right] .$$ (8)

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 the matrix equation to  

$$\left[ P(x_j,z_j,\omega) \right] \;=\; \left[ \begin{array}... ...,\omega) \\ P(x_{i=3},z_{i=3},\omega) \\ \end{array}\right] .$$ (9)

Equation () can be generalized 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 is  

$$\left[ \begin{array} {c} P(x_{j=1},z_{j=1},\omega) \\ P(x... ...,\omega) \\ P(x_{i=3},z_{i=3},\omega) \\ \end{array}\right] .$$ (10)

The time-shift matrices Z_ij correspond to equation () with ij corresponding to the time shift $\tau_{ij}$between input and output locations i and j. In actual implementation, more than three traces are required for summation to represent integration adequately. Three traces are used here merely for ease of exposition, and to establish a pattern which is trivial to extend to more traces.

Equation () represents the upward continuation of a wavefield between two arbitrarily-shaped surfaces. The downward continuation operator is derived by taking the complex conjugate transpose of the operator in equation (). This yields the adjoint operator:  

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

Equation () is the downward continuation operator which is depicted in Figure  for one output trace.

 

kweddown Figure 3 Schematic representation of Kirchhoff downward continuation corresponding to the operator of equation (). Each input trace Pj is filtered, time shifted according to $\tau_{ij}$, scaled, and summed into output trace Pi.

These matrix representations result in huge, albeit sparse, operator 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.

In order to understand why Kirchhoff datuming is much more efficient than Kirchhoff time migration, it is instructive to look at the kinematics of the two processes. Kirchhoff datuming is performed by summing along a trajectory given by

$$t_{\rm datum} = t_{\rm obs} + \sqrt{\frac{z_{\rm datum}^2}{v^2} + \frac{x^2}{v^2}}.$$

This is the equation of a time-shifted hyperbola. The variable $z_{\rm datum}$ is the vertical datuming distance, which remains fixed, so each input trace is shifted by a time-invariant amount. This is illustrated in Figure a for the simple case of extrapolation between two level datums. Time migration is performed by invoking the imaging condition (setting $t_{\rm obs} = 0$) and performing the summation for each image point $t_0 = z_{\rm depth}/v$.The resulting summation trajectory is given by

$$t_{\rm mig} = \sqrt{t_0^2 + \frac{x^2}{v^2}}.$$

The traces are shifted by a time-variant amount because t₀ changes with depth for every image point. This is illustrated in Figure b.

Most of the Kirchhoff datuming results in this dissertation are generated by implementing equation () on Thinking Machine's massively parallel CM5 (Bevc, 1993). The frequency domain is computationally attractive because the static time shift and the convolutional filtering are both efficiently applied by complex multiplication. This frequency domain implementation makes Kirchhoff datuming much more efficient than time domain Kirchhoff migration.

 

trajectory
Figure 4 Comparison of (a) datuming and (b) migration summation trajectories. Datuming is much less computationally intensive because each trace is shifted by a time-invariant amount.

Prestack datuming is performed by using reciprocity to extrapolate shot and receiver gathers (Berryhill, 1984). This can be represented mathematically by the double summation  

$$P({\bf r}_s,{\bf r}_g,t) \; = \; \sum_{\xi_s} \sum_{\xi_g} ... ...}_g, t-\tau({\bf r}_s,{\bf\xi}_s)-\tau({\bf r}_g,{\bf\xi}_g)).$$ (12)

This corresponds to a discrete form of Wiggin's KRK integral (1984). In this representation, the original data are recorded on surface S₁ with surface coordinates $({\bf\xi}_s,{\bf\xi}_g)$. The summations over $\xi_s$ and $\xi_g$ represent integration over the shot and receiver coordinates. $A({\bf r}_*,{\bf\xi}_*)$ and $\tau({\bf r}_*,{\bf\xi}_*)$ represent amplitude and traveltime delay in either shot (* = s) or receiver (* = g) gathers. Because both shot and group gathers are datumed, for 2-D the convolutional operator $D = D_{\frac{1}{2}} \ast D_{\frac{1}{2}}$ corresponds to filtering with a full derivative operator, and for 3-D data D corresponds to filtering with a second-order derivative operator. This sum is performed for all output shot ${\bf r}_s$ and receiver ${\bf r}_g$ coordinates.

So far, this development has been for constant velocity media. In variable velocity media, the calculation of the Green's functions becomes more complicated. The most common way of handling complex velocity structure with Kirchhoff methods is by calculating the traveltime shift $\tau$ using ray tracing or finite-difference solutions to the eikonal equation. I choose to use eikonal traveltimes calculated with van Trier and Symes' (1991) upwind finite-difference algorithm. In Chapter  I demonstrate that this gives excellent results, even for complex velocity models such as the Marmousi synthetic.