Yan Zhang , Weijian Ren,Guowei Tang and Can Zhao

(1.School of Computer and Information Technology, Northeast Petroleum University, Daqing 163318, Heilongjiang, China;2.School of Electrical Engineering and Information, Northeast Petroleum University, Daqing 163318,Heilongjiang, China)

Abstract: Recovery of under-sampled seismic data is a critical problem, in oil and gas exploration, therefore recovery algorithms with iterative shrinkage based on compressed sensing have been recently proposed. However most of these algorithms usually adopt a soft shrinkage function, which assumes that all of the sparse coefficients are independent of each other in curvelet or other domains, little attention has so far been devoted to the inter-dependencies of coefficients. In this paper, the dependencies of parent-child curvelet coefficients of seismic data are exploited by Bayesian estimation, moreover the new seismic data recovery algorithm via curvelet-based bivariate shrinkage function is proposed. First the respective parent-child curvelet coefficients joint distribution models of fully-sampled seismic data and noise signal caused by missing traces are established, then the bivariate shrinkage function according to the Bayesian maximum posterior probability estimation is obtained, finally the Landweber iterative shrinkage algorithm is used in the recovery process. When compared with existing recovery algorithms, it is proved that the proposed algorithm can obtain higher PSNR performance, and maintains the texture details better in events of seismic data recovery .

Keywords: Seismic data recovery; compressed sensing; iterative shrinkage; bivariate shrinkage function

1 Introduction

Under the influences of both the geological environment and acquisition conditions, seismic data is prone to dead traces, insufficient spatial sampling, and noise introduction, while demands for the integrity and regularity of seismic data are high in the subsequent processing and interpretation, hence the recovery of undersampled data is an extremely critical problem. The traditional seismic data recovery method[1-2]is restricted by the Nyquist sampling theory, which rules that the sampling frequency is at least two times of the highest frequency of the signal, otherwise it is very easy to result in aliased, consequently affecting the quality of recovery. Recently, the theory of compressed sensing as proposed by Candès[3]and Donoho[4]claims that under certain conditions, signals to be sampled at sub-Nyquist rates via linear projection onto a random basis while still could be accurately recovered[5-7], which provides a new idea for the problem of seismic data recovery. Curvelet transformation[8]with the significant advantages of directional sensitivity and anisotropy, can effectively represent multiple directional events, further making use of sparse seismic data, which has been used extensively in seismic data processing[9-10]. Herrmann and Hennenfent in Ref.[11] applied compressed sensing to the problem of seismic data recovery in the curvelet domain for the first time, further the curvelet-based recovery by the sparsity-promoting inversion[12]method was proposed the solving the sparse optimization problem. Shahidi et al.[13]used the cooling iterative threshold method to solve the L1 norm minimization problem in seismic data recovery, with the prior knowledge that the seismic data could be sparsely represented in the curvelet domain. The above solutions of the recovery problems depend critically on the seismic data’s sparsity in curvelet transformed domain, as well as the simple and effective soft threshold in the iterative shrinkage process based on compressed sensing. But the soft threshold method inherently assumes that all the curvelet coefficients are completely independent, hence the classical soft shrinkage function is used to determine the significances of all curvelet coefficients in each iteration, and however the dependencies between curvelet coefficients are entirely ignored. Naghizadeh[14]also reported a curvelet transform interpolation scheme that the curvelet coefficients were divided into two groups of scales: alias-free (coarse) scalse and alias-contaminated (high) scales, then the alias-free curvelet coefficients were upscaled to estimate a mask function that was used to constrain the inversion of the alias-contaminated scale coefficients, but the strong dependencies between a coefficient and its parent were not captured as well. Bayesian estimation techniques were widely used in the field of data processing[15-17], moreover Sendur[18]proposed the bivariate shrinkage functions using the bayesian maximum posterior probability estimation to denoise natural image in wavelet domain, which yielded performance superior to soft shrinkage function, in that it exploited the statistical dependencies of neighborhood and parent-child.

Amongst all these methods, it was observed that the dependencies between curvelet coefficients should be taken into account to calculate the shrinkage function instead of the simple soft shrinkage function, when the iterative shrinkage process was used to recover the seismic data based on compressed sensing. Inspired by Ref.[18], we focus on the parent-child curvelet coefficients dependencies of seismic data with the Bayesian estimation techniques in order to address this issue. Our contribution lies in that we model the joint parent-child distribution of curvelet coefficients, furthermore obtaining the bivariate shrinkage function based on bayesian maximum posterior probability estimation, in order to predict the true coefficients of observed seismic data (undersampled), and then the distribution characteristics of the finest-scale and neighbor curvelet coefficients of seismic data are used to determine the parameters, finally the Landweber iterative shrinkage algorithm is used in the recovery process. In experimental simulations, we find that the proposed algorithm usually matches or exceeds PSNR performance of the soft shrinkage function, and maintains the texture details better in the events of the recovery of seismic data.

