next up previous [pdf]

Next: Implementation of methodology with Up: Barak: Implicit helical finite-difference Previous: Introduction

Review of methodology

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

$\displaystyle \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:

$\displaystyle \frac{P^{t+1}_x -2P^t_x + P^{t-1}_x}{\Delta t^2}$ $\displaystyle =$ $\displaystyle \frac{C^2}{4 \Delta x^2}
[ \left( P^{t+1}_{x+1} - 2P^{t+1}_x + P^{t+1}_{x-1} \right)$  
  $\displaystyle +$ $\displaystyle 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:

U_0 & U_1 & 0 & 0 \\
U_1 & U_0 &...
..._1 \\
P^{t+1}_2 \\
P^{t+1}_3 \\
\end{array}\right)\end{displaymath} $\displaystyle =$ \begin{displaymath}\left(
V_0 & V_1 & 0 & 0 \\
V_1 & V_0 &...
P^t_1 \\
P^t_2 \\
P^t_3 \\
  $\displaystyle +$ \begin{displaymath}\left(
W_0 & W_1 & 0 & 0 \\
W_1 & W_0 &...
...1 \\
P^{t-1}_2 \\
P^{t-1}_3 \\
\end{array}\right).\end{displaymath} (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:

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

The solution of this linear system is:

$\displaystyle 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:

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

$\displaystyle 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):

$\displaystyle 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:

$\displaystyle 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.

next up previous [pdf]

Next: Implementation of methodology with Up: Barak: Implicit helical finite-difference Previous: Introduction