Multichannel adaptive deconvolution based on streaming prediction-error filter

Theory

According to the theory of predictive deconvolution, time-varying filter coefficients are designed to predict nonstationary seismic data, which can be expressed as a nonstationary autoregressive equation,

$$r(t)=s(t)-\sum_{n=1}^N c_n(t)s(t-n+1-\alpha)=s(t)-\mathbf{S}^T\mathbf{C},$$ (1)

where $\alpha$ is the prediction step, $s(t)$ and $r(t)$ are the original seismic data and prediction error, respectively, $\mathbf{S}^T=[s(t-\alpha),s(t-1-\alpha),\cdots,s(t-N+1-\alpha)]$ , which represents the causal translation of $s(t)$ , all prediction coefficients $\mathbf{C}=[c_1(t),c_2(t),\cdots,c_N(t)]$ are estimated in a time variant manner, and $N$ denotes the length of the predictive filter. The filter coefficients are obtained by solving the least-squares problem:

$$\min_{\mathbf{C}}\parallel{s(t)-\mathbf{S}^T\mathbf{C}}\parallel_2^2,$$ (2)

which is a classical underdetermined problem. The local smoothness constraints with streaming computation (Fomel and Claerbout, 2016) avoids the problem of high computational costs associated with iterative approaches. One can assume that the filter coefficients corresponding to two adjacent data stay close, and the autoregressive equation can be expressed by an overdetermined equation as

$$\begin{aligned}\begin{bmatrix}s(t-\alpha) & s(t-1-\alpha) & \cd... ...1) \cr \vdots \cr \epsilon_tc_N(t-1) \end{bmatrix} \end{aligned}.$$

Equation 3 can be written in terms of a shortened block-matrix notation

$$\begin{bmatrix}\mathbf{S}^T \cr \epsilon_t\mathbf{I} \end{bmatrix... ...bf{C}= \begin{bmatrix}s(t) \cr \epsilon_t\overline{\mathbf{C}}_t \end{bmatrix},$$ (4)

where $\overline{\mathbf{C}}_t=[c_1(t-1),c_2(t-1),\cdots,c_N(t-1)]$ , which represents the previous filter coefficient on time axis, $\epsilon_t$ is the constant scale parameter controlling the variability of two adjacent filter coefficients in the time axis, and $\mathbf{I}$ is the identity matrix. Consider a multichannel seismic data, an extra spatial constraint can introduce to the smoothness of the time-varying filter coefficients along the space axis. The new prediction filter can be expressed in the form of a block matrix as

$$\begin{bmatrix}\mathbf{S}_m^T \cr \epsilon_t\mathbf{I} \cr \epsil... ...n_t\overline{\mathbf{C}}_t \cr \epsilon_x\overline{\mathbf{C}}_x \end{bmatrix},$$ (5)

where $s_m(t)$ and $\mathbf{S}_m^T=[s_m(t-\alpha),s_m(t-1-\alpha),\cdots,s_m(t-N+1-\alpha)]$ represent the multichannel data with space index $m$ and the causal translation of the $m$ th seismic trace, $\epsilon_x$ is the scalar regularization parameter in the space axis, and $\mathbf{C}=[c_1^m(t),c_2^m(t),\cdots,c_N^m(t)]$ is the time-varying filter coefficient with space-varying index $\emph{m}$ . The previous filter $\overline{\mathbf{C}}_t=[c_1^m(t-1),c_2^m(t-1),\cdots,c_N^m(t-1)]$ on the time axis and the previous filter on space axis $\overline{\mathbf{C}}_x=[c_1^{m-1}(t),c_2^{m-1}(t),\cdots,c_N^{m-1}(t)]$ are similar to the current filter $\mathbf{C}=[c_1^m(t),c_2^m(t),\cdots,c_N^m(t)]$ .

Assuming that adjacent filter coefficients are similar, the current filter coefficients at a certain point can be constrained by the adjacent filter coefficients in both time and space directions, and the regularization constraint terms are $\epsilon_t^2\parallel{\mathbf{C}-\overline{\mathbf{C}}_t}\parallel_2^2$ and $\epsilon_x^2\parallel{\mathbf{C}-\overline{\mathbf{C}}_x}\parallel_2^2$ . $\epsilon_t$ and $\epsilon_x$ are weights for regularization constraint terms along time and space directions, respectively, which control the similarity of the adjacent filter coefficients. In this case, the underdetermined problem is transformed into an overdetermined problem. The filter coefficients are calculated by solving the regularized autoregression problem

$$\min_{\mathbf{C}}\parallel{s_m(t)-\mathbf{S}_m^T\mathbf{C}}\paral... ...el_2^2+ \epsilon_x^2\parallel{\mathbf{C}-\overline{\mathbf{C}}_x}\parallel_2^2,$$ (6)

the least-squares solution of equation 6 is

$$\mathbf{C}=(\mathbf{S}_m\mathbf{S}_m^T+\epsilon_t^2\mathbf{I}+ \e... ...S}_m+\epsilon_t^2 \overline{\mathbf{C}}_t+\epsilon_x^2\overline{\mathbf{C}}_x).$$ (7)

