Seismic Blind Deconvolution Based on Self-Supervised Machine Learning

Abstract

Seismic deconvolution is a useful tool in seismic data processing. Classical non-machine learning deconvolution methods usually apply quite a few constraints to both wavelet inversion and reflectivity inversion. Supervised machine learning deconvolution methods often require appropriate training labels. The existing self-supervised machine learning deconvolution methods need a given wavelet, which is a non-blind process. To overcome these issues, we propose a blind deconvolution method based on self-supervised machine learning. This method first estimates an initial zero-phase wavelet by smoothing the amplitude spectrum of averaged seismic data. Then, the loss function of self-supervised machine learning is taken as the error between the observed seismic data and the reconstructed seismic data that come from the convolution of phase-rotated wavelet and reflectivity generated by the network. We utilize a residual neural network with long skip connections as the reflectivity inversion network and a fully connected convolutional neural network as the wavelet phase inversion network. Numerical experiments on synthetic data and field data show that the proposed method can obtain reflectivity inversion results with higher resolution than the existing self-supervised machine learning method without given wavelet.

1. Introduction

As a common technique for seismic data processing, deconvolution is important for applications like impedance inversion and high-resolution processing [1,2,3,4]. According to the specific implementation, the present deconvolution methods can be generally divided into three categories: non-machine learning methods, supervised machine learning methods, and self-supervised machine learning methods.

Classical non-machine learning deconvolution methods have gained wide applications in the past few decades. Peacock and Treitel [5] established a theoretical basis for predictive deconvolution. Ulrych [6] proposed homomorphic deconvolution which utilized logarithmic spectrum to obtain an average wavelet and then attained the seismic wavelet and reflectivity. Wiggins [7] proposed minimum entropy deconvolution, which was one of the first studies to explore the concept of sparsity in seismic deconvolution. Based on the prior research findings, a number of novel and influential deconvolution techniques have been proposed in recent years. Wang et al. [8] proposed a blind deconvolution method based on Toeplitz-sparse matrix factorization. This method simultaneously estimated seismic wavelet and reflectivity, and can better maintain the lateral continuity of the seismic profile. Wang et al. [9] developed a high-resolution seismic deconvolution method based on joint sparse representation using the stratigraphic information from well-logging reflectivity and observational seismic data. Xu et al. [10] proposed a blind deconvolution method that exhibited improved fidelity and simplified parameter selection through the application of a compact smoothness constraint on the wavelet and a relative sparsity constraint on the reflectivity.

In the classical methods widely used in the past, scholars often impose L1 or L2 norms on the reflectivity, or apply constraints such as smoothness and finite support on the wavelet. However, these constraints render the objective function challenging to solve, and the parameters of the regularization terms are also difficult to determine. As a type of artificial intelligence, machine learning algorithms enable computers to learn from experience and data without being explicitly programmed. These algorithms provide the key benefits of increased efficiency, accuracy, and scalability [11]. In recent years, machine learning methods have gradually demonstrated promising results in geophysical fields such as denoising, inversion, and fault detection [12,13,14,15,16]. In the field of deconvolution, researchers have also proposed many machine learning-based algorithms. Chen et al. [17] implemented a sequential iterative approach that sequentially inverted the initial reflectivity profile, seismic wavelet, and final reflectivity profile using deep learning techniques. Gao et al. [18] proposed a generalized convolution model based on deep learning and achieved high-fidelity deconvolution results in complex situations. Gao et al. [19] used the information of well data to generate reflectivity sequences, and then constructed labels consistent with the relevant data to obtain well-performing networks. Min et al. [20] proposed a channel attention U-Net to establish a nonlinear mapping relationship between seismic data and reflectivity.

However, most of the aforementioned methods are based on supervised machine learning, which usually relies on a large and representative training dataset comprising low-resolution seismic data and corresponding high-resolution reflectivity labels [21]. If labelled data are limited or unavailable, supervised machine learning typically cannot generalize well [22]. Therefore, the self-supervised machine learning deconvolution method, which utilizes the inherent physical prior contained within the data to guide machine learning towards yielding physically consistent results, is a valuable topic. Chai et al. [23] established a basic framework of self-supervised learning method for deconvolution problem with given wavelet (the SSL-D method for short). This framework incorporates the convolution model into a loss function, then the loss function guides the neural network to produce physically consistent results by minimizing the error between the observed data and the reconstructed data generated by a deep neural network. Wang et al. [21] integrated the sparsity of reflectivity and lateral structural constraint as specific priors into the training process, thereby enhancing the performance of the self-supervised machine learning network to recover more structure details.

As we know, seismic deconvolution is typically a blind process. Seismic blind deconvolution means that both wavelet and reflectivity are unknowns, which leads to high ill-posedness. To overcome this issue, the existing classical non-machine learning studies usually apply quite a few constraints to both wavelet inversion and reflectivity inversion. Excessive constraints make the objective function challenging to solve, and the parameters of the regularization terms are also difficult to determine. As for the machine learning studies, the SSL-D methods mentioned above both need a given wavelet. To our knowledge, self-supervised machine learning has not been explored to solve the problem of seismic blind deconvolution. Therefore, in this paper, we propose a new blind deconvolution method based on self-supervised learning aiming to reduce the reliance on a given wavelet and improve the resolution of reflectivity inversion results. Blind deconvolution is an ill-posed and underdetermined system. We take two main measures to overcome this issue. The first one is estimating a zero-phase wavelet from the seismic data. The second one is applying a relative sparsity process to the reflectivity during the joint inversion of wavelet phase and reflectivity. The steps of the proposed method are briefly described as follows. First, we smooth the spectrum of the averaged seismic data and perform inverse Fourier transform on the smoothed spectrum to obtain an initial zero-phase wavelet. Then, we use a residual neural network and a fully connected convolutional neural network to simultaneously invert reflectivity and wavelet phase (the amplitude spectrum of the wavelet has been obtained at the first step). After that, the loss function for self-supervised optimization is taken as the error between the observed data and the reconstructed data that come from the convolution of phase-rotated wavelet and reflectivity predicted by the network.

