Hybrid sparse blind deconvolution: an implementation of SOOT algorithm to real data

The purpose of this study is to evaluate the reflection series to provide seismic data with a higher resolution. Attenuation processes are parts of the Earth property and depend on frequency components. The earth is a low-pass filter (Sheriff and Geldart 1995), which means that when the wave passes through it, the amplitude of the high frequencies has been reduced rapidly compared to the lowest frequencies. Thus, as a side effect of attenuation, the dominant frequency change to lower frequencies and the wavelength is increased and therefore, the vertical resolution has decreased. In order to recover and improve time resolution, the frequency content and bandwidth must be increased, which provides a compressed signal (Yilmaz 2001). Removing the wavelet from the recorded trace is called spiking deconvolution. Because of the complicated state of these types of inverse problems, the solution is often ill-posed. This type of approach has been studied in seismic (Kaaresen and Taxt 1998, Pham et al 2014) and image processing (Stockham et al 1975, Kundur and Hatzinakos 1996, Sroubek et al 2007, Amizic et al 2012, Sroubek and Milanfar 2012, Ahmed et al 2014). To solve deconvolution problems, iterative algorithms must be implemented (Ayers and Dainty 1988). The sparse spiking deconvolution algorithm recovers the reflectivity coefficient by minimizing the l1 norm cost function.

