Stability analysis

The starting point for wavefield extrapolation algorithms is an equation that governs the evolution of the wavefield in depth,  

$$\frac{dq}{dz}=-Rq.$$ (50)

The Crank-Nicolson finite-difference scheme makes a numerical approximation of the depth derivative, by equating the two sides of equation () at a point midway between depth steps n and n+1:  

$$\frac{1}{\Delta z} \left( {\bf q}_{n+1} - {\bf q}_{n} \right... ...{1}{2} \, {\bf R}\, \left( {\bf q}_{n+1} + {\bf q}_{n} \right).$$ (51)

Rearranging terms gives the implicit system,  

$$\left( {\bf I} + \frac{\Delta z}{2} {\bf R} \right) {\bf q}_... ...( {\bf I} - \frac{\Delta z}{2} {\bf R} \right) {\bf q}_{n}. \\$$ (52)

An extrapolator in equation () will be stable if $q' \, q$ either decreases or remains constant with depth. Therefore for stable extrapolation,

$$\begin{eqnarray} \frac{d}{dz}(q' \, q) & \leq & 0 \nonumber \\ q' \, \frac{dq}... ...ac{dq'}{dz} \, q & \leq & 0 \nonumber \\ q' (R + R') q & \geq & 0.\end{eqnarray}$$ (53)

This implies that if R+R' is symmetric non-negative definite, then extrapolation with equation () will be unconditionally stable. With a similar proof, Godfrey et al. (1979) showed the same condition applies to the extrapolation matrix ${\bf R}$ under Crank-Nicolson extrapolation with equation (). If ${\bf R}+{\bf R}'$ is symmetric non-negative definite, then stable extrapolation is guaranteed.

While the helical factorization begins with equation (), we really solve the implicit system,

$$\begin{eqnarray} {\bf L}\,{\bf L}^T \; {\bf q}_{n+1} & = & \left( {\bf L}\,{\bf ... ...}, \\ {\bf A} \; {\bf q}_{n+1} & = & {\bf A}^\ast \; {\bf q}_{n},\end{eqnarray}$$ (54) (55)

where ${\bf A} = {\bf L}\,{\bf L}^T \approx {\bf I} + \frac{\Delta z}{2} {\bf R}$.To relate this to equation (), we can premultiply that equation by an invertible matrix ${\bf B}$ to provide a slightly more general Crank-Nicolson extrapolation system with the same stability requirements:  

$${\bf B} \, \left( {\bf I} + \frac{\Delta z}{2} {\bf R} \rig... ...left( {\bf I} - \frac{\Delta z}{2} {\bf R} \right) {\bf q}_{n}.$$ (56)

Equating equations () and () produces a formula for ${\bf R}$ in terms of the helical factorization, ${\bf A}$:

$$\begin{eqnarray} 2 \, {\bf B} & = & {\bf A} + {\bf A}^\ast \\ \Delta z \, {\bf B... ...+ {\bf A}^\ast \right)^{-1} \left( {\bf A} - {\bf A}^\ast \right).\end{eqnarray}$$ (57) (58) (59)

Bulletproof stability requires ${\bf R}+{\bf R}'$ to be symmetric non-negative definite. For the helical factorization this matrix is given by,  

$$\frac{\Delta z}{2} \left( {\bf R}+{\bf R}' \right) = \left( ... ...\bf A}^\ast \right) \left( {\bf A} + {\bf A}^\ast \right)^{-1}.$$ (60)

In the constant velocity case, the matrix ${\bf A}$ represents a stationary filtering operation. Therefore the composite matrices, ${\bf A} + {\bf A}^\ast$ and ${\bf A} - {\bf A}^\ast$ commute with each other. Under this scenario, the matrix ${\bf R}+{\bf R}'$ becomes the zero matrix, which clearly satisfies the non-negative definite criterion required for stability. Constant velocity extrapolation with equation () is therefore unconditionally stable.

Unfortunately, however, if the velocity varies laterally, ${\bf A}$and ${\bf A}^\ast$ no longer commute with each other, and so the composite matrices ${\bf A} + {\bf A}^\ast$ and ${\bf A} - {\bf A}^\ast$ do not commute either. Consequently stable extrapolation cannot be guaranteed. Furthermore, there are no obvious steps we can take to ensure that equation () remains non-negative definite in areas of strong lateral velocity variations. In practice, equation () does indeed encounter stability problems in some areas. Section  illustrates this problem with some examples.

Essentially the problem revolves around the fact that I factor ${\bf I}+ \Delta z/2 \; {\bf R}$ and ${\bf I}- \Delta z/2 \; {\bf R}$, rather than ${\bf R}$ itself. We can ensure our factorization is symmetric non-negative definite, but not the extrapolator itself.