Efficient path-integral diffraction migration

Velocity continuation describes vertical and lateral shifts of time-migrated events under the change of migration velocity. The continuous process in a zero-offset isotropic case is described by equation (Claerbout, 1986; Fomel, 2003b):

$$\frac{\partial^2 P}{\partial t \partial v} + vt\frac{\partial^2 P}{\partial x^2} = 0.$$ (1)

Applying double-Fourier transform after time warpping $\sigma=t^2$, the solution for equation 1 can be expressed as (Fomel, 2003a):

$$\tilde P(k,\Omega,v) = \tilde P_0(k,\Omega)e^{\frac{-ik^2v^2}{16\Omega}},$$ (2)

where $v$ is migration velocity, $\Omega$ is the Fourier dual of $\sigma$ and $k$ is wavenumber. The input zero-offset stack $\tilde P_0(k,\Omega)$ is transformed to a constant velocity time migrated image $\tilde P(k,\Omega,v)$. Burnett and Fomel (2011) provide an extension to the 3-D anisotropic case.

The path-integral formulation creates velocity independent images (Landa et al., 2006) in time domain; the formulation of velocity-weighted path-integral of VC images is:

$$I(t,x) = \int_{v_a}^{v_b} W(v)P(t,x,v)\,\mathrm{d}v,$$ (3)

where $W(v)$ is the velocity-weighting function used to fine-tune velocity constraints formed by $v_a$ and $v_b$.

Efficient workflow of the integral above, proposed by Merzlikin and Fomel (2015), is integrating velocity analytically in the double-Fourier domain. For example, an unweighted integral takes the form:

$$\begin{aligned} \tilde I(k,\Omega) &= \tilde P_0(k,\Omega) \int... ...frac{kv}{4\sqrt{\Omega}}) \right\rvert_{v_a}^{v_b}, \end{aligned}$$

which turns path-integral into an analytical filter in double-Fourier domain. Velocity-weighting function with analytical forms can be included also in this integral.

The efficient path-integral time migration workflow for passive seismic data imaging can be summarized as:

1. Apply different time shifts $\tau$ to passive seismic data to cancel onsets; all operations afterwards are done for each constant $\tau$ slice.
2. Apply $\sigma=t^2$ time warpping to $t$-$x$-$\tau$ domain data.
3. Double-FFT transform from $\sigma$-$x$-$\tau$ domain to $\Omega$-$k$-$\tau$ domain.
4. Apply equation 4 to build the filter in the double-Fourier domain.
5. Apply the filter in last step to $\Omega$-$k$-$\tau$ data.
6. Inverse FFT and inverse time warpping into $t$-$x$-$\tau$ image.

Because the path-integral filtering is implemented in $\Omega$-$k$-$\tau$ domain independently without data communication, all computations can be performed efficiently and in parallel.