next up previous [pdf]

Next: Examples Up: Liu et al.: Noise Previous: Review of - domain

$f$-$x$ domain regularized nonstationary autoregression

Nonstationary autoregression has been developed and used in signal processing (AllenAboutajdine et al., 1996; Bakrim et al., 1994; Izquierdo et al., 2006) and seismic data processing (Sacchi and Naghizadeh, 2009). Fomel (2009) developed a general method of nonstationary autoregression using shaping regularization technology and applied it to multiples subtraction. In this paper, we extend the RNA method to $f$-$x$ domain for complex numbers and apply it to seismic random noise attenuation.

Consider two adjacent seismic traces ${{S}_{n}}(f)$ and ${{S}_{n-1}}(f)$ in $f$-$x$ domain of a seismic section that consists of a single nonlinear event with the slope ${{p}_{n}}$ and varying amplitude ${{B}_{n}}(f)$. Similar to equation 1, we can write

\begin{displaymath}
[{{S}_{n}}(f)=\frac{{{B}_{n}}(f)}{{{B}_{n-1}}(f)}{{e}^{j2\p...
...elta x{{p}_{n}}}}{{S}_{n-1}}(f)={{a}_{1}}(f){{S}_{n-1}}(f).
\end{displaymath} (5)

From equation 5 we can find that the coefficients of AR is the function of space index $n$ and frequency index $f$. Therefore, we can use nonstationary autoregression to describe this problem. Equation 5 describes the relation between two traces. If we consider multiple traces, the nonstationary autoregression can be defined as (Fomel, 2009)

\begin{displaymath}
{{\varepsilon }_{n}}(f)={{S}_{n}}(f)-\sum\limits_{i=1}^{M}{{{a}_{n,i}}(f){{S}_{n-i}}(f)},
\end{displaymath} (6)

where $n$ and $f$are the coordinate of space and frequency, respectively. If considering the situation of non-causal nonstationary autoregression, we can rewrite the Nonstationary autoregression models (equation 6) as
\begin{displaymath}
{{\varepsilon }_{n}}(f)={{S}_{n}}(f)-\sum\limits_{i=1}^{M}{...
...(f)}-\sum\limits_{i=-1}^{-M}{{{a}_{n,i}}(f){{S}_{n-i}}(f)}.
\end{displaymath} (7)

Equation 7 indicates that one trace noise-free in $f$-$x$ domain can be estimated by weighted stacking adjacent traces with the weights ${{a}_{n,i}}(f)$, which is varying along the space and frequency. Note that the difference between equations 4 and 6 is that the coefficients are varying with space coordinate $n$ in equation 7. To obtain the coefficients ${{a}_{n,i}}(f)$ from equation 7, we can transform equation 7 to the following least square problem:
\begin{displaymath}
\min_{a_{n,k}(f)}\vert\vert{{S}_{n}}(f)-\sum\limits_{i=1}^{...
...=-1}^{-M}{{{a}_{n,i}}(f){{S}_{n-i}}(f)}\vert\vert _{2}^{2},
\end{displaymath} (8)

where $\left\Vert {} \right\Vert _{2}^{2}$denotes the squared L-2 norm. Note that both the data ${{S}_{n}}(f)$ and the coefficients ${{a}_{n,i}}(f)$ are in complex numbers domain.

The problem of the minimization in equation 7 is ill-posed because it has more unknown variables than constraints. To obtain the spatial-varying coefficients in nonstationary autoregression, several methods can be employed (AllenAboutajdine et al., 1996). Some of these are related to the expansion of the spatial-varying coefficients in terms of a given sets of orthogonal basis functions and estimation of the coefficients of the expansion by the least-square method (Izquierdo et al., 2006). Naghizadeh and Sacchi (2009) used exponentially weighted recursive least square (EWRLS) to solve the adaptive problem and applied it to seismic trace interpolation. Fomel (2009) proposed regularized nonstationary autoregression in which shaping regularization technology (Fomel, 2007) is used to constrain the smoothness of the coefficients of nonstationary autoregression. In this paper, we also adopt shaping regularization to solve the under-constrained problem equation 8.

With the addition of a regularization term, equation 8 can be written as

\begin{displaymath}
\min_{a_{n,i}{f}}\vert\vert{{S}_{n}}(f)-\sum\limits_{i=-M,i...
...i}}(f){{S}_{n-i}}(f)}\vert\vert _{2}^{2}+R[{{a}_{n,i}}(f)],
\end{displaymath} (9)

