Implicit finite difference in time-space domain with the helix transform

Explicit finite difference in 2D using spectral factorization and helical coordinates

In order to test the general methodology and the programing modules, I first tested whether an explicit 2D finite-difference approximation in time-space domain done by spectral factorization and the helix transform properly emulates ``standard'' finite difference (i.e. using a differencing ``star'' which traverses over the data). The difference approximation was 2nd order in time and 2nd order in space:

$$\frac{P^{t+1}_{x,z} -2P^t_{x,z} + P^{t-1}_{x,z}}{\Delta t^2} = C^... ...\Delta x^2} + \frac{P^t_{x,z+1} - 2P^t_{x,z} + P^t_{x,z-1}}{\Delta z^2} \right)$$ (1)

where $P$ is the pressure wavefield, $t$ , $x$ and $z$ are the time and space coordinate indices, and $\Delta t$ , $\Delta x$ , and $\Delta z$ are the temporal and spatial step sizes.

The 2D laplacian representing the spatial derivative's coefficients is:

$$\nabla^2 \approx \left( \begin{array}{ccc} & 1 & 1 & -4 & 1 & 1 & \end{array} \right)$$ (2)

In helical coordinates, the laplacian is represented as a sparse 1-D array:

$$\nabla^2 \approx \left( \begin{array}{ccccccc} 1 & \cdots & 1 & -4 & 1 & \cdots & 1 \end{array} \right)$$ (3)

The dots represent the number of elements required for the wrap-around of the helix along the ``fast'' axis of the 2D grid (see Figure 1). The greater the number of coefficients used by the spectral factorization, the closer the convolution of those coefficients will be to the desired operator, in this case the laplacian in Equation 3. In this experiment, I used the SEPlib wilson module to factorize the laplacian into 28 coefficients. The following is a partial selection of those coefficients:

$$h = \left( \begin{array}{cccccccccc} & & & & 1.0 & -.363 & -.0024... ...ots \cdots & -.0024 & -.0048 & -.113 & -.3114 & & & & & \end{array} \right).$$

Convolving this series with it's adjoint results in the following coefficients:

$$h' * h = \left( \begin{array}{ccccccc} 0.9999771 & \cdots & 0.999... ...000000 & 0.9999771 & \cdots & 0.9999771 \end{array} \right) \approx \nabla^2$$

The remaining coefficients at close proximity to the laplacian coefficients were of 4 or more magnitudes smaller than the ones displayed. The rest were all zeros. This led me to the conclusion that 28 coefficients are sufficient.

I used the SEPlib helicon module to convolve the coefficients in (4) with the data, coiling first in one direction, and then uncoiling in the other direction. The helicon module implements convolution by the linear operator:

$$y_k = x_k + \sum_{i=1}^{N_a} a_i x_{k-i},$$ (4)

where $x$ is the input data vector, $a$ is a causal filter of length $N_a$ , and $y$ is the output vector. The module also implements the adjoint:

$$x_k = y_k + \sum_{i=1}^{N_a} a_i' y_{k+i},$$ (5)

where $a'$ is the time reverse of filter $a$ .

The application of the forward and reverse convolution replicates the act of using the standard finite differencing ``star'' to approximate the spatial derivative. The spatial derivative is then used to estimate the temporal derivative, and forward propagate in time the values of the acoustic wavefield.

With reference to Equation 1, the propagation kernel is:

$$P^{t+1}_{x,z} = 2P^t_{x,z} - P^{t-1}_{x,z} + (h' * h) P^t_{x,z}.$$ (6)

A comparison of propagation of a source in the center of a 2D acoustic wavefield with regular finite differencing and with the helix derivative is shown in Figure 2.

Although qualitatively the two results are similar, the wavefield created by the helix spatial derivative displays some form of dispersion both leading and trailing the main wavelet. Furthermore, the wavefront itself and the dispersion appear to be propagating anisotropically, with a slightly faster component along the positive diagonal. A possible explanation of this outcome is that the filter operator which convolves the data does not have a symmetric angular coverage. Refering to Figure 3, the better coverage of updip angles is apparent in comparison with downdip angles, as the filter is coiled around the data (in the direction of the arrow). If this is indeed the cause of dispersion, one way it can be mitigated is by using a cube finite differencing stencil (Haohuan et al. (2009)).

Another problem which is not apparent in Figure 2 is the wrap-around effect of the helix operator on the data being propagated. Since the helix wraps the filter convolution around one of the dimensions, it will invariably mix the information from one side of the wavefield with that from the other side. The longer and more accurate the filter will be, the more pronounced this effect will be as well. There is no immediate remedy for this inherent property of the helix, except for the standard model-padding solution. However, since the general intention is to use the helix for 2-way wavefield propagation in reverse time migration, one way the wrap around effect can be disregarded is by utilizing random boundaries in the wavefield, as shown in Clapp (2009).

exp-vs-hel-snap
Figure 2. comparison of 2nd order in time and space finite difference Vs. helix derivative after spectrally factorizing the coefficients of the 2nd order spatial derivative and convolving them with the data using helical boundaries. The propagation here is with constant velocity $v = 3000 m/s$ .

helical-aniso
Figure 3. Angular anisotropy of helical filter operator

Implicit finite difference in time-space domain with the helix transform