Fourier finite-difference migration

Three-dimensional FFD extrapolation is based on the equation Ristow and Ruhl (1994),  

$$\frac{\partial P}{\partial z}=i \left[ \sqrt{\frac{\omega^2}... ...a_{x,y}^2}{a +b\frac{v^2}{\omega^2} \nabla_{x,y}^2} \right] P,$$ (1)

where v=v(x,y,z) is the medium velocity, c is a reference velocity ($c\leq v$), and a and b are coefficients subject to optimization. The first term describes a simple Gazdag phase-shift that must be applied in the ($\omega, {\bf k}$) domain; the second term describes the first-order split-step correction Stoffa et al. (1990), applied in ($\omega, {\bf x}$); and the third term describes an additional correction that can be applied as an implicit finite-difference operator Claerbout (1985), also applied in ($\omega, {\bf x}$).

In areas with strong lateral velocity variations ($c/v\approx0$), FFD reduces to a finite-difference migration, while in areas of weak lateral velocity variations ($c/v\approx1$), FFD retains the steep-dip accuracy advantages of phase-shift migration. As a full-wave migration method, FFD also correctly handles finite-frequency effects.

For constant lateral velocity, the finite-difference term in equation (1) can be rewritten as the following matrix equation,

$$\begin{eqnarray} \left({\bf I} + \alpha_1 {\bf D} \right) {\bf q}_{z+1} & = & \l... ... {\bf q}_{z} \\ {\bf A}_1{\bf q}_{z+1} & = & {\bf A}_2{\bf q}_{z}\end{eqnarray}$$ (2) (3)

where ${\bf D}$ is a finite-difference representation of the x,y-plane Laplacian, $\nabla_{x,y}^2$, and ${\bf q}_{z}$ and ${\bf q}_{z+1}$ represent the diffraction wavefield at depths z and z+1 respectively . Scaling coefficients, $\alpha_1$ and $\alpha_2$, are complex and depend both on the ratio, $\omega/v$, and the ratio c/v.

The right-hand-side of equation (3) is known. The challenge is to find the vector ${\bf q}_{z+1}$ by inverting the matrix, ${\bf A}_1$.For 2-D problems, only a tridiagonal matrix must be inverted; whereas, for 3-D problems the matrix becomes blocked tridiagonal. For most applications, direct inversion of such a matrix is prohibitively expensive, and so approximations are required for the algorithm to remain cost competitive with other migration methods.

A partial solution is to split the operator to act sequentially along the x and y axes. Unfortunately this leads to extensive azimuthal operator anisotropy, and necessitates expensive additional phase correction operators.