2 Principle of Seismic Data Recovery Based on Compressed Sensing

Setf∈RNt×Nrthe fully-sampled seismic data, whereNris the number of sampling coordinates, andNtis the number of time coordinates, thenfcan be vectorized asx∈RNtNr. Similarly the observed seismic data is set toq∈RNt×Mr,Mr

y=Rx

(1)

whereRis a randomly sampling matrix, which is used to extractyfromx. Though the Eq. (1) is a set of underdetermined equations, in whichxcan not be solved according toy. But as we all known that the coefficients of seismic data in the curvelet domain is sparse, letαas the curvelet coefficients of the seismic data,Cas the curvelet transform, andCTthe inverse curvelet transform, then the Eq.(1) can be written as Eq.(2).

y=Rx=RCTα=ACSα

(2)

According to the theory of compressed sensing, if the sensing matrixACSmeets the Isometric Property Restricted (RIP)[19], the fully-sampled seismic data could be recovered by the1norm optimization problem as Eq.(3).

min‖α‖1s.t.ACSα=RCTα=y

(3)

Asαis sparse, it can be solved by the sparsity-promoting inversion algorithm of the formula (4).

(4)

α(k+1)=T(α(k)+RT(y-Rα(k)))

(5)

(6)

3 Multi-scale Analysis of Seismic Data in Curvelet Domain

A curvelet is a geometric and directional wavelet, which has anisotropic needle-shape elements, and makes an optimal sparse representation of signals with continuities along smooth curves. Where the curvelet basis and signal edges are in the same direction, it will produce large coefficiencies; otherwise the random noise usually corresponds to small coefficients in the cruvelet domain. Just for this reason the curvelet is more appropriate for the two-dimensional signals with curves and edges, and it also will achieve higher approximation accuracy and better sparse representation. The seismic data is gathered in the process of artificial seismic exploration, thus they are of a huge dynamic range with a wealth of low-intermediate frequency information, in addition problems such as an insufficient sampling rate and low accuracy often occur due to the effect of technical, and/or geological conditions, and the effects of human action. However the wave fronts in the seismic data are usually composed of curved lines, which are very similar to the curve shape primitives of curvelet transform as shown in Fig.1. Therefore the curvelet transform has extensively been used to approximate the wave front curve in seismic data, as aforementioned, the stronger events directivity are, the greater curvelet coefficients, and vice versa. The curvelet coefficients of Fig.1 as shown in Fig.2(by 5 scales decomposition). The white subbands are the isolated strips of curvelet coefficients, and they are not the real coefficients, just enough to satisfy the needs of the display. The black subbands are the curvelet coefficients regions, in which the brighter the pixels are in these regions, the greater the value of the coefficients, and conversely the smaller. The coefficients in the red boxes represents the same direction of different scales, where the scales vary from coarse to fine as the regions move from the center to outwards. Take the regions surround by red boxes as an example, according to the multi-scale feature of curvelet coefficients, the parent-child relationship of curvelet coefficients are established from above to below, as shown in Fig.3(For the purpose of convenient observation, the coefficients of all directions are mapped in the same size). The distribution of current coefficients are very similar to the parents, if the paternal coefficients are important, the current coefficients are likely to be important, so it can be used to estimate the importance of current coefficients.

Fig. 1 Fully-sampled Seismic data

Fig.2 Curvelet coefficients

Usually, the noise affected by random missing traces are corresponding to the relatively small coefficients in the curvelet domain, where the soft shrinkage function that main idea is to retain from all coefficients larger than the threshold, and further set all other coefficients to zero, is effectively promoting sparsity under the relatively high sampling rate, and it will gain a better recovery effect. However when the sampling rate is becomes lower, the curved feature of events in observed seismic data are destroyed, as a result, without the obvious directional information, there will be a large number of small curvelet coefficients. Furthermore, the true signal coefficients and noise coefficients can not be distinguished, if we determine the significance of coefficients only based on amplitude like the soft shrinkage function, consequently it will cause recovery problems of due to discontinuous events or serious blocking artifacts.

Fig.3 Parent-child relationship of curvelet coefficients

4 Multi-Scale Analysis of Seismic Data in Curvelet Domain