The seismic deconvolution algorithm first was used in 1973 (Dossal and Mallat 2005). This method minimizes the optimization function by the l1 norm (Selesnick 2012). Then many researchers studied the production of images according to the norm l1 (O'Brien et al 1994, Donoho and Elad 2003, Donoho 2006, Goldstein and Osher 2009). If the wavelet and reflectivity are not available, it is called a blind deconvolution, and it is a nonlinear reversible problem and is ill-posed that the wavelet should be thin (Lopes and Barry 2001, Nose-Filho et al 2015, Repetti et al 2015).

Blind deconvolution is applied for image recovery, medical and biological studies, GIS, astronomy, non-destructive testing, and so on (Campisi and Egiazarian 2016). It also applied to seismic data (Kaaresen and Taxt 1998, Luo and Li 1998). The semi-blind deconvolution methods based on φHL regularization (SBD-HL) and based on adaptive φHL regularization (SBD-AHL) are proposed to preserve details (Zhu and Deng 2015, Zhu et al 2015).

The Smoothed-One-Over-Two algorithm (SOOT) was introduced in 2015 (Repetti et al 2015). This method uses the wavelet and reflectivity series as initial values and then solves the deconvolution problem by optimizing the conditions. The mismatched blind deconvolution method uses l1 and l2 norms as regularization parameters. The computation costs of this method are low and provide a reliable vertical resolution (Kazemi and Sacchi 2014, Nose-Filho et al 2015). The minimum entropy algorithm is used to generate a filter to reduce the wavelet effect (Takahata et al 2012). The MM algorithm was introduced in 2012 (Selesnick 2012). This algorithm solves the inverse problem by optimization condition and convex function (Figueiredo et al 2006). During the solution, the algorithm monitors to minimize the difference between the observed and desired data and simultaneously applies the control parameter to the seismic data. The SOOT algorithm is only applied to synthetic seismic data because it requires RC and wavelet coefficients (Pham 2015), which are unknown to real explosive source data. In this research, for the first time, the algorithm used for real seismic data. The reflectivity is obtained by the MM algorithm, and then the wavelet is extracted from the real data and finally, the compressed seismic data with better resolution have been obtained.

The blind seismic deconvolution (or inversion) is provided in Osman and Robinson (1996), Ulrych and Sacchi (2005), where the optimization problem of Repetti et al (2015) solves the problem. The estimation of the original signal can be obtained by solving the following optimization problem:

$$\begin{eqnarray*}&&\hat{h}={\rm{a}}{\rm{r}}{\rm{g}}{\rm{m}}{\rm{i}}{{\rm{n}}}_{x}\{\parallel y-hx\parallel {}_{2}^{2}+\gamma \parallel x\parallel {}_{1}\}.\end{eqnarray*}$$

If the observed data are $y=y(n)$ where $1\leqslant n\leqslant N$ belongs to the space of R^N then:

$$\begin{eqnarray}&&y=\bar{h}\,* \,\bar{x}+w,\end{eqnarray} \tag{ 1 }$$

where $\bar{x}$ is an unknown sparse signal, $\bar{h}$ is the unknown wavelet and w is random noise (Repetti et al 2015).

The goal is to obtain $\hat{x}$ and $\hat{h}$ which is calculated from $\bar{x}$ and $\bar{h}.$ For this matter, the optimization function is introduced.

The optimization function based on equation (1) (Repetti et al 2015) is as follows:

$$\begin{eqnarray}&&F(x,h)=\rho (x,h)+\varphi (x,h)+g(x,h).\end{eqnarray} \tag{ 2 }$$

Then $\rho (x,h)$ is;

$$\begin{eqnarray}&&\rho (x,h)=\displaystyle \frac{1}{2}{\parallel h* x-y\parallel }^{2}.\end{eqnarray} \tag{ 3 }$$

This is a smoothed function, the value of $\varphi (x,h)+g(x,h)$ being used as a regularization parameter (Repetti et al 2015). The value of $\varphi (x,h)$ is written in terms of ${l}_{1}$ and ${l}_{2}$ norms:

$$\begin{eqnarray}&&{l}_{1}(x)=\displaystyle {\sum }_{n=1}^{N}| {x}_{n}| ,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{l}_{2}(x)={\left(\displaystyle {\sum }_{n=1}^{N}{x}_{n}^{2}\right)}^{\displaystyle \frac{1}{2}}.\end{eqnarray} \tag{ 5 }$$

The result of regulation is non-convex and non-smooth (Repetti et al 2015).

For smooth consideration of $\varphi (x,h)$ the constant values $(\alpha ,\eta )$ is added to the equations (Pham 2015):

$$\begin{eqnarray}&&{l}_{1,\alpha }(x)=\displaystyle \sum _{n=1}^{N}(\sqrt{{x}_{n}^{2}+{\alpha }^{2}-\alpha }),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{l}_{2,\eta }(x)=\displaystyle \sum _{n=1}^{N}\sqrt{\displaystyle {\sum }_{n=1}^{N}{x}_{n}^{2}+{\eta }^{2}}.\end{eqnarray} \tag{ 7 }$$

If $\alpha =\eta =0$ then equations (4) and (5) are same as equations (6) and (7).

For transformation of non-convex function to a convex function, the logarithmic properties are used. If $\upsilon$ is assumed a constant the Lagrangian function is defined as:

$$\begin{eqnarray*}&&\mathrm{log}\left(\displaystyle \frac{{l}_{1,\alpha }(x)+\beta }{{l}_{2,\eta }(x)}\right)\leqslant \,\mathrm{log}\,\upsilon .\end{eqnarray*}$$

If $\beta$ is a small value, and the logarithm is not considered, the ${l}_{1}$ to ${l}_{2}$ ratio in this equation is equal or smaller than epsilon.

We assume that the noisy data represented below:

$$\begin{eqnarray*}&&y(n)=(h* x)(n)+w(n),\end{eqnarray*}$$

where $h(n)$ is, the impulse response of an LTI system, $x(n)$ is sparse signal and $w(n)$ is white Gaussian noise. The * denotes convolution. We assume that an LTI system could be described by a recursive equation that has a matrix form as below:

$$\begin{eqnarray*}&&Ay=BX,\end{eqnarray*}$$

where A and B are banded matrix. For example, if the equation is the first order:

$$\begin{eqnarray*}&&y(n)=b(0)x(n)+b(1)x(n-1)-a(1)y(n-1).\end{eqnarray*}$$

Then A and B are:

$$\begin{eqnarray*}{A}=\left[\begin{array}{ccccc}1 & & & & \\ {a}_{1} & 1 & & & \\ & \ddots & \ddots & & \\ & & {a}_{1} & 1 & \\ & & & {a}_{1} & 1\end{array}\right],\\ {B}=\left[\begin{array}{ccccc}{b}_{0} & & & & \\ {b}_{1} & {b}_{0} & & & \\ & \ddots & \ddots & & \\ & & {b}_{1} & {b}_{0} & \\ & & & {b}_{1} & {b}_{0}\end{array}\right].\end{eqnarray*}$$

The output of the system can be rewritten:

$$\begin{eqnarray*}&&y={A}^{-1}Bx=Hx,\end{eqnarray*}$$

where H is equal to:

$$\begin{eqnarray*}&&H={A}^{-1}B.\end{eqnarray*}$$

However A and B are banded matrices; H is not banded matrix (in general inversion of a banded matrix is not banded). Then the model is presented as below:

$$\begin{eqnarray*}&&y=Hx+w.\end{eqnarray*}$$

H and x are inputs to the SOOT algorithm, while the convolution of H and x is input, and the output is a sparse spike signal. Figure 1 shows a comparison between the MM, SOOT and the least squares method applied to the synthetic seismic trace. Figure 2 compares the amplitude spectrum of the MM and SOOT algorithms. As shown in figure 1, the output from the SOOT algorithm provides better compression.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As it is obvious, we can conclude that the MM algorithm provided the spiky result, but compared to the SOOT algorithm, the noise is destroying the signal and cannot provide the desired results such as the SOOT algorithm. To execute the SOOT algorithm, we must obtain source wavelets and reflectivity series as initial values to perform a deconvolution that is not accessible in most real data cases. There is, therefore, an idea to estimate the initial values to perform SOOT deconvolution so that the output from the MM algorithm can be considered as initial values for the SOOT algorithm. The output of MM algorithm is spiky so it can be used as input for the SOOT algorithm. The advantage of combining the SOOT and MM algorithm is that the resolution of the layers is improved. Figure 3 shows the synthetic input data.

Zoom In Zoom Out Reset image size

Figures 3 and 4 give the synthetic results for all the methods, which contain thin layers, a fault plane, and wedges. As shown in figure 3, it demonstrates that the SOOT algorithm can provide a well-compressed section. Figure 4 provides the power spectrum results of figure 3. It illustrates that all frequency values exist after deconvolution. For higher frequencies, the power of the SOOT algorithm shows an improved bandwidth compared to the MM algorithm. This finding shows that SOOT in combination with the MM algorithm could provide satisfactory results for synthetic examples.

Zoom In Zoom Out Reset image size

According to tuning effect, the distance between two adjacent layers must be less than $\displaystyle \frac{\lambda }{4}$ to be distinguishable (Gochioco 1991). When the wavelet is compressed, the frequency bandwidth is increased, and the tuning effect is decreased, so the resolution improves.

To test the ability of the proposed methods, the post-stack data of a reservoir zone are selected. First, MM algorithm applied to the data set and the estimated reflectivity series and a minimum phase wavelet are provided for each trace. The reflectivity is obtained by the MM algorithm, and the wavelets extracted by deconvolution (Bo et al 2013). These values are considered as initial values for the SOOT algorithm. In figure 5 (middle part) the SOOT algorithm is applied to the real data, and the resolution of the layers is improved compared to the input data. In the yellow box and the marked parts are demonstrated that the resolution of these two algorithms improved, and the SOOT algorithm provided best signal compression and hence improved resolution. We can conclude that: first, for a data to be well processed in order to achieve a better, well-compressed interpretation, data could be needed and obtained by the SOOT method. Secondly, the output from the MM algorithm shows a well-compressed data and seems that more compression is not needed; however, looking at the central part of figure 5, (as indicated), there are significant differences between the two methods that illustrate the superiority of the SOOT algorithm. The marked part shows that the tuning effect is reduced due to the SOOT implementation. Thirdly, a disadvantage of the proposed algorithms compared to other regular algorithms is the computational costs. Because this algorithm is a combination of two different deconvolution methods in a chain and at least the run time of the SOOT algorithm is twice the MM method. It is therefore, easy to conclude that the SOOT algorithm is a powerful method with a higher computational cost for real examples, and must be applied to post-stack data to provide a reliable interpretation.

Zoom In Zoom Out Reset image size

In figure 6, the frequency spectrum of the output sections of figure 5 for the SOOT algorithm shows a significant increase over the MM algorithm.

Zoom In Zoom Out Reset image size

The recorded seismic traces are the result of the convolution of the wavelet and reflectivity series; these are the input parameters for the SOOT algorithm and should be sparse and spiky. The results indicate that the SOOT algorithm recovers higher frequencies better than the MM algorithm. The frequency spectrum of the SOOT algorithm with respect to the input data and the output of the MM algorithm show the best resolution improvement. Since the MM algorithm can handle zero phase and minimum phase wavelets, the combination of MM algorithm and the SOOT algorithm is effective. The SOOT algorithm bandwidth of the frequency spectrum increases considerably so that it greatly improves the resolution. The reflection coefficient of the SOOT output is sparser and spiky, consequently it is suitable for thin layer identifications. The SOOT algorithm is a powerful method with a high computational cost for real examples and should only be applied to post-stack data to provide reliable interpretation.