High-order plane-wave destructors

The mathematical basis of the plane-wave destructor filters is the local plane differential equation  

$$\frac{\partial P}{\partial x} + s\,\frac{\partial P}{\partial t} = 0\;,$$ (1)

where P(t,x) is the wave field, and s is the local slope, which may also depend on t and x. In the case of a constant slope, equation (1) has the simple general solution  

$$P(t,x) = f(t - s x)\;,$$ (2)

where f(t) is an arbitrary waveform. Equation (2) is nothing more than a mathematical description of a plane wave.

If the slope s does not depend on the t coordinate, we can transform equation (1) to the frequency domain, where it takes the form of the ordinary differential equation  

$$\frac{d \hat{P}}{d x} + i \omega\,s\, \hat{P} = 0$$ (3)

and has the general solution  

$$\hat{P} (x) = \hat{P} (0)\,e^{i \omega\,s x}\;,$$ (4)

where $\hat{P}$ is the Fourier transform of P. The complex exponential term in equation (4) simply represents a shift of a t-trace according to the slope s and the trace separation x. In the frequency domain, the operator for transforming the trace at position x-1 to the neighboring trace at position x is a multiplication by $e^{i \omega\,s}$. In other words, a plane wave can be perfectly predicted by a two-term prediction-error filter in the F-X domain:  

$$a_0 \, \hat{P} (x) + a_1\, \hat{P} (x-1) = 0\;,$$ (5)

where a₀ = 1 and $a_1 = - e^{-i \omega\,s}$. The goal of predicting several plane waves can be accomplished by cascading several two-term filters. In fact, any F-X prediction-error filter, represented in the Z-transform form as  

$$A(Z_x) = 1 + a_1 Z_x + a_2 Z_x^2 + \cdots + a_N Z_x^N\;,$$ (6)

can be factored into a product of two-term filters:  

$$A(Z_x) = \left(1 - \frac{Z_x}{Z_1}\right)\left(1 - \frac{Z_x}{Z_2}\right) \cdots\left(1 - \frac{Z_x}{Z_N}\right)\;,$$ (7)

where $Z_1,Z_2,\ldots,Z_N$ are the zeroes of polynomial (6). According to equation (5), the phase of each zero corresponds to the slope of a local plane wave multiplied by the frequency. Zeroes that are not on the unit circle carry an additional amplitude gain not included in equation (3).

In order to incorporate time-varying slopes, we need to return back to the time domain and look for an appropriate analog of the phase-shift operator (4) and the plane-prediction filter (5). An important property of plane-wave propagation across different traces is that the total energy of the transmitted wave stays invariant throughout the process. This property is assured in the frequency-domain solution (4) by the fact that the spectrum of the complex exponential $e^{i \omega\,s}$ is equal to one. In the time domain, we can reach an equivalent effect by using an all-pass digital filter. In the Z-transform notation, convolution with an all-pass filter takes the form  

$$\hat{P}_{x+1}(Z_t) = \hat{P}_{x} (Z_t) \frac{B(Z_t)}{B(1/Z_t)}\;,$$ (8)

where $\hat{P}_x (Z_t)$ denotes the Z-transform of the corresponding trace, and the ratio B(Z_t)/B(1/Z_t) is an all-pass digital filter, approximating the time-shift operator (5). In finite-difference terms, equation (8) represents an implicit finite-difference scheme for solving equation (1) with the initial conditions at a constant x. The coefficients of filter B(Z_t) can be determined, for example, by fitting the filter frequency response at small frequencies to the response of the phase-shift operator. The Taylor series technique (equating the coefficients of the Taylor series expansion around zero frequency) yields the expression  

$$B_3(Z_t) = \frac{(1-s)(2-s)}{12}\,Z_t^{-1} + \frac{(2+s)(2-s)}{6} + \frac{(1+s)(2+s)}{12}\,Z_t$$ (9)

for a three-point centered filter B₃(Z_t) and the expression

$$\begin{eqnarray} B_5(Z_t) & = & \frac{(1-s)(2-s)(3-s)(4-s)}{1680}\,Z_t^{-2} + ... ...)(3+s)(4+s)}{420}\,Z_t + \frac{(1+s)(2+s)(3+s)(4+s)}{1680}\,Z_t^2\end{eqnarray}$$ (10)

for a five-point centered filter B₅(Z_t). It is easy to generalize these expressions to longer filters. Figure 1 shows the phase of the all-pass filters B₃(Z_t)/B₃(1/Z_t) and B₅(Z_t)/B₅(1/Z_t) for two values of the slope s in comparison with the exact linear function of equation (4). As expected, the phases fit the exact line at low frequencies, and the accuracy of the approximation increases with the length of the filter.

 

phase
Figure 1 Phase of the implicit finite-difference shift operators in comparison with the exact solution. Left plot corresponds to s=0.5. Right plot; s=0.8.

In two dimensions, equation (8) transforms to the prediction equation analogous to (5) with the 2-D prediction filter  

$$A(Z_t,Z_x) = 1 - Z_x \frac{B(1/Z_t)}{B(Z_t)}\;.$$ (11)

In order to characterize several plane waves, we can cascade several filters of the form (11) in a manner similar to equation (7). In all examples of this paper, I used a modified version of the filter A(Z_t,Z_x), namely the filter  

$$C(Z_t,Z_x) = A(Z_t,Z_x) B(Z_t) = B(Z_t) - Z_x B(1/Z_t)\;,$$ (12)

which avoids the need for polynomial division. In case of the 3-point filter (9), the 2-D filter (12) has exactly six coefficients, with the second t column being a reversed copy of the first column. When filter (12) is used in interpolation problems, it can occasionally cause undesired high-frequency oscillations in the solution, resulting from the near-Nyquist zeroes of the polynomial B(Z_t). The oscillations are easily removed in practice with an appropriate low-pass filtering.

In the next section, I address the problem of estimating the local slope s with the filters of form  (12). Estimating the slope is a necessary step for applying the finite-difference plane-wave filters on real data.