Plane wave prediction in 3-D

Factorizing plane waves

Let us denote the coordinates of a three-dimensional space by $t$, $x$, and $y$. A theoretical plane wave is described by the equation

$$U(t,x,y) = f (t - p_x x - p_y y)\;,$$ (1)

where $f$ is an arbitrary function, and $p_x$ and $p_y$ are the plane slopes in the corresponding direction. It is easy to verify that a plane wave of the form (1) satisfies the following system of partial differential equations:

$$\left\{\begin{array}{rcl}\displaystyle \left(\frac{\partial}... ...rac{\partial}{\partial t}\right) U & = & 0 \end{array}\right.$$ (2)

The first equation in (2) describes plane waves on the $\{t,x\}$ slices. In is discrete form, it is represented as a convolution with the two-dimensional finite-difference filter $\mathbf{A}_x$. Similarly, the second equation transforms into a convolution with filter $\mathbf{A}_y$, which acts on the $\{t,y\}$ slices. The discrete (finite-difference) form of equations (2) involves a blocked convolution operator:

$$\left(\begin{array}{c}\displaystyle \mathbf{A}_x \\ \mathbf{A}_y \end{array}\right) \mathbf{U} = \mathbf{0}\;.$$ (3)

where $\mathbf{U}$ is the model vector corresponding to $U(t,x,y)$.

In many applications, we are actually interested in the spectrum of the prediction filter, which approximates the inverse spectrum of the predicted data. In other words, we deal with the square operator

$$\left(\begin{array}{cc}\displaystyle \mathbf{A}_x^T & \mathb... ... \mathbf{A}_x^T \mathbf{A}_x + \mathbf{A}_y^T \mathbf{A}_y\;.$$ (4)

If we were able to transform this operator to the form $\mathbf{A}^T \mathbf{A}$, where $\mathbf{A}$ is a three-dimensional minimum-phase convolution, we could use the three-dimensional filter $\mathbf{A}$ in place of the inconvenient pair of $\mathbf{A}_x$ and $\mathbf{A}_y$.

The problem of finding $\mathbf{A}$ from its spectrum is known as spectral factorization. It is well understood for 1-D signals (Claerbout, 1976), but until recently it was an open problem in the multidimensional case. Helix transform maps multidimensional filters to 1-D by applying special boundary conditions and allows us to use the full arsenal of 1-D methods, including spectral factorization, on multidimensional problems (Claerbout, 1998b). A problem, analogous to (4), has already occurred in the factorization of the discrete two-dimensional Laplacian operator:

$$\Delta = \nabla^T \nabla = \left(\begin{array}{cc}\displays... ... \mathbf{D}_y \end{array}\right) = \mathbf{H}^T \mathbf{H}\;,$$ (5)

where $\mathbf{D}_x$ and $\mathbf{D}_y$ represent the partial derivative operators along the $x$ and $y$ dimensions, respectively, and the two-dimensional filter $\mathbf{H}$ has been named helix derivative (Claerbout, 1999; Zhao, 1999).

If we represent the filter $\mathbf{A}_x$ with the help of a simple first-order upwind finite-difference scheme

$$\mathbf{U}_{x+1}^t - \mathbf{U}_{x}^t + p_x \left(\mathbf{U}_{x+1}^{t+1} - \mathbf{U}_{x+1}^{t}\right) = 0\;,$$ (6)

then, after the helical mapping to 1-D, it becomes a one-dimensional filter with the $Z$-transform

$$\mathbf{A}_x (Z) = 1 - p_x Z^{N_t + 1} + (p_x - 1) Z^{N_t}\;,$$ (7)

where $N_t$ is the number of samples on the $t$-axis. Similarly, the filter $\mathbf{A}_y$ takes the form

$$\mathbf{A}_y (Z) = 1 - p_y Z^{N_t N_x + 1} + (p_y - 1) Z^{N_t N_x}\;.$$ (8)

The problem is reduced to a 1-D spectral factorization of

$$\mathbf{A}_x (1/Z) \mathbf{A}_x (Z) + \mathbf{A}_y (1/Z) \mathbf{A}_y (Z) = - p_y \frac{1}{Z^{N_t N_x + 1}} + (p_y - 1) \frac{1}{Z^{N_t N_x}} -$$

$$p_x \frac{1}{Z^{N_t + 1}} + (p_x - 1) \frac{1}{Z^{N_t-1}} + \left[p_x (1 - p_x) + p_y (1 - p_y)\right] \frac{1}{Z} +$$

$$2 + p_x (p_x - 1) + p_y (p_y - 1) + \left[p_x (1 - p_x) + p_y (1 - p_y)\right] Z +$$

$$(p_x - 1) Z^{N_t-1} - p_x Z^{N_t + 1} + (p_y - 1) Z^{N_t N_x} - p_y Z^{N_t N_x + 1}\;.$$ (9)

After a minimum-phase factor of (9) has been found, we can use it for 3-D forward and inverse convolution.

All examples in this paper actually use a slightly more sophisticated formula for 2-D plane-wave predictors:

$$\mathbf{A}_x (Z) = 1 + \frac{p_x}{2} (1 - p_x) Z^{N_t - 1} + (p_x^2 - 1) Z^{N_t} - \frac{p_x}{2} (1 + p_x) Z^{N_t + 1}\;.$$ (10)

Formula (10) corresponds to the Lax-Wendroff finite difference scheme (Clapp et al., 1997). It provides a valid approximation of the plane-wave differential equation for $-1 \le p_x \le 1$.

eplane
Figure 1. 3-D plane wave prediction with a 402-point filter. Left: $p_x=0.7$, $p_y=0.5$. Right: $p_x=-0.7$, $p_y=0.5$.

shape Figure 2. Schematic filter shape for a 26-point 3-D plane prediction filter. The dark block represents the leading coefficient. There are 9 blocks in the first row and 17 blocks in the second row.

tplane
Figure 3. 3-D plane wave prediction with a 26-point filter. Left: $p_x=0.7$, $p_y=0.5$. Right: $p_x=-0.7$, $p_y=0.5$.

Figure 1 shows examples of plane-wave construction. The two plots in the figure are outputs of a spike, divided recursively (on a helix) by $\mathbf{A}^{T} \mathbf{A}$, where $\mathbf{A}$ is a 3-D minimum-phase filter, obtained by Wilson-Burg factorization. The factorization was carried out in the assumption of $N_t=20$ and $N_x=20$; therefore, the filter had $N_t N_x +2 = 402$ coefficients. Using such a long filter may be too expensive for practical purposes. Fortunately, the Wilson-Burg method allows us to specify the filter length and shape beforehand. By experimenting with different filter shapes, I found that a reasonable accuracy can be achieved with a 26-point filter, depicted in Figure 2. Plane-wave construction for a shortened filter is shown in Figure 3. The predicted plane wave is shorter and looks more like a slanted disk. It is advantageous to deal with short plane waves if the filter is applied for local prediction of non-stationary signals.

Clapp (2000) has proposed constructing 3-D plane-wave destruction (steering) filters by splitting. In Clapp's method, the two orthogonal 2-D filters $\mathbf{A}_x$ and $\mathbf{A}_y$ are simply convolved with each other instead of forming the autocorrelation (4). While being a much more efficient approach, splitting suffers from induced anisotropy in the inverse impulse response. Figure 4 illustrates this effect in the 2-D plane by comparing the inverse impulse responses of plane-wave filters obtained by spectral factorization and splitting. The splitting response is evidently much less isotropic.

bob
Figure 4. Two-dimensional inverse impulse responses for filters constructed with spectral factorization (left) and splitting (right). The splitting response is evidently much less isotropic.

In the next sections, I address the problem of estimating plane-wave slopes and show some examples of applying local plane-wave prediction in 3-D problems.

Plane wave prediction in 3-D