4.1 Curvelet Bivariate Shrinkage Function

There are strong correlations among the curvelet coefficients which are between adjacent scales, so it is possible to predict the true value of the observed seismic data, according to the respective joint parent-child curvelet coefficients distribution models of fully-sampled seismic data and noise signal, coupled with the Bayesian maximum a posteriori estimator.

Setγis the curvelet coefficient of observed seismic data,wis the fully-sampled seismic data curvelet coefficient, andεis the noise signal curvelet coefficient, so there is:γ=w+ε.

According to Bayesian theory, under the condition thatγis known, the maximum posteriori probability is estimated as the formula (7) .

(7)

The formula (8) is obtained by using Bayesian rule:

(8)

It is equivalent to Eq.(9),

(9)

wherepwis the probability distribution function of curvelet coefficients of full-sampled seismic data,pεis the probability distribution function of noise signal curvelet coefficients. Letγ2represent the parent ofγ1. (γ2is the curvelet coefficient at the same direction asγ1, but at the adjacent coarser scale) , Similarly letw2represent the parent ofw1, andε2represent the parent ofε1, so one can getw=(w1,w2),γ=(γ1,γ2),ε=(ε1,ε2).

In this paper, the respective distributions in various scales of curvelet coefficients are obtained according to a lot of experiments, as a result the former is shown as Fig.4, which is approximate to the Gaussian distribution, while the latter is shown as Fig.5, which is approximate to the Laplasse distribution.

Fig.4 Distribution of noise

We model the noise caused by missing traces in each scale curvelet domain by a zero mean Gaussian distribution, so the joint parent-child coefficients probability distribution function can be defined as Eq.(10).

Fig. 5 Distribution of fully-sampled Seismic data

(10)

whereσεis the mean variance of curvelet coefficients of noise signal.

In that the curvelet coefficients of the fully-sampled seismic data are more approximate to the Laplasse distribution, the joint probability distribution function of the parent-child coefficients can be defined as the formula (11):

(11)

whereσis the mean variance of curvelet coefficients of the fully-sampled seismic data.

Bring formula (10), (11) into (9), obtained formula (12):

(12)

Ifpw(w) is convex and differentiable, after finishing one can get the formula (13):

(13)

4.2 Bivariate Shrinkage Function Based Iterative Recovery Algorithm

In the recovery process based on compressed sensing, we use the Landweber iterative algorithm with the iterative shrinkage strategy, then the bivariate shrinkage function is used to promote sparsity in each iteration such as Eq.(14):

(14)

The proposed algorithm of seismic data recovery based on compressed sensing is described as follows:

Algorithm1The seismic data recovery algorithm with curvelet bivariate shrinkage function based on compressed sensing.

Algorithm steps:

Input: Given the iteration stop parametersτ, the maximum number of iterationsK, the seismic data observation matrixR, observed seismic datay.

Output: Recovery of seismic datax.

2) while (k≤K)

α(k)=C(x(k))

fori=2 tos

xk+1=CT(α(k+1))

if (‖x(k+1)-x(k)‖2≤τ)

break;

end if

end while

5 Experimental Results and Analysis

In order to validate the feasibility and effectiveness of the proposed algorithm, we compare the proposed algorithm with the existing algorithms by using synthetic seismic data and the Marmousi model seismic data as examples. Our approach was that first we got both of the fully-sampled seismic data sets, and then they are extracted by randomly sampling matrix to generate the observed seismic data (with random missing traces), finally the proposed algorithm is used to recover the observed seismic data. The objective measurement of the recovery data is based on the peak signal to noise ratio (PSNR), as shown in the formula (15):

(15)

wherefis the vector of the fully-sampled seismic data, MSE is the mean square error offand the vector of the recovered seismic data.

5.1 Experimental Analysis of Synthetic Seismic Data

The partial single shot seismic synthetic records is shown in Fig.6(a) and the observed data under 70% sampling rate is shown in Fig.6(b). The Fig.6(c) shows the recovery effect of soft shrinkage function iterative algorithm, and the Fig.6(d) shows the recovery effect of proposed algorithm, as we can see the events in elliptical selection area of Fig.6(c) are fuzzy and have an obvious block effect, while the corresponding region in Fig.6(d) is relatively continuous and smooth. Fig.6(e) illustrates the comparison of PSNR evolution curves with respect to an iteration number with an under 70% sampling rate, with the increasing number of iterations, the PSNR of two algorithms are gradually improved, but the proposed algorithm is first to achieve convergence, and the convergence value is higher than that of the soft shrinkage function algorithm. Fig.6(f) illustrates the PSNR comparison of the two algorithms under different sampling rates, and we can see that the proposed algorithm in different sampling rates gains higher PSNR than the soft shrinkage function algorithm.

