OVERMIGRATION CAUSED BY THE CRANK-NICOLSON METHOD

Equation (5) can be transformed into the space-frequency domain, resulting in the 15-degree $\omega-x$ time migration equation

$${\partial P\over \partial \tau } = -i\omega P - {v^2\over {-i\omega 2}} {{\partial}^2 P\over \partial x^2}$$ (6)

This partial differential equation can be solved by splitting it (Claerbout, 1985) into the diffraction term

$${\partial P\over \partial \tau } = -i {v^2\over {\omega 2}} {{\partial}^2 P\over \partial x^2}$$ (7)

and the lens term

$${\partial P\over \partial \tau } = -i\omega P$$ (8)

The lens equation can be solved accurately by analytic exponentiation. Here I consider the diffraction equation, which is a partial differential equation (PDE). A PDE can be readily converted to a system of ordinary differential equations (ODEs) by using finite-difference approximations for derivatives in all but one of the dimensions. To do so, I discretize coordinate x with N+1 uniformly spaced grid points, as follows:

$$x_{j}=x_{j-1}+\Delta x, j=0,1,2, \ldots, N$$ (9)

Using the second-order central-difference scheme for the second derivative in x results in

$${\partial P_{j}\over \partial \tau } = \alpha {{ P_{j+1} -2P_{j}+P_{j-1}}\over \Delta x^2}, j=1,2, \ldots, N-1$$ (10)

where $P_{j}=P(\tau,x_{j})$ and $\alpha = -i { v^2 \over {\omega 2}}$. With a further assumption of the boundary condition $P(\tau,0)=P(\tau,L)=0$, I obtain a system of ODEs that can be written in matrix form as

$${d \vec{P}\over d \tau } = A \vec{P}$$ (11)

where P_j's are the elements of the vector $\vec{P}$, and A is the (N-1)*(N-1) tri-diagonal matrix  

$$A = {\alpha \over \Delta x^2} \pmatrix{-2&1&0&\ldots&0&0&0\c... ...dots&\vdots\cr 0&0&0&\ldots&1&-2&1\cr 0&0&0&\ldots&0&1&-2\cr}$$ (12)

This system of ODEs can be solved using any of the numerical methods, such as the Runge-Kutta formulas or multi-step methods.

However, when dealing with systems, we need to be concerned about stiffness. To investigate numerical stiffness, we must study the eigenvalue structure of the matrix A. The eigenvalues of A obtained from a known formula for the eigenvalues of a tri-diagonal matrix with constant entries are

$${\lambda}_{j}={\alpha \over \Delta x^2} (-2+2\cos {{\pi j}\over N}), j=1,2, \ldots, N-1$$ (13)

It is important to keep in mind that all the eigenvalues of A are non-zero and purely imaginary, since $\alpha$ is purely imaginary. The eigenvalue with the smallest magnitude is

$${\lambda}_{1}={\alpha \over \Delta x^2} (-2+2\cos {{\pi }\over N})$$ (14)

For large N, the series expansion for $\cos {{\pi}\over N}$,

$$\cos {{\pi}\over N} = 1 - {1\over {2!}} {({{\pi}\over N})}^2 + {1\over {4!}}{{\pi}\over N}^4 + \ldots$$ (15)

converges rapidly. Retaining the first two terms in the expansion results in

$${\lambda}_{1} \approx -{ {{\pi}^2 \alpha}\over {N^2 \Delta x^2}}$$ (16)

Also, for large N

$${\lambda}_{N-1} \approx -{ {4 \alpha}\over {\Delta x^2}}$$ (17)

Therefore, the ratio of the eigenvalue with the largest modulus to the eigenvalue with the smallest modulus is

$$\left\vert {\lambda}_{N-1} \over {\lambda}_{1} \right\vert \approx { {4N^2}\over {\pi}^2}$$ (18)

Thus for seismic migration, typically $N\geq 500$ traces, this is an extremely stiff system of ODEs. There is severe contradiction in choosing the pseudodepth (time $\tau$) step size of downward extrapolation, so that explicit methods, for example, explicit Euler and Runge-Kutta methods are unworkable.

Standard decoupling using the matrix of eigenvectors S yields

$${d \vec{Q}\over d \tau } = \Lambda \vec{Q}$$ (19)

where $\Lambda=SAS^{-1}$ is a diagonal matrix with the eigenvalues of A on the diagonal, S is the matrix whose columns are the eigenvectors of A, and $Q=S\vec{P}$. Since the equations are uncoupled, the following solution can be readily obtained:

$$Q_{j}(\tau)=Q_{j}(0) e^{{\lambda}_{j}\tau}$$ (20)

Thus, purely imaginary eigenvalues result in the oscillation of the solution that is expected for the wave equation.