next up previous [pdf]

Next: Synthetic data tests Up: Liu and Li: - Previous: introduction


In the $ t$ -$ x$ domain, seismic data always display native nonstationarity. However, one can assume that the new prediction filter $ \mathbf{a}$ stays close to the previous filter on time axis $ \bar{\mathbf{a}}_t$ and the previous filter on space axis $ \bar{\mathbf{a}}_x$ . The autoregression equation of the $ t$ -$ x$ streaming PF can be written in a shortened block-matrix notation (Fomel and Claerbout, 2016) as

$\displaystyle \left[ \begin{array}{c} \mathbf{d}^T\\ \lambda_t \mathbf{I}\\ \la...
...\lambda_t \mathbf{\bar{a}}_t\\ \lambda_x \mathbf{\bar{a}}_x \end{array}\right],$ (1)

where $ \mathbf{I}$ is the identity matrix and $ d(t,x)$ is the given data sample at the point $ (t,x)$ . $ \lambda_t$ and $ \lambda_x$ are the scale parameters controlling the filter variability on the two axes.

Consider a 2D noncausal prediction filter with 20 prediction coefficients $ a_{i,j}$ :

\begin{displaymath}\begin{array}{ccccc} a_{-2,-2}(t,x) & a_{-2,-1}(t,x) & \cdot ...
...a_{2,-1}(t,x) & \cdot & a_{2,1}(t,x) & a_{2,2}(t,x) \end{array}\end{displaymath} (2)

where ``$ \cdot$ '' denotes zero, and the data vector $ \mathbf{d}$ and the filter vector $ \mathbf{a}$ in equation 1 are shown as follows:

$\displaystyle \mathbf{d}=\left[ \begin{array}{c} d_{-M,-N}(t,x)\\ d_{-M+1,-N}(t...
...\ a_{-M+1,-N}(t,x)\\ \ldots\\ a_{M-1,N}(t,x)\\ a_{M,N}(t,x) \end{array}\right],$ (3)

where $ d_{i,j}(t,x)$ represents the translation of $ d(t,x)$ in both time and space directions with time shift $ i$ and space shift $ j$ . $ 2M+1$ and $ 2N$ are the temporal and spatial lengths of the prediction filter, e.g., $ M$ =2 and $ N$ =2 in equation 2.

The least-squares solution of equation 1 is

$\displaystyle \mathbf{a}=(\mathbf{d}\mathbf{d}^T+\lambda^2\mathbf{I})^{-1}(-d(t,x) \mathbf{d}+\lambda^2\mathbf{\bar{a}}),$ (4)


$\displaystyle \lambda^2=\lambda_t^2+\lambda_x^2, \quad\mathbf{\bar{a}}=\frac{\l...
...t^2\mathbf{\bar{a}}_t+ \lambda_x^2\mathbf{\bar{a}}_x}{\lambda_t^2+\lambda_x^2}.$ (5)

Fomel and Claerbout (2016) propose the Sherman-Morrison formula to directly transform the inverse matrix in equation 4 without iterations. The Sherman-Morrison formula is an analytic method for solving the inverse of a special matrix (Hager, 1989). If both matrices $ \mathbf{A}$ and $ \mathbf{I}-\mathbf{V}\mathbf{A}^{-1}\mathbf{U}$ are invertible, then $ \mathbf{A}-\mathbf{U}\mathbf{V}$ is invertible and

$\displaystyle (\mathbf{A}-\mathbf{U}\mathbf{V})^{-1}=\mathbf{A}^{-1}+ \mathbf{A...
...\mathbf{I}- \mathbf{V}\mathbf{A}^{-1}\mathbf{U})^{-1}\mathbf{V}\mathbf{A}^{-1}.$ (6)

In the special case where matrix $ \mathbf{U}$ is a column vector $ \mathbf{u}$ and matrix $ \mathbf{V}$ is a row vector $ \mathbf{v}$ , equation 6 can be rewritten as

$\displaystyle (\mathbf{A}-\mathbf{u}\mathbf{v})^{-1}=\mathbf{A}^{-1}+ \mathbf{A...
...athbf{u}(1- \mathbf{v}\mathbf{A}^{-1}\mathbf{u})^{-1}\mathbf{v}\mathbf{A}^{-1}.$ (7)

The direct derivation of the Sherman-Morrison formula (equation 7) is described in Appendix A.

Applying the Sherman-Morrison formula to equation 4, the $ t$ -$ x$ streaming PF coefficients and prediction error can be calculated as

$\displaystyle \mathbf{a}=\mathbf{\bar{a}}-\left(\frac{d(t,x)+\mathbf{d}^T\mathbf{\bar{a}}} {\lambda^2+\mathbf{d}^T\mathbf{d}}\right)\mathbf{d},$ (8)


$\displaystyle r(t,x)=\frac{\lambda^2(d(t,x)+\mathbf{d}^T\mathbf{\bar{a}})} {\lambda^2+\mathbf{d}^T\mathbf{d}}.$ (9)

For seismic random noise attenuation, we assume the residual of prediction filtering $ r(t,x)$ is the random noise at the point $ (t,x)$ . For calculating 2D $ t$ -$ x$ streaming PFs, we need to store one previous time-neighboring PF, $ \mathbf{\bar{a}}_t$ , and one previous space-neighboring PF, $ \mathbf{\bar{a}}_x$ , both $ \mathbf{\bar{a}}_t$ and $ \mathbf{\bar{a}}_x$ will be used when the stream arrives at its adjacency.

One can compare a streaming PF with a stationary PF. The autoregression equation for a traditional PF takes the following form:

$\displaystyle d(t,x)+\sum_{j=-N,j\neq0}^{j=N}\sum_{i=-M}^{i=M}a_{i,j}d_{i,j}(t,x) \approx 0.$ (10)

The least-squares solution of equation 10 at each point $ (t,x)$ is

$\displaystyle \mathbf{a}=-d(t,x)(\mathbf{d}\mathbf{d}^T)^{-1}\mathbf{d}.$ (11)

The matrix in equation 11 is similar to that in equation 4. The comparison of equation 4 and equation 11 indicates that the results of the streaming PFs become gradually more accurate as the scale parameter $ \lambda$ decreases. However, according to equation 9, a small $ \lambda$ may cause the residual $ r(t,x)$ to also be small, which means that there is too much noise in the signal section. To solve this problem, we use a two-step strategy. First, we choose a relatively large $ \lambda$ to get a large residual $ r(t,x)$ . The first step amounts to an ``over-filtering'', which generates an approximately ``clean'' signal. Next, the signal leakage in the noise section can be extracted by applying signal-and-noise orthogonalization.

We derive the definition of the streaming orthogonalization weight (SOW) from the global orthogonalization weight (GOW) (Chen and Fomel, 2015). Assume that the leaking signal $ \mathbf{s}_1$ has a linear correlation with the estimated signal section $ \mathbf{s}_0$ in the first step,

$\displaystyle \mathbf{s}_1=\omega\mathbf{s}_0,$ (12)

where $ \omega$ is the GOW. The ideal signal $ \mathbf{s}$ and noise $ \mathbf{n}$ can then be estimated by

$\displaystyle \mathbf{s}=\mathbf{s}_0+\omega\mathbf{s}_0,$ (13)

$\displaystyle \mathbf{n}=\mathbf{n}_0-\omega\mathbf{s}_0.$ (14)

Under the assumption that the noise $ \mathbf{n}$ is orthogonal to the signal $ \mathbf{s}$ ,

$\displaystyle \mathbf{n}^T\mathbf{s}=0.$ (15)

Substituting equation 13 and 14 into equation 15, one can get the GOW as

$\displaystyle \omega=\frac{\mathbf{n}_0^T\mathbf{s}_0}{\mathbf{s}_0^T\mathbf{s}_0}.$ (16)

To get the orthogonalization weight for each data value, one possible definition of the SOW is:

$\displaystyle \omega(t,x)=\frac{s(t,x)n(t,x)}{s(t,x)^2}=\frac{n(t,x)}{s(t,x)},$ (17)

where $ \omega(t,x)$ is the SOW for the data point $ d(t,x)$ . $ s(t,x)$ and $ n(t,x)$ represent the signal and noise values at the point $ (t,x)$ , respectively. In the SOPF, $ s(t,x)$ and $ n(t,x)$ correspond to the predictable part $ \mathbf{d}^T\mathbf{a}$ and the prediction residual $ r(t,x)$ produced in the first step. Direct point-by-point division of the values may produce unstable values. One common method to solve this problem is the iterative least-squares method (Chen and Fomel, 2015). Here, we propose a streaming method to calculate the SOW.

Suppose that the SOW gets updated with each new data point $ d(t,x)$ . The new SOW, $ \omega$ , should stay close to the previous time-neighboring SOW $ \bar{\omega}_t$ and the previous space-neighboring SOW $ \bar{\omega}_x$ . Equation 17 can be rewritten as

$\displaystyle \left[ \begin{array}{c} s(t,x)\\ \gamma_t\\ \gamma_x \end{array}\...
...c} n(t,x)\\ \gamma_t\bar{\omega}_t\\ \gamma_x\bar{\omega}_x \end{array}\right],$ (18)

where $ \gamma_t$ and $ \gamma_x$ are the scale parameters controlling the variability on the two axes.

The least-squares solution of equation 18 is

$\displaystyle \omega=\frac{s(t,x)n(t,x)+\gamma^2\bar{\omega}}{s(t,x)^2+\gamma^2},$ (19)


$\displaystyle \gamma^2=\gamma_t^2+\gamma_x^2, \quad\bar{\omega}=\frac{\gamma_t^2\bar{\omega}_t+\gamma_x^2\bar{\omega}_x} {\gamma_t^2+\gamma_x^2}.$ (20)

Applying equation 13, one can get the denoised data value $ \hat{s}(t,x)$ after the SOPF

$\displaystyle \hat{s}(t,x)=s(t,x)+\frac{s(t,x)r(t,x)+\gamma^2\bar{\omega}} {s(t,x)^2+\gamma^2}s(t,x),$ (21)

where $ s(t,x)$ and $ r(t,x)$ are the predictable signals and the prediction residuals calculated in the first step.

We implement the two-step strategy within the streaming method and obtain the denoised data value as each new noisy data value arrives. Table 1 compares the computational cost between $ f$ -$ x$ deconvolution, $ f$ -$ x$ regularized nonstationary autoregression (RNA) (Liu et al., 2012), and the proposed method. In general, the cost of SOPF is minimal.

Method Filter storage Cost
$ f$ -$ x$ deconvolution $ O(N_a)$ $ O(N_a N_f N_x N_{iter})$
$ f$ -$ x$ RNA $ O(N_a N_f N_x)$ $ O(N_a N_f N_x N_{iter})$
$ t$ -$ x$ SOPF $ O(N_a N_t)$ $ O(N_a N_t N_x)$

Table 1. Rough cost comparison between $ f$ -$ x$ deconvolution, $ f$ -$ x$ RNA, and $ t$ -$ x$ SOPF. $ N_a$ is the filter size, $ N_f$ is the frequency size of the data in frequency domain, $ N_x$ is the spatial size of the data, $ N_t$ is the temporal size of the data, and $ N_{iter}$ is the number of iterations.

next up previous [pdf]

Next: Synthetic data tests Up: Liu and Li: - Previous: introduction