Exploring three-dimensional implicit wavefield extrapolation with the helix transform

Inverse filter factorization

The conventional way of applying implicit finite-difference schemes reduces to solving a system of linear equations with a sparse matrix. For example, to apply the scheme of equation (11), we can put the filter denominator on the other side of the extrapolation equation, writing it as

$$\left(\mathbf{I} - \frac{a-\beta}{2}\,\mathbf{D}_2\right)... ...bf{I} + \frac{a+\beta}{2}\,\mathbf{D}_2\right) \mathbf{T}_t\;,$$ (14)

where $\mathbf{I}$ is the identity matrix, $\mathbf{D}$ is the convolution matrix for filter (10), and $\mathbf{T}_t$ is the vector of temperature distribution at time level $t$. In the case of two-dimensional extrapolation, the matrix on the left side of equation (14) takes the tridiagonal form

$$\mathbf{A} = \left(\mathbf{I} -c\,\mathbf{D}_2\right) = \l... ...n-1} \\ 0 & & & & -c_{n} & 1 + 2c_{n} \end{array}\right]\;,$$ (15)

where $c = \frac{a-\beta}{2}$, and where, for simplicity, we assume zero-slope boundary conditions. Like any positive-definite tridiagonal matrix, matrix $\mathbf{A}$ can be inverted recursively by an $LU$ decomposition into two bidiagonal matrices. The cost of inversion is directly proportional to the number of vector components. The same conclusion holds for the case of depth extrapolation [equation (13)] with the substitution $c = \frac{\beta + i a}{1-4\,a^2}$.

In the case of a laterally constant coefficient $a$, we can take a different point of view on the tridiagonal matrix inversion. In this case, the matrix $\mathbf{A}_2$ represents a convolution with a symmetric three-point filter $1-c\,D_2(k)$. The $LU$ decomposition of such a matrix is precisely equivalent to filter factorization into the product of a causal minimum-phase filter with its adjoint. This conclusion can be confirmed by the easily verified equality

$$1 + c (Z^{-1} - 2 + Z) = \frac{(1+b)^2}{4}\, \left(1 + \fra... ...b}{1+b} Z\right) \,\left(1 + \frac{1-b}{1+b} Z^{-1}\right)\;,$$ (16)

where $b = \sqrt{1+ 4\,c}$. The inverse of the causal minimum-phase filter $1 + \frac{1-b}{1+b} Z$ is a recursive inverse filter. Correspondingly, the inverse of its adjoint pair, $1 + \frac{1-b}{1+b} Z^{-1}$, is the same inverse filtering, performed in the adjoint mode (backwards in space). In the next subsection, we show how this approach can be carried into three dimensions by applying the helix transform.

Exploring three-dimensional implicit wavefield extrapolation with the helix transform