Forward interpolation

Asymptotically pseudo-unitary operators as orthonormal bases

It is interesting to note that many integral operators routinely used in seismic data processing have the form of operator (25) with the Green's function

$$G (t,\mathbf{y};z,\mathbf{x}) = \left\vert\frac{\partial}... ...) \delta \left(z - \theta(\mathbf{x};t,\mathbf{y}) \right)\;.$$ (30)

where we have split the variable $x$ into the one-dimensional component $z$ (typically depth or time) and the $m$-dimensional component $\mathbf{x}$ (typically a lateral coordinate with $m$ equal $1$ or $2$). Similarly, the variable $y$ is split into $t$ and $\mathbf{y}$. The function $\theta$ represents the summation path, which captures the kinematic properties of the operator, and $A$ is the amplitude function. In the case of $m=1$, the fractional derivative $\left\vert\frac{\partial}{\partial t}\right\vert^{m/2}$ is defined as the operator with the frequency response $(i\,\omega)^{m/2}$, where $\omega$ is the temporal frequency (Samko et al., 1993).

The impulse response (30) is typical for different forms of Kirchhoff migration and datuming as well as for velocity transform, integral offset continuation, DMO, and AMO. Integral operators of that class rarely satisfy the unitarity condition, with the Radon transform (slant stack) being a notable exception. In an earlier paper (Fomel, 1996), I have shown that it is possible to define the amplitude function $A$ for each kinematic path $\theta$ so that the operator becomes asymptotically pseudo-unitary. This means that the adjoint operator coincides with the inverse in the high-frequency (stationary-phase) approximation. Consequently, equation (28) is satisfied to the same asymptotic order.

Using asymptotically pseudo-unitary operators, we can apply formula (29) to find an explicit analytic form of the interpolation function $W$, as follows:

$$W (t, \mathbf{y}; t_n, \mathbf{y}_n) = \int\!\!\int G (t, \mathbf{y}; z,\mathbf{x})\, G(t_n,\mathbf{y}_n;z,\mathbf{x})\, d z \, d \mathbf{x} =$$

$$\left\vert\frac{\partial}{\partial t}\right\vert^{m/2} \left\ver... ...\mathbf{y} ) - \theta(\mathbf{x};t_n,\mathbf{y}_n) \right) \, d \mathbf{x}\;.$$ (31)

Here the amplitude function $A$ is defined according to the general theory of asymptotically pseudo-inverse operators as

$$A = {\frac{1}{\left(2\,\pi\right)^{m/2}}}\, {\left\vert F\,\... ...ert\frac{\partial \theta}{\partial t}\right\vert^{(m+2)/4}}\;,$$ (32)

where

$$F = \frac{\partial \theta}{\partial t}\, \frac{\partial^2 \theta}{\pa... ...ial \mathbf{y}}\, \frac{\partial^2 \theta}{\partial \mathbf{x}\, \partial t}\;,$$ (33)

$$\widehat{F} = \frac{\partial \widehat{\theta}}{\partial z}\, \frac{\partial^2 \... ...f{x}}\, \frac{\partial^2 \widehat{\theta}}{\partial \mathbf{y}\, \partial z}\;,$$ (34)

and $\widehat{\theta} (\mathbf{x};t,\mathbf{y})$ is the dual summation path, obtained by solving equation $z=\theta(x;t,y)$ for $t$ (assuming that an explicit solution is possible).

For a simple example, let us consider the case of zero-offset time migration with a constant velocity $v$. The summation path $\theta$ in this case is an ellipse

$$\theta(\mathbf{x};t,\mathbf{y}) = \sqrt{t^2 - \frac{(\mathbf{x}-\mathbf{y})^2}{v^2}}\;,$$ (35)

and the dual summation path $\widehat{\theta}$ is a hyperbola

$$\widehat{\theta}(\mathbf{y};z,\mathbf{x}) = \sqrt{z^2 + \frac{(\mathbf{x}-\mathbf{y})^2}{v^2}}\;.$$ (36)

The corresponding pseudo-unitary amplitude function is found from formula (32) to be (Fomel, 1996)

$$A = {\frac{1}{\left(2\,\pi\right)^{m/2}}}\, {\frac{\sqrt{t/z}}{v^m z^{m/2}}}\;.$$ (37)

Substituting formula (37) into (31), we derive the corresponding interpolation function

$$W (t, \mathbf{y}; t_n, \mathbf{y}_n) = \frac{1}{\left(2\... ...t\,t_n}}{v^{2m} z^{m+1}}\, \delta (z - z_n) \,d \mathbf{x}\;,$$ (38)

where $z = \theta(\mathbf{x};t,\mathbf{y})$, and $z_n = \theta(\mathbf{x};t_n,\mathbf{y}_n)$. For $m=1$ (the two-dimensional case), we can apply the known properties of the delta function to simplify formula (38) further to the form

$$W = {\frac{v}{\pi}}\, {\left\vert\frac{\partial}{\partial... ...v^2 (t + t_n)^2 - (\mathbf{y}-\mathbf{y}_n)^2\right] }}}\;.$$ (39)

The result is an interpolant for zero-offset seismic sections. Like the sinc interpolant in equation (21), which is based on decomposing the signal into sinusoids, equation (39) is based on decomposing the zero-offset section into hyperbolas.

While opening a curious theoretical possibility, seismic imaging interpolants have an undesirable computational complexity. Following the general regularization framework of Chapter , I shift the computational emphasis towards appropriately chosen regularization operators discussed in Chapter . For the forward interpolation method, all data examples in this dissertation use either the simplest nearest neighbor and linear interpolation or a more accurate B-spline method, described in the next section.

Forward interpolation