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 Nx +1 as follows
The operator, , 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, ,means is symmetric, but not Hermitian, so the filter, a1, 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 2Nx+1 can be factored into a pair of (minimum-phase) causal and (maximum-phase) anti-causal filters, each of length Nx+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 and c/v, we precompute the factors of the 1-D helical filters, a1 and a2. 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 a2. The polynomial division is with the filter pair corresponding to a1. 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.