where $R$ denotes the shaping regularization operator. Fomel (2009) compared classic Tikhonov’s regularization with shaping regularization in RNA problem. Shaping regularization Fomel (2007) provides a particularly convenient method of enforcing smoothness in iterative optimization schemes. Shaping regularization has clear advantages of a more intuitive selection of regularization parameters and a faster iterative convergence (Fomel, 2009). In the shaping regularization, we assume the initial value for the estimated model ${{a}_{n,i}}(f)$ is zero. If we choose a more appropriate initial value, we can have a fast iterative convergence of the conjugate-gradient iteration in shaping regularization.

Note that the RNA equation in this paper (equation 9) is in complex number domain while the RNA used in multiples subtraction (Fomel, 2009) is in real number domain. Analogous to RNA in real number domain, we force the complex coefficients ${{a}_{n,i}}(f)$ in equation 9 to have a desired behavior, such as smoothness. In shaping regularization technology, we need to choose a shaping operator $\mathbf{S}$ (Fomel, 2007). In this paper, we choose shaping operator $\mathbf{S}$ as Gaussian smoothing with adjustable radius $r$. Fomel (2007) indicated Gaussian smoothing can be implemented by repeated triangle smoothing operator. Both the data ${{S}_{n}}(f)$ and the coefficients ${{a}_{n,i}}(f)$ are complex numbers, but the shaping operator $\mathbf{S}$ is real. Therefore, shaping operator $\mathbf{S}$ is operated in real and imaginary parts of the complex coefficients respectively and the L-2 norm in equation 9 is the norm of complex numbers. When using the algorithm of conjugate-gradient iterative inversion with shaping regularization proposed by Fomel (2007) to solve the complex RNA, we only need to replace transpose of real number by conjugate transpose of complex number.

Once we obtain the complex coefficients of RNA, we can achieve an estimation of signal

\begin{displaymath}
{{\tilde{S}}_{n}}(f)=\sum\limits_{i=-M,i\ne 0}^{M}{{{a}_{n,i}}(f){{S}_{n-i}}(f)}.
\end{displaymath} (10)

When we use $f$-$x$ RNA to noise attenuation, we first select a time window in $t$-$x$ domain and transform the data to $f$-$x$ domain. The usage of time window is to guarantee that the data is approximately stationary in time. The $f$-$x$ NAR can deal with spatial nonstationary data in $f$-$x$ domain. Then we use equation 9 with shaping regularization to compute the coefficients ${{a}_{n,i}}(f)$ and use equation 10 to estimate the signal in $f$-$x$ domain. Finally, we transform data back to the $t$-$x$ domain and repeat for the next time window.

The computational cost of the proposed method is $O\left( {{N}_{d}}{{N}_{f}}{{N}_{iter}} \right)$, where ${{N}_{d}}$ is data size, ${{N}_{f}}$is filter size, ${{N}_{iter}}$is the number of conjugate-gradient iterations. If ${{N}_{iter}}$and ${{N}_{f}}$ are small, this is a fast method. In practical implementation, we can choose the range of computed frequency to reduce the computation cost. In addition, if we can simplify the coefficients not to be frequency-dependent, we can apply RNA in frequency slice. In that case the different-frequency slices can be processed in parallel. And the computation cost and memory requirements can be further reduced.

sin nsin tpefpatch fxpatch npre1 npre2
sin,nsin,tpefpatch,fxpatch,npre1,npre2
Figure 1.
(a) Synthetic data. (b) Noisy data (SNR=1.53). (c) Result of $f$-$x$ domain prediction (SNR=2.53). (d) Result of $t$-$x$ domain prediction (SNR=3.12). (e) Result of $f$-$x$ RNA without constraint of smoothness along the frequency axis (SNR=4.87). (f) Result of $f$-$x$ RNA with ${{r}_{f}}=3$ (SNR=5.06). We decimate the data for display purpose.
[pdf] [pdf] [pdf] [pdf] [pdf] [pdf] [png] [png] [png] [png] [png] [png] [scons]

npar2
npar2
Figure 2.
The real part of coefficients at a given shift ${{a}_{n,i=1}}(f)$.
[pdf] [png] [scons]

fxdiff mpapatch ndiff2
fxdiff,mpapatch,ndiff2
Figure 3.
Difference sections of $f$-$x$ domain prediction (a), $t$-$x$ domain prediction (b), and $f$-$x$ RNA (c).
[pdf] [pdf] [pdf] [png] [png] [png] [scons]