Fig.6 The experiment with synthetic data

Fig.7 shows the effect of the soft threshold algorithm and the proposed algorithm in 60, and 120 iterations. It can be seen that under the same number of iterations, the proposed algorithm can maintain events better, and the curve feature is more smooth and continuous.

Fig.7 The experiment of different number of iterations

5.2 Experimental Analysis of the Marmousi Model

The Marmousi model is used to test the proposed algorithm further, the fully-sampled seismic data is shown in Fig.8(a), and the observed data with an under 60% sampling rate is shown in Fig.8(b). The Fig.8(c) shows the recovery effect of the soft shrinkage function iterative algorithm, and the Fig.8(d) shows the recovery effect of the proposed algorithm, as we can see that the events in elliptical selection area of Fig.8(c) have a more prominent block effect, while the corresponding region in Fig.8(d) is relatively smooth. Fig.8(e) illustrates the comparison of PSNR evolution curves with respect to an iteration number under with a 60% sampling rate, with the increase in the number of iterations, we can see that the PSNR value of the proposed algorithm is lower than that of the former at the beginning, but after 20 iterations, the convergence speed is apparently accelerated, eventually converging to the higher PSNR. Because the dual tree complex wavelet transformation (DDWT) and contourlet transformation have a good ability for sparse representation with texture and contour information, Fig.8(f) illustrates the PSNR comparison of four algorithms under different sampling rates, it shows that the reconstruction effect of curvelet transformation is higher than that of the DDWT and contourlet, and the reconstruction effect of the bivariate threshold algorithm being better than that of the soft threshold reconstruction algorithm in different sampling rates. The above experiments have proved that the proposed algorithm is more effective and adaptable.

Fig.8 The experiment with marmousi model data

Table 1 demonstrates the run time comparison of the proposed algorithm with the soft threshold algorithm, DDWT and contourlet in the different sampling rates, under the 200 iterations. In the same sampling rate, and the run time of the proposed algorithm is higher than that of the soft threshold algorithm for 20-50 s, because in each iteration of the threshold processing, the threshold is fixed in soft threshold algorithm, thus it just needs to compare the value of the current coefficients with the threshold, while in the proposed algorithm, it needs to not only to search for parent coefficients of current coefficients, but also generate the respective thresholds through the variance of neighborhood, and therefore it increases the computation requirements. As the efficiency of DDWT and contourlet transformation is higher than the curvelet transformation, the running time is relatively short.

Fig.9 shows the recovered effect of the soft threshold algorithm and the proposed algorithm in 50, 100 iterations. It can be seen that with the increasing of the iterative number, the subjective effect of the two algorithm is enhanced, and the quality of the two recovered images is not so far different in the smooth region, but the reconstruction effect of the proposed algorithm is better than that of the former in the complex region of the events.

Table 1 Comparison of run time under different sampling rate (s)

Fig.9 The experiment of different number of iterations

Fig.10 illustrates the comparison of the recovered results between the soft threshold algorithm and the proposed algorithm, the two algorithms show some fuzzy phenomena at sampling rates of 0.5, 0.7 and 0.8. When the sampling rate is 0.5, the effect of the reconstructions of the complex region of the events. This is because the unified threshold processing is adopted to all scales in the soft threshold algorithm, while the proposed algorithm employs the correlation of the curvelet coefficients, thus the latter subjective effect is improved. With an increasing of sampling rate, such as the sampling rates of 0.7, 0.8, the advantages of the proposed algorithm in this paper are more obvious.

Fig.10 Experiment of different sampling rates

6 Conclusions

In the idea of the compressed sensing recovery coupled with the iterative shrinkage method based on curvelet transform, the recovery algorithm using the curvelet bivariate shrinkage function is proposed. We first model the joint distribution of parent-child curvelet coefficients of fully-sampled seismic data by the Laplasse distribution, and the noise signal by zero mean Gauss distribution, and then the curvelet bivariate shrinkage function based on Bayesian maximum posteriori estimation is obtained. Furthermore, the distribution of finest-scale coefficients and the interscale dependencies are exploited to estimate the parameters of the curvelet bivariate shrinkage function. Finally, the Landweber iterative algorithm is used in the recovery process. The experiments on the synthetic and Marmousi model seismic data prove that the proposed algorithm can obtain higher recovery PSNR and has a better effect in maintaining the texture than the existing algorithms.