Helical factorization

The blocked-tridiagonal matrix of the 3-D extrapolation, ${\bf A}_1$, represents a two-dimensional convolution operator. Following Rickett et al.'s 1998 approach to factoring the 45 equation, I apply helical boundary conditions Claerbout (1998b) to simplify the structure of the matrix, reducing the 2-D convolution to an equivalent problem in one dimension.

For example, helical boundary conditions allow a two-dimensional 5-point Laplacian filter to be expressed as an equivalent one-dimensional filter of length 2 N_x +1 as follows

$$d = \left[ \begin{array} {ccc} & 1 & \\ 1 & -4 & 1\\ & 1 &... ...tions} (1,\,0, \; ... \; 0,\,1,\,-4,\,1,\,0, \; ... \; 0,\,1).$$

The operator, ${\bf D}$, in equation (2) could represent convolution with this filter; however, I use a more accurate, but equivalent, 9-point filter.

Unfortunately, the complex scale-factor, $\alpha_1$,means ${\bf A}_1$ is symmetric, but not Hermitian, so the filter, a₁, is not an autocorrelation function, and standard spectral factorization algorithms will fail. Fortunately, however, the Kolmogoroff method can be extended to factor any cross-spectrum into a pair of minimum phase wavelets and a delay Claerbout (1998a).

With this algorithm, the 1-D convolution filter of length 2N_x+1 can be factored into a pair of (minimum-phase) causal and (maximum-phase) anti-causal filters, each of length N_x+1. Fortunately, filter coefficients drop away rapidly from either end, and in practice, small-valued coefficients can be safely discarded.

By reconsituting the matrices representing convolution with these filters, I obtain an approximate LU decomposition of the original matrix. The lower and upper-triangular factors can then be inverted efficiently by recursive back-substitution.

While we have only described the factorization for v(z) velocity models, the method can also be extended to handle lateral variations in velocity. For every value of $\omega/v$ and c/v, we precompute the factors of the 1-D helical filters, a₁ and a₂. Filter coefficients are stored in a look-up table. We then extrapolate the wavefield by non-stationary convolution, followed by non-stationary polynomial division. The convolution is with the spatially variable filter pair corresponding to a₂. The polynomial division is with the filter pair corresponding to a₁. The non-stationary polynomial division is exactly analogous to time-varying deconvolution, since the helical boundary conditions have converted the two-dimensional system to one-dimension.