Consider a simple section (Figure 1(a)), which includes one event, 501 traces. The event is obtained by convolution with Ricker wavelet. Both the dip and amplitude of the event are space-varying. The travel time is a sine function and the amplitude of the Ricker wavelet is multiplied by $B(x)=0.2{{(x-2.5)}^{2}}+0.5$. This event is obviously nonstationary in space. We add some random noise to it (Figure 1(b)). The event is greatly contaminated by random noise, especially in the middle part, because the signal of the middle part is poorer than the sideward and the noise levels are the same in the whole section. We use three methods, $f$-$x$ domain prediction, $t$-$x$ domain prediction, and $f$-$x$ domain RNA, to suppress random noise. The $f$-$x$ domain and $t$-$x$ domain prediction methods we used in this paper are discussed by (Abma and Claerbout, 1995). The $f$-$x$ domain prediction is implemented over a sliding window of 20 traces width with 50% overlap and the filter length is 4, $M=2$ in equation 4 and the $t$-$x$ domain prediction is implemented over the same sliding window and the filter length in space and time are 4 and 5 respectively. The $f$-$x$ RNA is implemented with the parameters: the filter length is 4, $M=2$ in equation 8; the smoothing radiuses in space and frequency axes of shaping regularization are respectively 20 and 3, ${{r}_{x}}=20$, ${{r}_{f}}=3$. Comparing the results of $f$-$x$ RNA (Figure 1(f)) with other prediction methods (Figures 4(c) and 4(d), we find that the $f$-$x$ RNA is more effective in random noise attenuation for this simple nonstationary data. From the difference sections (Figure 2), we find that all the three methods remove some effective signals. However, the proposed method removes fewest signals than other two methods. Note that the removed signals by $t$-$x$ and $f$-$x$ prediction methods are bigger in the sideward part than middle part, because the filters are the same in a sliding time window while the amplitude and slope of the event are different.

In order to quantitatively evaluate the effect of denoising between different methods, we use signal-to-noise ratio (SNR) to compare these three methods. The SNR is estimated by

\begin{displaymath}
SNR=10{{\log }_{10}}(\frac{\sum\limits_{t,n}{{{[{{d}_{n}}(t...
...}{\sum\limits_{t,n}{{{[{{d}_{n}}(t)-{{r}_{n}}(t)]}^{2}}}}),
\end{displaymath} (11)

where ${{r}_{n}}(t)$ is the noisy section and ${{d}_{n}}(t)$is the noisy-free section or the section after noise attenuation. In this simple example, the SNR of the input data without processing is 1.53. After noise attenuation, the SNR of $f$-$x$ RNA is 5.06 dB, and the SNRs of $f$-$x$ domain and $t$-$x$ domain prediction are respectively 2.53 dB and 3.12 dB. The $f$-$x$ RNA can improve SNR more greatly.

To test the sensitivity to the constraint of smoothness along the frequency axis, we give the result of $f$-$x$ RNA without constraint of smoothness along the frequency axis in this simple example (Figure 1(e)). The result without constraint along frequency axis (Figure 1(e)) is a little worse than that with constraint along frequency. But both of them are better than the results of conventional $f$-$x$ domain and $t$-$x$ domain prediction. Therefore, to reduce the computation cost, we can simplify the coefficients not to be frequency-dependent when the input seismic data is huge. It is a tradeoff between computation cost and effect of noise attenuation.

Because the dip and amplitude of the event are varying smoothly, $f$-$x$ RNA can predict the signal with smooth coefficients ${{a}_{n,i}}(f)$. To specify the coefficients, we display the real parts of the complex coefficients at a given shift ${{a}_{n,i=1}}(f)$ (Figure 2). The reason of displaying real parts not imaginary parts is that the signs of real parts of coefficients are the same for forward and backward prediction. We find that the real parts of the complex coefficients are smooth. The smoothing radius controls the smoothness of the coefficients. If we only use one adjacent trace to predict the trace, the coefficients should have the expression

\begin{displaymath}
{{a}_{n,i=1}}(f)=\frac{{{B}_{n}}(f)}{{{B}_{n-1}}(f)}{{e}^{j2\pi f\Delta x{{p}_{n}}}}.
\end{displaymath} (12)

From equation 12 we can note that if the dip and amplitude are smoothly varying, the coefficients are smooth. Therefore, we use Gaussian shaping regularization to constrain the coefficients when solving the least squares equation 9.


next up previous [pdf]

Next: Examples Up: Liu et al.: Noise Previous: Review of - domain

2013-11-13