Accelerated plane-wave destruction

Accelerated PWD

Gauss-Newton's iteration searches the solutions for nonlinear equation 6 as follows: Let $p_k$ be the estimated slope at step $k$, with estimating error (or destructive error) $e_k=C(p_k)U(Z_x,Z_t)$. In order to find the correct solution $p_{k+1}$ that minimizes $e_{k+1}$, we need to find the increment $\Delta p_k$ from the local linearization:

$$e_{k+1} = C(p_k)U(Z_x,Z_t)+C'(p_k)U(Z_x,Z_t)\Delta p_k \approx 0\;,$$ (7)

where $C'(p_k)$ is the derivative of $C(p)$ at $p_k$ with respect to $p$.

The iterative algorithm stops when a stationary point or a root of $C(p)U$ is reached. They are:

1. Points where $C'(p_k)U=0$: When $p_k$ satisfies $C'(p)U=0$, then $e_{k+1}=e_k$, and the $\Delta p_k$ dependency in equation 7 is removed, stopping further iterations on $p_k$.
2. Points where $C(p_k)U=0$ and $C'(p)U\neq 0$: In this case, $\Delta p_k=0$, thus $p_{k+1}=p_k$, eliminating the need for further improvements on $p_k$.

The iterative algorithm for equation 6 may converge at different points, depending on the initial point that we chose; $p_0=0$ is a common practical choice for the initial solution. In this case, the iterative algorithm may converge to the least absolute root, which denotes the event with smallest dip angle.

In order to analyze the convergence results, the maximally flat fractional delay filter (Thiran, 1971; Zhang, 2009) is designed with polynomial coefficients:

$$b_k(p)= \frac{(2N)!(2N)!}{(4N)!(N+k)!(N-k)!} \prod_{m=0}^{N-1-k}(m-2N+p) \prod_{m=0}^{N-1+k}(m-2N-p).$$ (8)

Details on how to design the filter can be found in the Appendix.

Since $b_k(p)$ is a polynomial of $p$, expanding it, we get $b_k(p)=\displaystyle{\sum_{i=0}^{2N}c_{ki}p^i}$ and

$$B(Z_t,p)=\sum_{k=-N}^N\sum_{i=0}^{2N}c_{ki}Z_t^{-k}p^i\;.$$ (9)

From equation 8, it is obvious that $b_k(p)=b_{-k}(-p)$, therefore $c_{ki}=(-1)^ic_{-k,i}$ and

$$B(Z_t,p)=B(\frac{1}{Z_t},-p).$$ (10)

Substituting the above two equations, the nonlinear equation 6 becomes a $2N$-th degree polynomial equation for $p$:

$$\sum_{i=0}^{2N}a_ip^i=0,$$ (11)

and the coefficients of the polynomial plane-wave destruction can be expressed as

$$a_i=[1-(-1)^iZ_x]\sum_{k=-N}^Nc_{ki}Z_t^{-k}U,$$ (12)

which says that the coefficients of the polynomial PWD can be obtained by applying a 2D filter on the wavefield $u$. Moreover, the 2D filter can be decoupled into the cascade of two 1D directional filters: the temporal filter $${\sum_{k=-N}^Nc_{ki}Z_t^{-k}}$$ and the spatial filter $1-(-1)^iZ_x$.

In the special case of $N=1$, we get a three-point approximation of $B(Z_t)$. It takes the following form (Fomel, 2002):

$$B(Z_t) = \frac{(1+p)(2+p)}{12}Z_t^{-1}+\frac{(2+p)(2-p)}{6} +\frac{(1-p)(2-p)}{12}Z_t\;.$$ (13)

The plane-wave destruction equation 11 is a quadratic equation. The coefficients $a_i(i=0,1,2)$ can be solved for and expressed as

$$a_i=\frac{1}{12}[1-(-1)^iZ_x]v_i\;,$$ (14)

where $v_i$ are outputs of the following three-point temporal filters:

$$v_0 = 2(Z_t^{-1}+4+Z_t)U\;,$$ (15)

$$v_1 = 3(Z_t-Z_t^{-1})U\;,$$ (16)

$$v_2 = (Z_t^{-1}-2+Z_t)U\;.$$ (17)

In this case, the quadratic plane-wave destruction equation 11 has one stationary point and two roots, which can be analytically expressed as: $\left\{ \frac{-a_1}{2a_2},\quad \frac{-a_1\pm\sqrt{D}}{2a_2} \right\}$, where $D=a_1^2-4a_0a_2$.

The plots in Figure 1 show the convergence process of the iterative algorithm when we choose $p_0=0$ as the starting value. Geometrically when $D\leq 0$, the iteration converges to the stationary point $${\frac{-a_1}{2a_2}}$$, as shown in Figure 1a. When $D>0$, it converges to the least absolute solution of equation 11. Figure 1b and 1c shows the convergence process to the least absolute solution in different cases. We can summarize the convergence result of the iterative algorithm as follows:

$$p=\left\{\begin{array}{ll} \displaystyle{\frac{-a_1}{2a_2}} ... ...rac{-2a_0}{a_1+\sqrt{D}}} & D>0, ~a_1>0\\ \end{array}\right.$$ (18)

As $${\frac{-2a_0}{a_1\pm\sqrt{D}}= \frac{-a_1\pm\sqrt{D}}{2a_2}}$$ in the above equation, we use $${\frac{-2a_0}{a_1\pm\sqrt{D}}}$$ instead of $${\frac{-a_1\pm\sqrt{D}}{2a_2}}$$ to obtain better numerical stability.

When the data is polluted by noise, in order to obtain a robust slope estimation, we can combine the equations in a local window into the following equation set:

$$\tensor F \mathbf p \approx \mathbf g,$$ (19)

where $\tensor F$ is a normalized diagonal matrix and $\mathbf g$ is a vector. Their elements are denuminators and numerators of equation 18 respectively. When we are solving the above equation set by least squares, we can use Tikhonov's regularization (Fomel, 2002) or the shaping regularization (Fomel, 2007a, equation 13) to obtain a smooth solution as follows

$$\mathbf p = \tensor H [\tensor I+\tensor H^T (\tensor F^T\te... ...-\tensor I)\tensor H]^{-1} \tensor H^T\tensor F^T \mathbf g\;,$$ (20)

where $\tensor H$ is an appropriate smoothing operator. In this case, the regularization runs only once, therefore the computational cost is reduced to $O(N_dN_fN_l)$.

In 3D applications, there are two polynomial PWD equations for inline and crossline slopes separately. Note that, using the decoupling, inline and crossline slope estimations can share the temporal filtering results in equations 15$-$17. We can obtain the coefficients of the crossline plane-wave destruction equation as

$$a_i=\frac{1}{12}[1-(-1)^iZ_y]v_i$$ (21)

The five-point or longer approximations of $B(Z_t)$ can achieve higher accuracy. Equation 6 in this case becomes a higher-order polynomial equation (see details in the Appendix), which can be solved numerically. However, there are multiple stationary points, and it is difficult to determine the right one analytically. For applications that need five-point or higher accuracy, we suggest obtaining an initial slope estimation by the proposed three-point method and using it to make the iterative algorithm converge faster (to decrease $N_n$).

Accelerated plane-wave destruction