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... ...\frac{\partial}{\partial t}\right)\,U & = & 0\end{array}\right.$$ (2)

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

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

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 \bold{A}_x^T & \bold{A... ...\right) = \bold{A}_x^T \bold{A}_x + \bold{A}_y^T \bold{A}_y\;.$$ (4)

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

The problem of finding $\bold{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}\display... ...d{D}_x \\ \bold{D}_y\end{array}\right) = \bold{H}^T \bold{H}\;,$$ (5)

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

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

$$\bold{U}_{x+1}^t - \bold{U}_{x}^t + p_x \left(\bold{U}_{x+1}^{t+1} - \bold{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  

$$\bold{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 $\bold{A}_y$ takes the form  

$$\bold{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 spectrum factorization of

$$\begin{eqnarray} \nonumber & \bold{A}_x (1/Z) \bold{A}_x (Z) + \bold{A}_y (1/Z)... ...p_x Z^{N_t + 1} + (p_y - 1) Z^{N_t N_x} - p_y Z^{N_t N_x + 1}\;. &\end{eqnarray}$$ (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:  

$$\bold{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 $\bold{A}^T \bold{A}$, where $\bold{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.

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.