Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 427153, 11 pages

http://dx.doi.org/10.1155/2015/427153

## A Seismic Blind Deconvolution Algorithm Based on Bayesian Compressive Sensing

^{1}School of Electrical Engineering & Automation, Tianjin University, Tianjin 300072, China^{2}Department of Disaster Prevention Equipment, Institute of Disaster Prevention, Beijing 101601, China

Received 19 March 2015; Accepted 20 April 2015

Academic Editor: Wanquan Liu

Copyright © 2015 Yanqin Li and Guoshan Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Compressive sensing in seismic signal processing is a construction of the unknown reflectivity sequence from the incoherent measurements of the seismic records. Blind seismic deconvolution is the recovery of reflectivity sequence from the seismic records, when the seismic wavelet is unknown. In this paper, a seismic blind deconvolution algorithm based on Bayesian compressive sensing is proposed. The proposed algorithm combines compressive sensing and blind seismic deconvolution to get the reflectivity sequence and the unknown seismic wavelet through the compressive sensing measurements of the seismic records. Hierarchical Bayesian model and optimization method are used to estimate the unknown reflectivity sequence, the seismic wavelet, and the unknown parameters (hyperparameters). The estimated result by the proposed algorithm shows the better agreement with the real value on both simulation and field-data experiments.

#### 1. Introduction

Digital seismic signal receiving devices have been developed for many years. In traditional seismic signal receiving system, the accurate signal is obtained by Nyquist sampling theorem which claims that the seismic signal sampling rate should be at least twice the seismic signal’s highest frequency. Measurement reduction, such as compressive sensing (CS), can make the seismic record devices be simpler, smaller, and cheaper. CS is a signal processing method which can reconstruct the original signal from a small number of random incoherent liner measurements [1, 2]. CS uses the property of the signal’s sparseness or compressibility in some domain, which achieves the entire signal from relatively few measurements [3]. CS has become very popular in recent years because it can reduce the number of data transmissions and give potential aid in many practical applications [4, 5].

The CS problem is to compress data from original data while the inverse CS problem is to reconstruct original data from the compressed data. There exist lots of effective methods to solve the CS problem. Some of these methods are based on energy limit [6, 7] and the others are based on the Bayesian framework [8, 9]. Bayesian compressive sensing gives a Bayesian framework for solving the inverse problem of compressive sensing (CS). The basic Bayesian CS algorithm exploits the relevance vector machine and then it marginalizes the noise variance with improved robustness [10, 11].

Deconvolution, as a seismic data processing technique, has been widely used to yield reflectivity sequence from the seismic records [12–15]. The Bernoulli-Gaussian model introduced in [16], as the reflectivity sequence model, is popular. It is very intuitive and its statistical properties can easily apply in Bayesian methods [17–19] such as Gibbs sampling and MCMC (Markov Chain Monte Carlo) techniques. In statistical models, an expectation maximization (EM) algorithm is an iterative method for finding maximum likelihood (ML) or maximum a posteriori (MAP) estimates of parameters, where the model depends on unobserved latent variables [20]. So, EM algorithm is usually used to estimate the unknown parameters (the wavelet, parameters of the BG model and additive noise variance) in blind seismic deconvolution. The EM algorithm iteration alternates between performing an expectation (E) step and a maximization (M) step. E step creates a function for the expectation of the log-likelihood while M step computes parameters maximizing the expected log-likelihood found on the E step. These parameter-estimates are then used to determine the distribution of the latent variables in the next E step.

The Bayesian method is model-based method which is applied to estimate unknown parameters using the prior information. So it can handle observation noise and missing data problems [21]. Moreover, sparsity, as the important characteristic of a signal, has been exploited in many application areas, such as seismic deconvolution; see [22–24]. It could be an alternative method for the deconvolution of sparse signals that model a sparse signal with a Gaussian random variable whose variance is an inverse-Gamma random variable [25–29]. It exploits variational Bayesian techniques for a closed-form approximate posterior inference [30].

In this paper, the reflectivity sequence and the seismic wavelet are recovered through CS measurement of the seismic records. Bayesian framework is used in CS problem in seismic signal processing. Blind seismic deconvolution and optimization method are combined to estimate the unknown reflectivity sequence, the seismic wavelet, and the hyperparameters. The organization of this paper is as follows. In Section 2, the Bayesian model with CS measurement of the blind seismic restoration problem is given. In Section 3, the observation model and the prior model on the signals coefficients are presented. In Section 4, the Bayesian inference is presented. Optimization method is used to get the solution of the unknown reflectivity sequence and the parameters. In Section 5, comparison between ML and Bayesian deconvolution with CS measurement is given based on the both simulation and field data experiments. In the last section, we draw the conclusions.

#### 2. Problem Formulation

Seismic record is the convolution with the seismic wavelet and the reflectivity sequence. A standard seismic record image model is given in a matrix-vector form by where denotes a seismic observation vector, denotes an unknown seismic reflectivity sequence vector, and is an additive white Gaussian noise vector, respectively. , which represents seismic wavelet, is a block Toeplitz matrix obtained by a block-circulant matrix. The seismic wavelet and reflectivity vector are both unknown, so the deconvolution is a blind problem.

Considering the sparsity of the seismic reflectivity vector, is compressible by a linear wavelet basis ; that is, , where is a sparse signal. Reference [18] proposed sparse schemes that capture the convolution structure of pulse streams by defining (2) which resorts to a CS matrix to complete measurements.

Equation (1) can be rewritten in the following form: where . Based on the sparsity of , the norm of is used in the inverse problem. The estimation of the original signal can be obtained by solving the following optimization problem [31]: where represents the norm.

#### 3. Bayesian Modeling

The prior distribution of unknown signal is assumed to be . The observation is the random process with distribution , where denotes the unknown noise variance. These seismic record and reflectivity depend on the unknown hyperparameters and . Following the Bayesian framework of the CS inverse problem, the following joint posterior distribution can be expressed by

In Bayesian estimation, all unknown parameters are assumed to be random variables with specific distribution. Our goal is to estimate the parameters , , , and from . The parameter can be estimated by the Stochastic expectation maximization algorithm of Rosec et al. [32]. The following gives how to get the other parameters.

##### 3.1. Seismic Observation Modeling

By assuming that the seismic data noise in (2) is Gaussian distribution with zero mean and variance equal to , the conditional distribution of the seismic record is

with a Gamma prior distribution on as the following:

The Gamma distribution for the hyperpriors is defined as where represents hyperparameter, is the scale parameter, and is the shape parameter. These parameters are also assumed to be random variables. We will show how they are calculated in the following subsection. The Gamma distribution has the following mean and variance:

##### 3.2. Signal Modeling

In [33], as the reflectivity sequence also has sparsity, can use the following prior model:where . Then, hyperprior of the hyperparameter is as the following:

Based on (9) and (10), we have the following equation:

At last, hyperparameter is modeled as the following Gamma hyperprior:

The modeling has three stage forms. The first stage is the prior distribution about in (9). The second stage is the prior distribution about in (10). Finally, (12) is used to calculate . So we can first sample a from the prior distribution on distribution to get and then sample a distribution times to get , , and at last sample a distribution to get .

Based on the above, the joint distribution can be given as the following:where , , , , and are defined in (5), (6), (9), (10), and (12), respectively.

The seismic record variables of the proposed model are . The parameter is while the hyperparameters are and . The parameter of the hyperparameters model is . The dependencies among the random variables that define the proposed Bayesian compressive sensing model are shown in Figure 1.