Time migration velocity analysis by velocity continuation

Fourier approach

The change of variable $\sigma = t^2$ transforms equation (1) to the form

$${2 \frac{\partial^2 P}{\partial v \partial \sigma}} + {v \frac{\partial^2 P}{\partial x^2}} = 0\;,$$ (6)

whose coefficients do not depend on the time variables. Double Fourier transform in $\sigma$ and $x$ further simplifies equation (6) to the ordinary differential equation

$$2 i\Omega {{d \hat{P}} \over {d v}} - v k^2 \hat{P} = 0\;,$$ (7)

where the ``frequency'' variable $\Omega$ corresponds to the stretched time coordinate $\sigma$, and $k$ is the wavenumber in $x$:

$$\hat{P}(k,\Omega,v) = \int\!\!\int P(x,t,v) e^{-i \Omega t^2 + i k x} d x dt$$ (8)

Equation (7) has an explicit analytical solution

$$\hat{P} (k,\Omega,v) = \hat{P}_0 (k,\Omega) e^{\frac{i k^2(v_0^2-v^2)}{4\Omega}}\;,$$ (9)

which leads to a very simple algorithm for the numerical velocity continuation. The algorithm consists of the following steps:

1. Input the zero-offset (post-stack) data migrated with velocity $v_0$ (or unmigrated if $v_0=0$).
2. Transform the input from a regular grid in $t$ to a regular grid in $\sigma$.
3. Apply Fast Fourier Transform (FFT) in $x$ and $\sigma$.
4. Multiply by the all-pass phase-shift filter $e^{\frac{i k^2(v_0^2-v^2)}{4\Omega}}$.
5. Inverse FFT in $x$ and $\sigma$.
6. Inverse transform to a regular grid in $t$.

t2
Figure 7. Synthetic seismic data before (left) and after (right) transformation to the $\sigma$ grid.

Figure 7 shows a simple synthetic model of seismic reflection data generated from the model in Figure 2 before and after transforming the grid, regularly spaced in $t$, to a grid, regular in $\sigma$. The left plot of Figure 8 shows the Fourier transform of the data. Except for the nearly vertical event, which corresponds to a stack of parallel layers in the shallow part of the data, the data frequency range is contained near the origin in the $\Omega-k$ space. The right plot of Figure 8 shows the phase-shift filter for continuation from zero imaging velocity (which corresponds to unprocessed data) to the velocity of 1 km/sec. The rapidly oscillating part (small frequencies and large wavenumbers) is exactly in the region, where the data spectrum is zero. It corresponds to physically impossible reflection events.

t2-fft
Figure 8. Left: the real part of the data Fourier transform. Right: the real part of the velocity continuation operator (continuation from 0 to 1 km/s) in the Fourier domain.

The described algorithm is very attractive from the practical point of view because of its efficiency (based on the FFT algorithm). The operation count is roughly the same as in the Stolt migration implemented with equation (4): two forward and inverse FFTs and forward and inverse grid transform with interpolation (one complex-number transform in the case of Stolt migration). The velocity continuation algorithm can be more efficient than the Stolt method because of the simpler structure of the innermost loop (step 4 in the algorithm).

Time migration velocity analysis by velocity continuation