The Sherman-Morrison formula in linear algebra (Hager, 1989) is able to directly transform the inverse matrix in equation 7 without iterations:

$$(\mathbf{S}_m\mathbf{S}_m^T+\epsilon_t^2\mathbf{I}+\epsilon_x^2 \... ...bf{S}_m\mathbf{S}_m^T}{\epsilon_t^2+\epsilon_x^2+ \mathbf{S}_m^T\mathbf{S}_m}),$$ (8)

where $\mathbf{S}_m$ is a column vector and $\mathbf{S}_m^T$ is the transpose of $\mathbf{S}_m$ . Substituting equation 8 into  7, the streaming PEF coefficients can be calculated as

$$\mathbf{C}=\overline{\mathbf{C}}+(\frac{s_m(t)-\mathbf{S}_m^T\overline{ \mathbf{C}}}{\epsilon^2+\mathbf{S}_m^T\mathbf{S}_m})\mathbf{S}_m,$$ (9)

where

$$\begin{cases}\epsilon^2=\epsilon_t^2+\epsilon_x^2 \\ \overlin... ...+ \epsilon_x^2\overline{\mathbf{C}}_x}{\epsilon^2} \end{cases}.$$ (10)

Equation 9 shows that the adaptive coefficients get updated by adding a scaled version of the data, and the scale is proportional to the streaming prediction error. Updating the filter coefficients requires only elementary algebraic operation without iteration.

According to the definition of prediction error (equation 1) and prediction coefficients (equation 9), the deconvolution result of streaming PEF can be expressed as

$$r_m(t)=s_m(t)-\mathbf{S}_m^T\mathbf{C}=\frac{\epsilon^2(s_m(t)-\mathbf {S}_m^T\overline{\mathbf{C}})}{\epsilon^2+\mathbf{S}_m^T\mathbf{S}_m}.$$ (11)

In this paper, we select the minimum-phase wavelet as the source wavelet to verify the effectiveness of the proposed deconvolution method. The minimum-phase wavelet can be expressed as

$$w(t)=e^{-(2\pi f_mt/25)^2}sin(2\pi f_mt),$$ (12)

where $f_m$ is the dominant frequency of the wavelet.

The minimum-phase wavelet (figure 1a) and the sequence of reflection coefficients (figure 1b) generate a simple 1D convolution model (figure 1c) , where the wavelet frequencies corresponding to every two reflection coefficients with opposite amplitude are selected to be 45 Hz, 35 Hz, 25 Hz, and 20 Hz, respectively. Figure 1d is the result by using the streaming PEF deconvolution ($N=10$ , $\epsilon_t=0.2$ , and $\alpha=1$ ). The streaming PEF deconvolution method effectively improve the time resolution, however, the relative amplitude relationship among different reflections has been changed, which occur more in predictive deconvolution methods. Notice that the amplitude distortion is related to the dominant frequency of the wavelet: the lower dominant frequency is, the worse amplitude fidelity is shown because the peak amplitude of minimum-phase wavelets is hard to happen in the first sample point.

wave,refl,in,spef1
Figure 1. Analysis of amplitude distortion for streaming PEFs. Minimum-phase wavelet (a), the reflectivity (b), the nonstationary synthetic seismic trace (c), and the deconvolution result by using streaming PEF with constant prediction step (d).

The red line in figure 2 represents the theoretical curve of sample number between the start point and peak-amplitude point of minimum-phase wavelets, which is the function of the dominant frequency. We select an empirical equation to fit the curve as follows:

$$gap=Round(\frac{b}{f\times\Delta t}),$$ (13)

where $Round$ represents a rounding function, $gap$ is the sample number between the start point and the peak-amplitude point, $\Delta t$ is the time interval, and $b$ is a constant. One can determine parameter $b$ from any point in the curve, e.g., $b$ is selected to be 0.232 in figure 2, where the dominant frequency is 1 Hz, the sample interval is 1 ms, and the $gap$ is 232. The blue dotted line calculated by equation 12 reasonably fits the red theoretical curve. Therefore, we further design the streaming PEF deconvolution method with the time-varying prediction step to preserve the relative amplitude relationship of different reflection coefficients. The adaptive prediction step is selected to be the parameter $gap$ , which makes the prediction error after deconvolution close to the maximum amplitude value of the original wavelet. The new streaming deconvolution residual is rewritten as follow:

$$r(t)=s(t)-\sum_{n=1}^N c_n(t)s(t-n+1-\alpha (t)),$$ (14)

where $\alpha(t)$ is the time-varying prediction step, which is determined by a time-varying $gap$ value

$$\alpha(t)=gap(t)=Round(\frac{b}{f(t)\times\Delta t}),$$ (15)

where $f(t)$ is the time-varying local frequency obtained by the shaping regularization (Fomel, 2007a). The instantaneous frequency calculated instantaneously at each signal point sometimes appears noisy and contains physically unreasonable negative frequency, so Fomel (2007a) defined the local frequency by using the shaping regularization (Fomel, 2007b) to constrain the linear inversion problem. When calculating the local frequency, the continuity and smoothness of the local frequency can be controlled only by adjusting a smooth radius parameter, and the local frequency shows better physical meaning than the instantaneous frequency.

gapdif
Figure 2. Parameter $gap$ curve of minimum phase wavelet with different frequencies.

Figure 3a shows the local frequency calculated from the synthetic trace (figure 1c), and the time-varying prediction steps (figure 3b) are obtained according to equation 14, where the parameter $b$ is 0.232 and the time interval is 1 ms. Figure 3c shows the result processed by the streaming PEF deconvolution with the time-varying predictionstep. Compared with original data (dotted line in figure 3c), the proposed method keeps the relative amplitude relationship in the deconvolution result (solid line in figure 3c) at the cost of a lower resolution improvement.

lfe,vvlag,dif
Figure 3. The deconvolution result by using the streaming PEF with time-varying prediction steps. The local frequency (a), time-varying prediction step (b), the deconvolution result with the proposed method (solid line), which is compared with the original trace (dotted line) (c).

The proposed method mainly includes four parameters: the filter length ($N$ ), time constraint factor ( $\epsilon_t$ ), spatial constraint factor ( $\epsilon_x$ ) and constant $b$ . According to the experiment, a reasonable deconvolution results can be obtained when $N\le10$ . The constraint factors $\epsilon_t$ and $\epsilon_x$ are the key parameters for the proposed method. The denominator in equation 9 suggests that $\epsilon_t^2$ and $\epsilon_x^2$ should have the same order of the magnitude as $\mathbf{S}_m^T\mathbf{S}_m$ . Too small a constraint factor would make the deconvolution results unstable, and too large a constraint factor would lead to the filter coefficients not being updated effectively. When the filter coefficients are constrained only in the time direction, it is only necessary to set the spatial constraint factor to zero ( $\epsilon_x=0$ ). In theory, the $b$ value is the ratio of the peak-amplitude time to the period of the wavelet. For minimum-phase wavelets, $b$ is typically less than 0.25.

Like the traditional predictive deconvolution, the relative amplitude relationship of the results that are generated by streaming PEF deconvolution with constant prediction step is not consistent with the original data, so we introduce the time-varying prediction step and derive its empirical formula. After adding the time-varying prediction step, the amplitude of the deconvolution result of the synthetic model (figure 3) tends to be consistent with the true amplitude of the original data. However, the actual seismic data is complex due to the earth absorption attenuation and other factors, so the amplitude of deconvolution results is not true amplitude, but its relative amplitude relationship is closer to the original data.

Since the proposed method can adaptively update the filter coefficients without iteration and characterize the nonstationary properties of the data, both the storage and computational cost of the filters are less than those of the adaptive predictive filtering methods based on the iterative algorithms. Table 1 compares the storage and computational cost of the different methods.

The proposed filter is constrained in both the time and space directions while the filter is still one-dimensional, that is, the multichannel adaptive deconvolution technology is based on the streaming PEF with one-dimensional structure and two-dimensional constraints. To ensure the stability of the calculation, the first boundary trace only uses the time constraint condition to compute the streaming PEF coefficients, and the constraints in both space and time directions are added to the other trace. An extension of the proposed method to 3D is straightforward and provides a fast adaptive multichannel deconvolution implementation for high-dimensional seismic data.

Method	Storage	Cost
Stationary PEF	$O(N_a)$	$O(N_a^2N_t)$
1D nonstationary PEF with iterative algorithm	$O(N_aN_t)$	$O(N_aN_tN_{iter})$
1D streaming PEF	$O(N_a)$	$O(N_aN_t)$
2D streaming PEF	$O(N_aN_t)$	$O(N_aN_tN_x)$

Table 1. Rough cost comparison among the different PEF estimation methods. $N_a$ is the the filter size, $N_t$ is the data length in the time direction, $N_x$ is the data length in the space direction, $N_{iter}$ is the number of iterations.

Multichannel adaptive deconvolution based on streaming prediction-error filter