The remainder of this paper is organized as follows. Section 2 describes the methodology, which includes the proposed algorithm framework, the network architectures, and the loss function. In Section 3, the experimental results of three datasets are presented. We discuss the limitations of the proposed method in Section 4. Conclusions are given in Section 5.

2. Methodology

Equation (1) describes the relationship between observed seismic data, seismic wavelet, and reflectivity:

(1)${\mathbf{Y}}_{\mathrm{obs}}={\mathbf{T}}_{\mathbf{w}}\mathbf{R},$

where ${\mathbf{Y}}_{\mathrm{obs}}$ is the multi-channel observed seismic data, ${\mathbf{T}}_{\mathbf{w}}$ represents the Toeplitz convolution matrix corresponding to the wavelet $\mathbf{w}$, $\mathbf{R}$ is the multi-channel reflectivity. In this work, we focus on the stationary seismic deconvolution and the wavelet is assumed to be stationary. For the i-th trace in the observed seismic data ${\mathbf{y}}_i$ and the reflectivity ${\mathbf{r}}_i$, Equation (1) can also be written as:

(2)$\left[\begin{array}{c}{y}_{i,1}\\ y_{i,2}\\ y_{i,3}\\ ⋮\\ y_{i,L}\end{array}\right]=\left[\begin{array}{ccccccc}{w}_0& w_1& \cdots & w_{-(m-1)}& 0& \cdots & 0\\ w_1& w_0& \cdots & w_{-(m-2)}& w_{-(m-1)}& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \cdots & 0& w_{m-1}& w_{m-2}& \cdots & w_0& w_1\\ 0& \cdots & 0& w_{m-1}& \cdots & w_1& w_0\end{array}\right]\left[\begin{array}{c}{r}_{i,1}\\ r_{i,2}\\ r_{i,3}\\ ⋮\\ r_{i,L}\end{array}\right]$

where $i=1,2,\dots ,N_{tr}$. $N_{tr}$ is the number of seismic traces, L is the number of time samples of seismic data, and the length of the wavelet is taken as $2m-1$. Within the blind deconvolution framework, the objective is to simultaneously invert both the reflectivity and wavelet from seismic data:

(3)$\mathbf{w},\mathbf{R}=\underset{\widehat{\mathbf{w}},\widehat{\mathbf{R}}}{argmin}{∥{\mathbf{Y}}_{\mathrm{obs}}-{\mathbf{T}}_{\widehat{\mathbf{w}}}\widehat{\mathbf{R}}∥}_2^2.$

2.1. Algorithm Framework

Figure 1 shows the framework of the proposed SSL-BD algorithm, which mainly includes three modules: zero-phase wavelet estimation, wavelet phase inversion network, and reflectivity inversion network. We will introduce these three modules in detail in the following subsections. Here, we briefly explain the steps of the proposed algorithm.

1. Input the observed seismic data ${\mathbf{Y}}_{\mathrm{obs}}$ into the zero-phase wavelet estimation module, wavelet phase inversion network, and reflectivity inversion network to obtain initial zero-phase wavelet ${\mathbf{w}}_0$, initial predicted wavelet phase ${\theta }_{\mathrm{pred}}$, and initial predicted reflectivity ${\mathbf{R}}_{\mathrm{pred}}$, respectively;

2. Rotate the initial zero-phase wavelet ${\mathbf{w}}_0$ with the predicted phase ${\theta }_{\mathrm{pred}}$ to obtain the seismic wavelet ${\mathbf{w}}_{\mathrm{pred}}$, and apply relative sparsity process to the reflectivity ${\mathbf{R}}_{\mathrm{pred}}$ to obtain the sparsity-improved reflectivity ${\tilde{\mathbf{R}}}_{\mathrm{pred}}$;

3. Calculate the convolution of seismic wavelet ${\mathbf{w}}_{\mathrm{pred}}$ and reflectivity ${\tilde{\mathbf{R}}}_{\mathrm{pred}}$ to obtain reconstructed seismic data ${\tilde{\mathbf{Y}}}_{\mathrm{pred}}$.

4. Calculate the mean square error of the input seismic data ${\mathbf{Y}}_{\mathrm{obs}}$ and reconstructed seismic data ${\tilde{\mathbf{Y}}}_{\mathrm{pred}}$ as the network loss, and then update the network’s parameters through gradient back-propagation algorithm to reduce the loss;

5. After the network converges, output the final seismic wavelet and reflectivity.

2.2. Zero-Phase Wavelet Estimation

Figure 2 illustrates the estimation for the initial zero-phase wavelet. We first calculate the averaged seismic data (shown in Figure 2a) and its absolute amplitude spectrum (shown in Figure 2b). Then we smooth the absolute amplitude spectrum using a 7-point local cubic polynomial [24,25]. The smoothed result is shown in Figure 2c. Finally, we perform inverse Fourier transform on the smoothed result to obtain the initial zero-phase wavelet (shown in Figure 2d). The whole process can be expressed by Equation (4) [10]:

(4)${\mathbf{w}}_0=IFFT\left({smooth}^{3L}\left(abs\left(FFT\left(\frac{1}{N_{tr}}{\displaystyle \sum _{i=1}^{N_{tr}}}{\mathbf{y}}_i\right)\right)\right)\right),$

where ${\mathbf{y}}_i$ is i-th trace of the observed seismic data ${\mathbf{Y}}_{\mathrm{obs}}$, $N_{tr}$ is the number of seismic traces, and ${smooth}^k(·)$ means applying k times smoothing. L is the number of time samples of seismic data. ${\mathbf{w}}_0$ is the initial zero-phase wavelet.

After that, we integrate the initial zero-phase wavelet ${\mathbf{w}}_0$ into the joint inversion of the the wavelet phase inversion network and the reflectivity inversion network. Then, we perform phase rotation on the initial zero-phase wavelet ${\mathbf{w}}_0$ according to Equation (5):

(5)${\mathbf{w}}_{\mathrm{pred}}=Re\left[\left({\mathbf{w}}_0+j^{*}hilbert\left({\mathbf{w}}_0\right)\right)·e^{-j{\theta }_{\mathrm{pred}}}\right],$

where $j=\sqrt{-1}$, ${\theta }_{\mathrm{pred}}$ is the wavelet phase predicted by the inversion network, $\mathrm{hilbert}(·)$ means performing the Hilbert transform on the signal, and $Re[·]$ means taking the real part of the complex signal.

2.3. Reflectivity Inversion Network

The reflectivity inversion network structure of 3D seismic data is shown in Figure 3. For the 2D case, we can directly replace the 3D convolution layer with a 2D convolution layer to construct a corresponding 2D network. The structure of the reflectivity inversion network is composed of residual blocks, and has long skip connections (indicated by the long red arrows in Figure 3). For the first convolution layer in the first residual block (indicate by the red arrow), the kernel size is set as $29\times 29\times 29$ to capture features on a larger scale. For the last convolution layer in the last block (indicate by the yellow arrow), the kernel size is set as $1\times 1\times 1$ to achieve cross-channel feature integration. For the other residual blocks, each has two convolution layers with the same kernel size of $3\times 3\times 3$, and each convolution layer is followed by a LeakyReLU activation function. The numbers above the colored blocks are the number of channels of the features extracted from the seismic data. The proposed reflectivity inversion network incorporates long skip connections that bypass certain layers and connect different blocks. These connections create alternative paths for the gradient during back-propagation, enhancing the network’s ability to incorporate fine details in the predictions. As shown in Figure 3, since there are no up-sampling or down-sampling operators in the network, there will be no mismatch in data size during concatenation, so there is no additional restriction for the size of the input data.

The proposed reflectivity inversion network takes 2D/3D post-stack seismic data as input and outputs reflectivity, which will be convolved with the corresponding seismic wavelet during the training process.

2.4. Wavelet Phase Inversion Network

Figure 4 shows the wavelet phase inversion network structure, which includes four convolution layers and two fully connected layers. For the first four convolution layers (indicated by a gray arrow and three orange arrows), the symbols below the colored blocks represent the length of the feature vectors extracted from the seismic data, where the first symbol L represents the length of the seismic trace. $\lfloor \rfloor$ means rounding the length down to the nearest integer. The numbers above the colored blocks represent the number of channels of the feature vectors. For the last two fully connected layers, the numbers below the colored dots mean the number of hidden neurons. The input of the network consists of all traces of seismic data, which are first passed through a convolution layer with a kernel-size of $29\times 1$ and a stride of 14 (indicated by the gray arrow). The subsequent convolution layers illustrated by the orange arrows in Figure 4 are characterized by a kernel size of $4\times 1$ and a stride of 2. The purple arrow represents flattening the feature maps into a one-dimensional vector of length $L_{\mathrm{sequence}}\times n_{\mathrm{chamel}}$ (set as $L/8\times 64$ in this paper). The subsequent fully connected layer further extracts this vector into a single value, which represents the phase of the wavelet ${\theta }_{\mathrm{pred}}$.

2.5. Loss Function

Before introducing the loss function, we first explain the relative sparsity process on the reflectivity predicted by the reflectivity inversion network. Due to the contamination of noise and the band-limited nature of seismic data, the blind deconvolution problem is ill-posed and underdetermined. Numerous deconvolution results can fit the same seismic data well. Hence, an additional process is generally needed to achieve the optimal solution. Inspired by [10,23], we impose a relative sparsity constraint with a variable relative sparsity that increases with the training epochs defined by Equation (7) to obtain the optimal reflectivity inversion result:

(6)${\tilde{\mathbf{R}}}_{\mathrm{pred}}=\left\{\begin{array}{c}{\mathbf{R}}_{\mathrm{pred}},\left|{\mathbf{R}}_{\mathrm{pred}}\right|\ge max\left(abs\left({\mathbf{R}}_{\mathrm{pred}}\right)\right)\times \mu ,\hfill \\ 0,\left|{\mathbf{R}}_{\mathrm{pred}}\right|<max\left(abs\left({\mathbf{R}}_{\mathrm{pred}}\right)\right)\times \mu ,\hfill \end{array}\right.$

where ${\mathbf{R}}_{\mathrm{pred}}$ is the output of the reflectivity inversion network, ${\tilde{\mathbf{R}}}_{\mathrm{pred}}$ is the sparsity-promoted reflectivity, and $\mu$ is defined as follows:

(7)$\mu =\left\{\begin{array}{c}0.01,r_{\mathrm{epoch}}<0.1,\hfill \\ \frac{r_{\mathrm{epoch}}-0.1}{0.9-0.1}\times 0.1+\frac{0.9-r_{\mathrm{epoch}}}{0.9-0.1}\times 0.01,0.1\le r_{\mathrm{epoch}}\le 0.9,\hfill \\ 0.1,r_{\mathrm{epoch}}>0.9,\hfill \end{array}\right.$

where $r_{\mathrm{epoch}}$ is the ratio of the current epoch to the number of all epochs during the training process.

Then, the sparsity-promoted reflectivity ${\tilde{\mathbf{R}}}_{\mathrm{pred}}$ is convolved with the phase-rotated seismic wavelet ${\mathbf{w}}_{\mathrm{pred}}$ to generate the reconstructed seismic data ${\tilde{\mathbf{Y}}}_{\mathrm{pred}}$. And then, ${\tilde{\mathbf{Y}}}_{\mathrm{pred}}$ further participates in the calculation of the loss function with observed seismic data ${\mathbf{Y}}_{\mathrm{obs}}$:

(8)$\begin{array}{cc}\hfill loss\left({\mathbf{Y}}_{\mathrm{obs}},{\tilde{\mathbf{Y}}}_{\mathrm{pred}}\right)& ={∥{\mathbf{Y}}_{\mathrm{obs}}-{\tilde{\mathbf{Y}}}_{\mathrm{pred}}∥}_2^2={∥{\mathbf{Y}}_{\mathrm{obs}}-{\mathbf{T}}_{{\mathbf{w}}_{\mathrm{pred}}}{\tilde{\mathbf{R}}}_{\mathrm{pred}}∥}_2^2.\hfill \end{array}$

It is worth noting that our proposed method obtains reflectivity and seismic wavelet from seismic data simultaneously, which is a completely blind process.

3. Experiments

One synthetic dataset and two field datasets were used in this study, which are denoted as the part of the Marmousi2 model, the CNPC data, and the iErsk3D data, respectively. For comparison, we conduct the geophysics-steered self-supervised learning method proposed by Chai et al. [23], which is a mature and peer-reviewed non-blind deconvolution method, also known as SSL-D. The experimental settings of the comparison method are consistent with the code published in [23]. We train our networks for a maximum of 10,000 epochs for both synthetic data and field data. We use the Adam optimizer with a constant learning rate of 10⁻⁵. The training process is based on the Pytorch framework and uses the GeForce RTX 3090 for GPU acceleration. For the three datasets used in this study, the runtimes of the proposed method and the comparison method are presented in Table 1. Although our proposed SSL-BD method is several times slower than the SSL-D method, the corresponding efficiency reduction should be reasonable because our method can avoid the requirement of a given wavelet for the SSL-D method.

We utilize the evaluation metrics presented in Table 2 to quantitatively compare the results obtained from different methods [26].

3.1. Synthetic Experiment

In this section, we demonstrate the effectiveness of the proposed SSL-BD method on synthetic data. We obtain the synthetic seismic data by convolving the reflectivity profile generated from a part of the Marmousi2 model [10] and the Ricker wavelet with a dominant frequency of 30 Hz and a phase of 30°. The reflectivity profile has 192 seismic traces, each trace has 800 time sampling points, and the sampling interval is 1 ms.

Figure 5 illustrates the decreasing trend of the loss value of the proposed SSL-BD method, indicating the convergence of the training process associated with the synthetic data. Figure 6 shows a comparison between the true wavelet and the wavelet estimated by the proposed SSL-BD method, in which the black solid line is the true wavelet and the red dotted line is the wavelet estimated by the SSL-BD method. The estimated wavelet is almost the same as the true wavelet. To further show how the wavelet estimation process behaves across the frequency band, we present the amplitude and phase spectra of the true wavelet and the estimated wavelet in Figure 7. Figure 8c,d show the predicted reflectivity profile of SSL-D [23] and the proposed SSL-BD method, respectively. It is evident that the predicted reflectivity profile exhibits a high degree of resemblance to the true reflectivity profile. Results of the evaluation metrics for the synthetic data are presented in Table 3, where the better indicator is highlighted in bold. From Table 3, we can observe that the proposed SSL-BD method achieves the best performance in all indicators compared to the SSL-D method. Figure 9 shows the detailed comparison of a single trace located in the 100-th trace of the true reflectivity profile and the predicted reflectivity profile. To investigate the sensitivity of the proposed SSL-BD method to noise, we conduct a series of reflectivity inversion experiments with varying noise levels. The corresponding evaluation metrics are tabulated in Table 4 and the specific results are depicted in Figure 10.

3.2. Field Experiment I

To demonstrate the practicality of the SSL-BD method, we present two field examples. The first field dataset comes from CNPC, as shown in Figure 11. It comprises 976 seismic traces, and each trace has 688 sampling points with an interval of 1 ms.

Figure 12a shows the wavelet extracted from the CNPC seismic data. Figure 12b shows the wavelet estimated by the proposed SSL-BD method. Figure 13 shows the comparisons of wavelet amplitude spectra and phase spectra between the wavelet extracted from the CNPC data and the wavelet estimated by the proposed SSL-BD method.

Figure 14 shows a detailed comparison among the impedance curve from the well, a reflectivity calculated using the impedance curve, and the reflectivity calculated by the SSL-D method [23] and the proposed SSL-BD method, respectively. In Figure 14, the reflectivity predicted by our proposed SSL-BD method has a similar trend to the well reflectivity. However, the resolution of the predicted reflectivity is much lower than the well reflectivity, which needs to be further improved in the future. The predicted reflectivity profiles of the SSL-D method [23] and the SSL-BD method are presented in Figure 15a and Figure 15b, respectively. The results demonstrate that both methods can enhance the resolution of seismic profiles, but the proposed SSL-BD method reveals detailed information between layers more clearly, as indicated by the arrows and ellipses. Figure 15c shows the reconstructed seismic profile of the SSL-D method [23]. It is the convolution result of the reflectivity profile (as shown in Figure 15a) and the wavelet extracted from seismic data (as shown in Figure 12a). Figure 15d shows the reconstructed seismic profile of the proposed SSL-BD mothod. It comes from the convolution of the predicted reflectivity profile (as shown in Figure 15b) and the estimated wavelet (as shown in Figure 12b). The reconstructed seismic profiles shown in Figure 15c,d are visually indistinguishable from the original seismic data (as shown in Figure 11), indicating that both methods yield physically consistent predictions. Figure 15e,f provide residuals between the original seismic profile and the reconstructed seismic profiles.

3.3. Field Experiment II

The second field dataset comes from central Alberta, Canada, and is a 3D volume. These seismic data have 16 × 128 traces and each trace has 320 time sampling points, with a sampling interval of 2 ms, as shown in Figure 16 [23].

Figure 17a shows the wavelet extracted from the iErsk3D seismic data [23]. Figure 17b shows the wavelet estimated by the proposed SSL-BD method. Figure 18 shows the comparisons of wavelet amplitude spectra and phase spectra between the wavelet extracted from the iErsk3D data and the estimated wavelet. Figure 19a,b display the reflectivity inversion results of the SSL-D method [23] and the proposed SSL-BD method, respectively. Consistent with the conclusions drawn from the 2D field example, both methods enhance the resolution of seismic data, but the proposed SSL-BD method yields a sparser result and thus reveals thin layers more clearly. Figure 20 shows the amplitude spectra of the original field seismic data (black curve), the amplitude spectra of the predicted reflectivity of the SSL-D method [23] (green curve), and the amplitude spectra of the predicted reflectivity of the proposed SSL-BD method (blue curve). Both methods are capable of extending the frequency bandwidth of the original data. However, the proposed SSL-BD method demonstrates two notable advantages; it not only recovers more high-frequency components, but also preserves the low-frequency components better.

4. Discussion

In this paper, we propose a self-supervised machine learning blind deconvolution method based on joint inversion of seismic wavelet and reflectivity. Through the estimation of the initial zero-phase wavelet and the sparsity-promoted process on predicted reflectivity, we mitigate the ill-posedness of the blind deconvolution under a more concise parameter selection. Though the SSL-BD method achieves impressive deconvolution performance, there are several aspects that can be further explored and improved in the future.

Firstly, for the joint inversion of reflectivity and wavelet, we utilize a residual neural network and a fully connected convolutional neural network. However, many efficient network structures and components have been proposed and achieve superior experimental results. Therefore, different network components can be explored to enhance the deconvolution performance. Secondly, the proposed SSL-BD method derives reflectivity and seismic wavelet solely from seismic data, and does not make full use of well-logging data. In this regard, adding stratigraphic information from well-logging data to our framework will be a fruitful topic. Thirdly, for each new seismic survey, the proposed SSL-BD method needs to retrain networks so as to obtain the corresponding reflectivity and wavelet. This method does not take full advantage of the generalization properties of neural networks. We consider taking strategies like transfer learning and pre-training into our framework to improve the computational efficiency in the future. Fourthly, due to the limited computational resources, the numerical examples presented in this paper are in small-scale or middle-scale. However, due to the scalability of neural networks, our method should be able to deal with large-scale seismic datasets with better computational resources. Fifthly, the wavelet is assumed to be stationary in this work. For future work, the proposed method has potential to be generalized to the nonstationary seismic blind deconvolution by combining with a quality-factor-based nonstationary convolution model or segmentation-based nonstationary convolution model.

5. Conclusions

We present a novel blind deconvolution method based on self-supervised machine learning in this paper. This method eliminates the necessity for generating labeled pairs. Moreover, to overcome the dependence on known wavelets in the existing self-supervised machine learning methods, we leverage a residual neural network and a fully connected convolutional neural network to invert reflectivity and wavelet simultaneously. To mitigate the ill-posedness of blind deconvolution, we employs two strategies. The first one is estimating an initial zero-phase wavelet by smoothing the spectrum of averaged seismic data. The second one is applying a relative sparsity process on the predicted reflectivity during the joint inversion. Numerical experiments on synthetic data and field data show that the proposed method can obtain reflectivity inversion results with higher resolution than the existing self-supervised machine learning method and eliminates the need for a given wavelet. Our proposed method has potential to be improved by using more advanced neural networks and is promising for generalizing to nonstationary blind deconvolution.

Author Contributions

Conceptualization, X.Y., W.X. and B.W.; methodology, X.Y. and W.X.; formal analysis, X.Y., W.X. and B.W.; investigation, X.Y.; data curation, X.Y., W.X., Z.Y. and B.W.; writing—original draft preparation, X.Y.; writing—review and editing, W.X., X.Y., B.W. and Z.Y.; visualization, W.X. and B.W.; supervision, W.X. and B.W.; funding acquisition, B.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The part of the Marmousi2 model comes from https://wiki.seg.org/wiki/Dictionary:Marmousi_model/zh (accessed on 19 August 2023). The iErsk3D field data are obtained from https://doi.org/10.5281/zenodo.7233751 (accessed on 14 September 2023). The CNPC field data are unavailable upon request due to restrictions, e.g., privacy and ethical restrictions.

Acknowledgments

We thank two anonymous reviewers for their constructive comments on this paper. We also would like to thank Xintao Chai et al. for the open-source code of their geophysics-steered self-supervised learning for deconvolution method, and thank Zhifang Yang et al. for the CNPC field data.

Conflicts of Interest

Author Wenhao Xu was employed by the company BGP INC. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

SSL-D	self-supervised learning method for deconvolution
SSL-BD	self-supervised learning method for blind deconvolution
MAE	mean absolute error
MSE	mean squared error
MARE	mean absolute relative error
MSRE	mean square relative error
SSIM	structural similarity
PSNR	peak signal-to-noise ratio

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Framework of the proposed SSL-BD algorithm.

Figure 2. Illustration for the estimation of initial zero-phase wavelet. (a) Averaged seismic data waveform. (b) Averaged seismic data amplitude spectrum. (c) Smoothed seismic data amplitude spectrum. (d) Estimated initial zero-phase wavelet.

Figure 3. Reflectivity inversion network structure of 3D seismic data.

Figure 4. Wavelet phase inversion network structure.

Figure 5. [Forumla omitted. See PDF.] with number of epochs for training process of the synthetic data.

Figure 6. Estimated wavelet and true wavelet of the synthetic data.

Figure 7. Two-dimensional synthetic test of part of Marmousi2 model. (a) Amplitude spectra of true wavelet. (b) Amplitude spectra of wavelet estimated by the proposed SSL-BD method. (c) Phase spectra of true wavelet. (d) Phase spectra of wavelet estimated by the proposed SSL-BD method.

Figure 8. Two-dimensional synthetic test of part of Marmousi2 model. (a) Input seismic data. (b) True reflectivity profile. (c) Predicted reflectivity profile of SSL-D [23]. (d) Predicted reflectivity profile of the proposed SSL-BD method.

Figure 9. Detailed comparison of a single trace. (a) True reflectivity sequence. (b) Predicted reflectivity sequence.

Figure 10. Reflectivity inversion results with regard to the different noise levels. (a–d) Reflectivity inversion results of SSL-D [23] when SNR equals 5, 10, 15, and 20 dB. (e–h) Reflectivity inversion results by SSL-BD when SNR equals 5, 10, 15, and 20 dB.

Figure 11. Two-dimensional field test of CNPC data.

Figure 12. Two-dimensional field test of CNPC data. (a) Wavelet extracted from the CNPC seismic data. (b) Wavelet estimated by the proposed SSL-BD method.

Figure 13. Two-dimensional field test of CNPC data. (a) Amplitude spectra of wavelet extracted from the CNPC data. (b) Amplitude spectra of wavelet estimated by the proposed SSL-BD method. (c) Phase spectra of wavelet extracted from the CNPC data. (d) Phase spectra of wavelet estimated by the proposed SSL-BD method.

Figure 14. Detailed comparison of CNPC well log data. (a) Impedance curve from the well. (b) Reflectivity calculated using the impedance curve. (c) Reflectivity calculated by the SSL-D method [23]. (d) Reflectivity calculated by the proposed SSL-BD method.

Figure 15. Two-dimensional field test of CNPC data. (a,b) Predicted reflectivity profiles of the SSL-D method [23] and the proposed SSL-BD method. The vertical red line represents a well and the black curve is the corresponding impedance curve. (c,d) Reconstructed seismic profiles of the SSL-D method and the SSL-BD method. (e,f) Residuals between the original seismic profile (Figure 11) and the reconstructed seismic profiles of the SSL-D method and the SSL-BD method.

Figure 16. Three-dimensional field test of iErsk3D data.

Figure 17. Three-dimensional field test of iErsk3D data. (a) Wavelet extracted from the iErsk3D seismic data. (b) Wavelet estimated by the proposed SSL-BD method.

Figure 18. Three-dimensional field test of iErsk3D data. (a) Amplitude spectra of wavelet extracted from the iErsk3D seismic data. (b) Amplitude spectra of wavelet estimated by the proposed SSL-BD method. (c) Phase spectra of wavelet extracted from the iErsk3D seismic data. (d) Phase spectra of wavelet estimated by the proposed SSL-BD method.

Figure 19. Three-dimensional field test of iErsk3D data. (a,b) Predicted reflectivity profiles of the SSL-D method [23] and the proposed SSL-BD method. (c,d) Reconstructed seismic profiles of the SSL-D method and the proposed SSL-BD method. (e,f) Residuals between original seismic profile (Figure 13) and the reconstructed seismic profiles of the SSL-D method and the proposed SSL-BD method.

Figure 20. Three-dimensional field test of iErsk3D data. Amplitude spectra of the original field seismic data, the predicted reflectivity of the SSL-D method [23], and the predicted reflectivity of the proposed SSL-BD method.

Table 1

Runtimes of the SSL-D and SSL-BD with regard to three datasets.

	Part of Marmousi2	CNPC Data	iErsk3D Data
SSL-D	589.2 s	1856.1 s	4424.3 s
SSL-BD	1965.7 s	5558.3 s	13,590.3 s

Table 2

Evaluation metrics for comparison of different results.

Evaluation Metric	Formula
Mean absolute error (MAE)	$MAE=\frac{1}{N_{tr}}\frac{1}{L}{\sum }_{j=1}^{N_{tr}}{\sum }_{i=1}^L|{\mathrm{Actual}}_{i,j}-{\mathrm{Predicted}}_{i,j}|$
Mean squared error (MSE)	$MSE=\frac{1}{N_{tr}}\frac{1}{L}{\sum }_{j=1}^{N_{tr}}{\sum }_{i=1}^L{({\mathrm{Actual}}_{i,j}-{\mathrm{Predicted}}_{i,j})}^2$
Mean absolute relative error (MARE)	$MARE=\frac{1}{N_{tr}}\frac{1}{L}{\sum }_{j=1}^{N_{tr}}{\sum }_{i=1}^L\left|\frac{{\mathrm{Actual}}_{i,j}-{\mathrm{Predicted}}_{i,j}}{{\mathrm{Actual}}_{i,j}}\right|$
Mean square relative error (MSRE)	$MSRE=\frac{1}{N_{tr}}\frac{1}{L}{\sum }_{j=1}^{N_{tr}}{\sum }_{i=1}^L{(\frac{{\mathrm{Actual}}_{i,j}-{\mathrm{Predicted}}_{i,j}}{{\mathrm{Actual}}_{i,j}})}^2$
Structure similarity (SSIM)	$SSIM=\frac{(2·{\mathrm{Actual}}_{\mathrm{avg}}·{\mathrm{Predict}}_{\mathrm{avg}}+c_1)(2·{\mathrm{Covariance}}_{\mathrm{actual},\mathrm{predict}}+c_2)}{({\mathrm{Actual}}_{\mathrm{avg}}^2+{\mathrm{Predict}}_{\mathrm{avg}}^2+c_1)({\mathrm{Actual}}_{\mathrm{std}}^2+{\mathrm{Predict}}_{\mathrm{std}}^2+c_2)}$
Peak signal-to-noise ratio (PSNR)	$PSNR=10lg(\frac{{\mathrm{Actual}}_{max}^2}{\left.\frac{1}{N_{tr}}\frac{1}{L}{\sum }_{j=1}^{N_t}{\sum }_{i=1}^L\right|{\mathrm{Actual}}_{i,j}-{\left.{\mathrm{Predict}}_{i,j}\right|}^2})$

Table 3

Metric values of the predicted reflectivity profiles of 2D synthetic test.

	MAE (${10}^{-5}$) ↓	MSE (${10}^{-5}$) ↓	MARE↓	MSRE ↓	SSIM ↑	PSNR ↑
SSL-D	1.078175	1.335004	0.641659	0.434459	0.741858	26.806973
SSL-BD	0.599912	0.449987	0.258171	0.172519	0.995130	31.529794

↓ means the lower the better and ↑ means the higher the better. The best values are in bold.

Table 4

Comparison of reflectivity inversion results with regard to different noise levels.

	MAE (${10}^{-5}$) ↓		MSE (${10}^{-5}$) ↓		MARE ↓
SNR	SSL-D	SSL-BD	SSL-D	SSL-BD	SSL-D	SSL-BD
5 dB	1.165915	1.146042	1.947228	1.912455	0.949872	0.846458
10 dB	1.160701	1.081864	1.587014	1.464261	0.812649	0.806687
15 dB	1.080636	1.036662	1.495140	1.423236	0.795928	0.793276
20 dB	1.102799	1.052075	1.460031	1.384961	0.788522	0.772822
	MSRE ↓		SSIM ↑		PSNR ↑
SNR	SSL-D	SSL-BD	SSL-D	SSL-BD	SSL-D	SSL-BD
5 dB	1.319454	0.732748	0.659032	0.637587	25.688322	25.245887
10 dB	0.667327	0.656488	0.684185	0.703965	26.055993	26.405615
15 dB	0.637562	0.633944	0.700601	0.712404	26.314981	26.529031
20 dB	0.626091	0.602983	0.708102	0.722129	26.418179	26.647425

References

1. Gholami, A.; Sacchi, M.D. A fast and automatic sparse deconvolution in the presence of outliers. IEEE Trans. Geosci. Remote Sens.; 2012; 50, pp. 4105-4116. [DOI: https://dx.doi.org/10.1109/TGRS.2012.2189777]

2. Kazemi, N.; Bongajum, E.; Sacchi, M.D. Surface-consistent sparse multichannel blind deconvolution of seismic signals. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 3200-3207. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2513417]

3. Wang, L.; Gao, J.; Zhao, Q.; Xu, Z.; Jiang, X. Relative acoustic impedance inversion via toeplitz-sparse matrix factorization. Chin. J. Geophys.; 2017; 60, pp. 639-654.

4. Zhao, X.; Wang, S.; Yuan, S.; Cheng, L.; Cai, Y. Multichannel sparse spike deconvolution based on dynamic time warping. Acta Geophys.; 2021; 69, pp. 783-793. [DOI: https://dx.doi.org/10.1007/s11600-021-00571-z]

5. Peacock, K.I.; Treitel, S. Predictive deconvolution-theory and practice. Geophysics; 1969; 34, pp. 155-169. [DOI: https://dx.doi.org/10.1190/1.1440003]

6. Ulrych, T.J. Application of homomorphic deconvolution to seismology. Geophysics; 1971; 36, pp. 650-660. [DOI: https://dx.doi.org/10.1190/1.1440202]

7. Wiggins, R.A. Minimum entropy deconvolution. Geoexploration; 1978; 16, pp. 21-35. [DOI: https://dx.doi.org/10.1016/0016-7142(78)90005-4]

8. Wang, L.; Zhao, Q.; Gao, J.; Xu, Z.; Fehler, M.; Jiang, X. Seismic sparse-spike deconvolution via Toeplitz-sparse matrix factorization. Geophysics; 2016; 81, pp. V169-V182. [DOI: https://dx.doi.org/10.1190/geo2015-0151.1]

9. Wang, Y.; Zhang, G.; Li, H.; Yang, W.; Wang, W. The high-resolution seismic deconvolution method based on joint sparse representation using logging–seismic data. Geophys. Prospect.; 2022; 70, pp. 1313-1326. [DOI: https://dx.doi.org/10.1111/1365-2478.13232]

10. Xu, W.; Gao, J.; Wang, L.; Tian, Y.; Qiu, J.; Jiang, X.; Wang, D. Compact Smoothness and Relative Sparsity Algorithm for High-Resolution Wavelet and Reflectivity Inversion of Seismic Data. IEEE Trans. Geosci. Remote Sens.; 2022; 60, pp. 1-11. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3223280]

11. Asgarkhani, N.; Kazemi, F.; Jankowski, R. Machine learning-based prediction of residual drift and seismic risk assessment of steel moment-resisting frames considering soil-structure interaction. Comput. Struct.; 2023; 289, 107181. [DOI: https://dx.doi.org/10.1016/j.compstruc.2023.107181]

12. Li, G.; Zhou, X.; Chen, C.; Xu, L.; Zhou, F.; Shi, F.; Tang, J. Multitype Geomagnetic Noise Removal via an Improved U-Net Deep Learning Network. IEEE Trans. Geosci. Remote Sens.; 2023; 61, pp. 1-12. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3307422]

13. Tao, L.; Ren, H.; Gu, Z. Acoustic Impedance Inversion from Seismic Imaging Profiles Using Self Attention U-Net. Remote Sens.; 2023; 15, 891. [DOI: https://dx.doi.org/10.3390/rs15040891]

14. Zheng, Q.; Wei, C.; Yan, X.; Ruan, H.; Wu, B. Seismic Elastic Parameter Inversion via a FCRN and GRU Hybrid Network with Multi-Task Learning. Appl. Sci.; 2023; 13, 519. [DOI: https://dx.doi.org/10.3390/app131810519]

15. Tang, Z.; Wu, B.; Wu, W.; Ma, D. Fault Detection via 2.5D Transformer U-Net with Seismic Data Pre-Processing. Remote Sens.; 2023; 15, 1039. [DOI: https://dx.doi.org/10.3390/rs15041039]

16. Wu, X.; Liang, L.; Shi, Y.; Fomel, S. FaultSeg3D: Using synthetic data sets to train an end-to-end convolutional neural network for 3D seismic fault segmentation. Geophysics; 2019; 84, pp. IM35-IM45. [DOI: https://dx.doi.org/10.1190/geo2018-0646.1]

17. Chen, H.; Gao, J.; Gao, Z.; Chen, D.; Yang, T. A sequential iterative deep learning seismic blind high-resolution inversion. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens.; 2021; 14, pp. 7817-7829. [DOI: https://dx.doi.org/10.1109/JSTARS.2021.3100502]

18. Gao, Z.; Hu, S.; Li, C.; Chen, H.; Jiang, X.; Pan, Z.; Gao, J.; Xu, Z. A deep-learning-based generalized convolutional model for seismic data and its application in seismic deconvolution. IEEE Trans. Geosci. Remote Sens.; 2021; 60, 4503117. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3076991]

19. Gao, Y.; Zhao, D.; Li, T.; Li, G.; Guo, S. Deep Learning Vertical Resolution Enhancement Considering Features of Seismic Data. IEEE Trans. Geosci. Remote Sens.; 2023; 61, 5900913. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3234617]

20. Min, F.; Tang, J.; Pan, S.; Song, G. CAUC: Combining Channel Attention U-Net and Convolution for Seismic Data Resolution Improvement. IEEE Geosci. Remote Sens. Lett.; 2023; 20, 3003105. [DOI: https://dx.doi.org/10.1109/LGRS.2023.3322263]

21. Wang, Y.; Xu, J.; Zhao, Z.; Gao, Y.; Zhang, H. Structurally-Constrained Unsupervised Deep Learning for Seismic High-Resolution Reconstruction. IEEE Trans. Geosci. Remote Sens.; 2024; 62, 5901115. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3340888]

22. Sun, J.; Yang, J.; Li, Z.; Huang, J.; Luo, X.; Xu, J. Intelligent ava inversion using a convolution neural network trained with pseudo-well datasets. Surv. Geophys.; 2023; 44, pp. 1075-1105. [DOI: https://dx.doi.org/10.1007/s10712-022-09766-5]

23. Chai, X.; Yang, T.; Gu, H.; Tang, G.; Cao, W.; Wang, Y. Geophysics-steered self-supervised learning for deconvolution. Geophys. J. Int.; 2023; 234, pp. 40-55. [DOI: https://dx.doi.org/10.1093/gji/ggad015]

24. Hildebrand, F.B. Introduction to Numerical Analysis; Courier Corporation: New York, NY, USA, 1987.

25. Cole, A.J.; Davie, A.J.T. Local smoothing by polynomials in n dimensions. Comput. J.; 1969; 12, pp. 72-76. [DOI: https://dx.doi.org/10.1093/comjnl/12.1.72]

26. Asgarkhani, N.; Kazemi, F.; Jakubczyk-Gałczyńska, A.; Mohebi, B.; Jankowski, R. Seismic response and performance prediction of steel buckling-restrained braced frames using machine-learning methods. Eng. Appl. Artif. Intell.; 2024; 128, 107388. [DOI: https://dx.doi.org/10.1016/j.engappai.2023.107388]