UBC Faculty Research and Publications

Curvelet-based primary-multiple separation from a Bayesian perspective Saab, Rayan; Wang, Deli; Yilmaz, Ozgur; Herrmann, Felix J. 2007

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-saab07seg.pdf [ 246.84kB ]
JSON: 52383-1.0107412.json
JSON-LD: 52383-1.0107412-ld.json
RDF/XML (Pretty): 52383-1.0107412-rdf.xml
RDF/JSON: 52383-1.0107412-rdf.json
Turtle: 52383-1.0107412-turtle.txt
N-Triples: 52383-1.0107412-rdf-ntriples.txt
Original Record: 52383-1.0107412-source.json
Full Text

Full Text

Curvelet-based primary-multiple separation from a Bayesian perspective ¨ ur Yılmaz, Math-UBC and Felix J. Herrmann, EOS-UBC Rayan Saab, ECE-UBC, Deli Wang, Jilin University. Ozg¨ SUMMARY In this abstract, we present a novel primary-multiple separation scheme which makes use of the sparsity of both primaries and multiples in a transform domain, such as the curvelet transform, to provide estimates of each. The proposed algorithm utilizes seismic data as well as the output of a preliminary step that provides (possibly) erroneous predictions of the multiples. The algorithm separates the signal components, i.e., the primaries and multiples, by solving an optimization problem that assumes noisy input data and can be derived from a Bayesian perspective. More precisely, the optimization problem can be arrived at via an assumption of a weighted Laplacian distribution for the primary and multiple coefficients in the transform domain and of white Gaussian noise contaminating both the seismic data and the preliminary prediction of the multiples, which both serve as input to the algorithm.  INTRODUCTION Predictive multiple suppression methods consist of two main steps: a prediction step, during which multiples are predicted from seismic data, and a primary-multiple separation step, during which the predicted multiples are ’matched’ with the true multiples in the data and subsequently removed. This second separation step, which we will call the estimation step, is crucial in practice. An incorrect separation will cause residual multiple energy in the result or may lead to a distortion of the primaries, or both. To reduce these adverse effects, a new transform-domain method was proposed by Herrmann et al. (2007), where primaries and multiples are separated in the curvelet domain rather than matched. This separation method is carried out on the basis of differences in the multiscale and multidirectional characteristics of the primaries and multiples. The sparsity in the curvelet domain of these two signal components is exploited by solving a weighted 1 -norm minimization problem that is designed to drive the two signal components apart. This separation is accomplished by using the (SRME) predicted primaries and multiples as weights. During the optimization, the two signal components are separated by enhancing sparseness (through weighted 1 -norms) in the transform domain subject to fitting the observed data as the sum of the separated components to within a user-defined tolerance level. Even though good results were obtained with this method, problems with the appropriate balancing of the weights were experienced. Lack of a good balance yielded a ’separation’ where all the signal energy ends up in one of the two signal components. To tackle this issue, we propose an alternative formulation for the primary-multiple separation where an additional constraint is imposed on the coefficients for the multiples. This formulation follows from a Bayesian approach where sparsity  for the two signal components as well as data misfits are taken into account. Our approach derives from techniques employed in the field of signal processing and geophysics where Bayesian approaches supplemented with sparsity-promoting Laplace/Cauchy probability density functions are frequently used. Our method differs from the earlier results in the literature by (i) using sparsity promoting curvelets (see, e.g., Candes et al., 2006); (ii) using inaccurate predictions for the to-be-separated signal components and (iii) providing a new iterative solution technique based on descent projections and soft thresholding. The outline of our abstract is as follows. First, we briefly discuss curvelets, our transform domain of choice, which has been shown to admit sparse coefficients for seismic data (Candes et al. (2006)). Next, Bayesian formulation of the primarymultiple separation problem is introduced, leading to the objective function to be minimized. Finally, a separation algorithm that provably converges to a minimizer of the objective function is given. Results of our algorithm applied to a synthetic shot with 3-D complexity are presented as well.  CURVELETS Curvelets are localized ’little plane-waves’ (see e.g. Hennenfent and Herrmann, 2006, and ancillary electronic material) that are oscillatory in one direction and smooth in the other direction(s). They are multiscale and multi-directional. Curvelets have an anisotropic shape – they obey the so-called parabolic scaling relationship, yielding a width ∝ length2 for the support of curvelets. This anisotropic scaling is optimal for detecting wavefronts and explains their high compression rates on seismic data and images (Candes et al., 2006; Hennenfent and Herrmann, 2006; Herrmann et al., 2007). Curvelets represent a specific tiling of the 2-D/3-D frequency plane into strictly localized multiscale and multi-angular wedges. Fig. 1(b)). Because the directional sampling increases everyother scale-doubling, curvelets become more anisotropic at finer scales. They become ’needle-like’ as illustrated in Fig. 1(a). The smoothness in the Fourier domain leads to a rapid decay in the physical domain. Their effective support is given by ellipsoids with a width ∝ 2 j/2 and length ∝ 2 j and an angle θ jl = 2πl2 j/2 with j the scale, l the angular index with a rounding off towards the nearest smaller integer. Curvelets are indexed by the multi-index µ := ( j, l, k) ∈ M with M the multi-index set that runs over the scales, j, angles, l, and positions k (see for details Candes et al., 2006; Ying et al., 2005). The number of wedges doubles every other scale-doubling and such that directional selectivity increases for finer scales (see Fig. 1(b)). Curvelets compose an arbitrary column vector f, with the reordered samples, according to f = CT Cf with C and CT , the forward/inverse discrete curvelet transform matrices (defined by the fast discrete curvelet transform, FDCT, with  Curvelet-based primary multiple separation wrapping Candes et al., 2006; Ying et al., 2005). CT denotes the transpose of C, which is equal to the pseudoinverse of C for our choise of discrete curvelet transform: the rows of C constitute a tight frame with frame bound 1 and with a moderate redundancy (a factor of roughly 8 for d = 2 and 24 for d = 3 with d the number of dimensions). Such tight frames (see e.g. Daubechies, 1992) are signal representations that preserve energy. Consequently, CT C = I with I the identity matrix. However, CCT = I .  BAYESIAN FORMULATION  0.50  0  100 0.25  k2  Samples  200 0  Here, we approach the primary-multiple separation problem from a probabilistic perspective, where we impose an a priori distribution for the unknown coefficients of the primaries and multiples. Estimates for the curvelet coefficients of the primaries and multiples are obtained by using a Laplaciantype probability density function for the curvelet coefficients and weights defined in terms of (SRME) predictions for the two signal components. The estimates are calculated within a Bayesian framework which leads to the minimization of a weighted sparsity-promoting functional. The choice of a Laplacian-type density function for the curvelet coefficients also serves as a sparsity promoting prior (see e.g. Zibulevsky and Pearlmutter, 2001; Li et al., 2004) and is consistent with the high compression rates that curvelets display for seismic data (Candes et al., 2006; Hennenfent and Herrmann, 2006; Herrmann et al., 2007).  300  More precisely, let -0.25  400  500 0  100  200  300  400  500  Samples  -0.50 -0.50  b = s1 + s2 + n  (1)  be the recorded total data set comprised of a noisy sum of primaries s1 and multiples s2 and let -0.25  0  0.25  k1  0.50  b2 = s2 + n2  (2)  (a)  be the predicted multiples, assumed to be erroneous versions of the true underlying multiples. Moreover, let the multiples and primaries admit (curvelet) transform coefficients satisfying  0.50  si = Axi , {i = 1, 2},  k2  0.25  where A is the inverse curvelet transform operator (defined by the fast discrete curvelet transform (FDCT by wrapping, see e.g. Candes et al., 2006; Ying et al., 2005)), which is linear. We can now rewrite the above equations as  0  300  400  500  -0.50 -0.50  = =  b2 b1  -0.25  ples  (3)  0  -0.25  k1  0.25  0.50  (b)  Figure 1: Curvelet examples. (a)-(b) spatial and frequency representation of four different curvelets in the spatial domain at three different scales and in the Fourier domain (b). The Fourier spectrum is overlaid with the dyadic partitioning of curvelets in the frequency. Each wedge corresponds to the compact frequency support of a curvelet. The curvelet transform itself consists for each wedge of a collection of shifted curvelets on a scale-dependent grid. Curvelets are not strictly localized in the space but they are of rapid decay. They oscillate in one direction and are smooth in the other(s).  Ax2 + n2 Ax1 + n − n2 .  (4)  Thus, our objective is to estimate the curvelet coefficients of the primaries x1 and multiples x2 , given the predictions b1 and b2 . These predictions could be given by SRME (see e.g. Berkhout and Verschuur, 1997) or by other means. From a Bayesian point of view, this means that our objective is to find x1 and x2 such that the conditional probability, P(x1 , x2 |b1 , b2 ), is maximized. In other words, we aim to maximize P(x1 , x2 |b1 , b2 )  = =  P(x1 ,x2 )P(b1 |x1 ,x2 )P(b2 |b1 ,x1 ,x2 ) P(b1 ,b2 ) P(x1 ,x2 )P(n)P(n2 ) . P(b1 ,b2 )  (5)  Since both b1 and b2 are known, we try to find the curvelet coefficients for the primaries and multiples that maximize the posterior probability in Eq. 5, i.e., arg max P(x1 , x2 |b1 , b2 ) = arg max P(x1 , x2 )P(n)P(n2 ). (6) x1 ,x2  x1 ,x2  If we assume independent identically distributed white Gaussian noise distributions for n and n2 , with possibly different  Curvelet-based primary multiple separation variances, and independent weighted Laplacian priors for x1 and x2 , then it can be shown that equation (6) reduces to arg maxx1 ,x2 P(x1 , x2 |b1 , b2 ) = arg minx1 ,x2 f (x1 , x2 ) with f (x1 , x2 )  =  x1 1,w1 + x2 1,w2 + Ax2 − b2 +η A(x1 + x2 ) − (b1 + b2 ) 22 ,  2 2  (7)  P where xi 1,wi = µ |wi,µ xi,µ |, µ ∈ M is a weighted 1 -norm. The weights are given by the predictions for the two signal components, i.e., w1 = λ1 |AT b2 | and w2 = λ2 |AT b1 | with AT the forward curvelet transform, and λ1 and λ2 parameters that control the sparsity of the coefficient vectors. The parameter η controls the tradeoff between fitting the total data and fitting the predicted multiples. Hence, the maximization of the posterior probability, P(x1 , x2 |b1 , b2 ), leads to a minimization problem involving the weighted 1 -norms of the coefficients and the following data misfits: the misfit between the predicted and estimated multiples and between the sum of the estimated primaries and multiples and the observed total data. As previously mentioned, the above formalism is similar to approaches in the field of signal processing where Laplacian priors have been used extensively (see e.g. Zibulevsky and Pearlmutter, 2001; Li et al., 2004). Long-tailed Laplacian or Cauchy-distributions also have extensively been used in (exploration) seismology (see e.g. Sacchi and Ulrych, 1996). What is new in our approach, is the use of the sparsity promoting curvelet transform and the use of weights to formulate the primary-multiple separation problem. Moreover, the above approach differs from earlier work (see e.g. Herrmann and Verschuur, 2005; Herrmann et al., 2007) since it adds additional constraints that guarantee the consistency of the estimated multiples with the predicted multiples which are the input to our algorithm. This additional constraint prevents all the signal’s energy from ending up in one of the two unknown coefficients vectors.  SEPARATION ALGORITHM The primary-multiple separation thus entails the minimization of the objective function f (x1 , x2 ) for appropriately chosen λ1 , λ2 , η. To minimize f (x1 , x2 ) in equation 7 above, we devised an iterative thresholding algorithm in the spirit of the one developed by Daubechies et al. (2003). The proposed algorithm provably converges to the minimizer of f (x1 , x2 ), provided that all the weights in the vectors wi , i ∈ {1, 2} are strictly positive. Starting from any initial estimates x01 and x02 of x1 and x2 , the iterations of the algorithm proceed as xn+1 1  =  xn+1 2  =  ` T ´ A b2 − AT Axn2 + AT b1 − AT Axn1 + xn1 h ´i η ` T S w2 A b1 − AT Axn1 AT b2 − AT Axn2 + xn2 + η+1 S w1 2η  (8)  2(1+η)  where Sα is the soft-thresholding operator acting elementwise as (9) Sαµ (vµ ) = sgn(vµ ) · max(0, |vµ | − |αµ |).  EXPERIMENT AND RESULTS To test the performance of our algorithm, it is applied on one shot of a three-dimensional data volume generated by a subsurface velocity model with two-dimensional inhomogeneities. This velocity model consists of a high-velocity layer, which represents salt, surrounded by sedimentary layers and a water bottom that is not completely flat. Applying an acoustic finite-difference modeling algorithm, 361 shots with 361 receivers are simulated on a fixed receiver spread with receivers located from 0 to 5400 m with steps of 15 m. The complete prestack dataset can be represented as a three-dimensional volume along the shot, receiver and time coordinates from which multiples are predicted and removed with the iterative matchedfilter method introduced by Berkhout and Verschuur (1997) (see figure 2). The results of the above method are compared to those obtained by using the method of Herrmann et al. (2007) as well as to those obtained by the proposed algorithm (applied here to one shot of the data). The results can be seen in figure 3. The comparison shows an improved removal of late arriving multiple energy by using the curvelet-based algorithms. However, the thresholded result of Herrmann et al. (2007) shows a large remaining multiple around 0.6 s and an offset of 2400 m which is no longer visible in the result obtained with our iterative method. Moreover, the artifacts that the thresholded method generates appear much less severe in the proposed method.  DISCUSSION AND CONCLUSION We have presented a primary-multiple separation algorithm which uses seismic data as well as predictions of the primaries and multiples as input and exploits the sparsity of the components in the curvelet domain to perform the separation. The algorithm, which minimizes an objective function involving weighted 1 -norms of the curvelet coefficients as well as misfits with the total data and the predicted multiples was derived from a Bayesian formulation of the separation problem. Moreover, the algorithm was applied to a synthetic data set and the results were compared to those of SRME and of the first iteration of curvelet-based method proposed by Herrmann et al. (2007). Our algorithm was shown to produce favorable results, by removing late arriving multiple energy when compared to SRME results. It was also shown to remove large multiples that remained upon applying the single-threshold algorithm and also to produce less artifacts. The improved performance relative to the SRME method can be partly attributed to the ability of the curvelet transform to represent seismic data, while the improved results of the proposed algorithm relative to the thresholded results can be attributed to the use of the new cost function f (x1 , x2 ) which mitigates the effect of artifacts and is simultaneously relatively robust to the choice of control parameters.  Curvelet-based primary multiple separation  (a)  (b)  Figure 2: Synthetic data and SRME-predicted multiples. (a) Total data. (b) SRME-predicted multiples .  (a)  (b)  (c)  Figure 3: Curvelet-based primary estimation. (a) SRME-predicted primaries (b) Estimated primaries by a single thresholding as in Herrmann et al. (2007). (c) Estimated primaries after 10 iterations of our separation algorithm (equation 8) with λ1 = 0.8, λ2 = 1.2 and η = 1.2. Comparing the SRME predicted primaries with the curvelet-based estimates show an improved removal of late arriving multiple energy. In addition, the thresholded result gives rise to a large remaining multiple around 0.6 s and an offset of 2400 m which is no longer visible in the result obtained with our iterative method. Also the artifacts in this region are less severe.  Curvelet-based primary multiple separation REFERENCES Berkhout, A. J. and D. J. Verschuur, 1997, Estimation of multiple scattering by iterative inversion, part I: theoretical considerations: Geophysics, 62, 1586–1595. Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899. Daubechies, I., 1992, Ten lectures on wavelets: SIAM. Daubechies, I., M. Defrise, and C. De Mol, 2003, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Hennenfent, G. and F. J. Herrmann, 2006, Seismic denoising with non-uniformly sampled curvelets: IEEE Comp. in Sci. and Eng., 8, 16–25. Herrmann, F. J., U. Boeniger, and D.-J. E. Verschuur, 2007, Nonlinear primary-multiple separation with directional curvelet frames: Geoph. J. Int. To appear. Herrmann, F. J. and D. J. Verschuur, 2005, Robust curvelet-domain primary-multiple separation with sparseness constraints: Presented at the EAGE 67th Conference & Exhibition Proceedings. Li, Y., A. Cichocki, and S. ichi Amari, 2004, Analysis of sparse representation and blind source separation: Neural Comput., 16, 1193–1234. Sacchi, M. and T. Ulrych, 1996, Estimation of the discrete fourier transform, a linear inversion approach: Geophysics, 61, 1128– 1136. Ying, L., L. Demanet, and E. J. Cand´es, 2005, 3D discrete curvelet transform: Wavelets XI, Expanded Abstracts, 591413, SPIE. Zibulevsky, M. and B. A. Pearlmutter, 2001, Blind source separation by sparse decomposition in a signal dictionary: Neural Computation, 13, 863–882.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items