Implementing implicit finite-difference in the time-space domain using spectral factorization and helical deconvolution

Review of methodology

The two-way acoustic wave equation in one dimension reads:

$$\frac{\partial^2 P}{\partial t^2} = C^2 \frac{\partial^2 P}{\partial x^2}.$$ (A-1)

The central implicit finite-difference approximation I used for the propagation tests was 2nd order in time and 2nd order in space:

$$\frac{P^{t+1}_x -2P^t_x + P^{t-1}_x}{\Delta t^2} = \frac{C^2}{4 \Delta x^2} [ \left( P^{t+1}_{x+1} - 2P^{t+1}_x + P^{t+1}_{x-1} \right)$$

$$+ 2 \left( P^t_{x+1} - 2P^t_x + P^t_{x-1} \right) + \left( P^{t-1}_{x+1} - 2P^{t-1}_x + P^{t-1}_{x-1} \right) ],$$ (A-2)

where $P$ is the pressure wavefield, $t$ and $x$ are the time and space coordinate indices, and $\Delta t$ and $\Delta x$ are the temporal and spatial step sizes. Note that this approximation is based on the Crank-Nicolson method, and so the spatial derivative is balanced between the three time steps: $t-1$ , $t$ , and $t+1$ , where the central time index $t$ has twice the weight of the other two time indices.
In order to propagate the wavefield, the pressure values at time $t+1$ must be equated to the values at times $t$ and $t-1$ . The linear system which must then be solved has the form:

$$\left( \begin{array}{ccccc} U_0 & U_1 & 0 & 0 \\ U_1 & U_0 &... ..._1 \\ P^{t+1}_2 \\ P^{t+1}_3 \\ P^{t+1}_4 \end{array}\right) = \left( \begin{array}{ccccc} V_0 & V_1 & 0 & 0 \\ V_1 & V_0 &... ...rray}{c} P^t_1 \\ P^t_2 \\ P^t_3 \\ P^t_4 \end{array}\right)$$

$$+ \left( \begin{array}{ccccc} W_0 & W_1 & 0 & 0 \\ W_1 & W_0 &... ...1 \\ P^{t-1}_2 \\ P^{t-1}_3 \\ P^{t-1}_4 \end{array}\right).$$ (A-3)

For simplicity, we can combine all the constants into one: $\alpha = \frac{C^2 \Delta t^2}{4 \Delta x^2}$ . The matrix coefficients in equation 3 (the finite-difference weights) are then:
$U_0 = 1 + 2 \alpha,\ \ \ \ \ \ U_1 = -\alpha;$
$V_0 = 2 - 4 \alpha,\ \ \ \ \ \ \ V_1 = 2 \alpha;$
$W_0 = -1 - 2 \alpha,\ \ \ \ W_1 = \alpha.$

In shorter notation, equation 3 reads:

$$U P^{t+1} = V P^t + W P^{t-1}.$$ (A-4)

The solution of this linear system is:

$$P^{t+1} = U^{-1} \left(V P^t + W P^{t-1}\right).$$ (A-5)

To solve this system, we must perform polynomial division. The system is tridiagonal (and easily solvable) only for one dimension. For multiple dimensions, matrix $U$ is block diagonal. Additional non-zero elements appear at a certain offset from the diagonal, making the solution process more complicated. However, using spectral factorization, the finite-difference weights of matrix $U$ (which pertain to time $t+1$ ) can be factorized into a set of causal filter coefficients $u$ and its time reverse $u^{'}$ . Using the helical approach to deconvolution, equation 5 can be recast as:

$$P^{t+1} = (u^{'} u)^{-1} \left(V P^t + W P^{t-1}\right);$$ (A-6)

$$P^{t+1} = u^{-1} (u^{'})^{-1} \left(V P^t + W P^{t-1}\right).$$ (A-7)

Polynomial division is comparable to deconvolution. This means that the polynomial division in equation 5 can be achieved by a set of two deconvolutions of the data by the spectrally factorized coefficients $u$ of matrix $U$ . One deconvolution is done along the data in the reverse direction (application of the adjoint of the filter):

$$y_k = x_k - \sum_{i=1}^{N_u} u_i' y_{k-i},$$ (A-8)

where $u'$ is the time reversed filter coefficients of $u$ . The other deconvolution is done in the forward direction:

$$x_k = y_k - \sum_{i=1}^{N_u} u_i x_{k-i}.$$ (A-9)

I used the SEPlib module polydiv, which uses the helical coordinates to perform the deconvolutions (the polynomial division) in equations 8 and 9.
The wavefield propogation is done by the following sequence:

1. Spectrally factorize the coefficients of matrix $U$ .
2. Multiply the saved wavefield at time $t-1$ by the coefficients of matrix $W$ .
3. Multiply the saved wavefield at time $t$ by the coefficients of matrix $V$ .
4. Sum the results of the previous 2 steps into a result vector.
5. Deconvolve the result vector by the time-reversed factorized filter coefficients $u^{'}$ (eq. 8).
6. Deconvolve the result vector by the factorized filter coefficients $u$ (eq. 9).

Steps 2 - 6 are repeated for each time step. The inputs of the spectral factorization are the finite-difference weights of the matrix $U$ (in Eq. 3), and the outputs are coefficients of a minumum phase filter $u$ . Since I used a constant velocity in all propagation tests, the finite-difference weights are constant, and the filter is stationary.

Implementing implicit finite-difference in the time-space domain using spectral factorization and helical deconvolution