UBC Faculty Research and Publications

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

You don't seem to have a PDF reader installed, try download the pdf

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 perspectiveRayan Saab, ECE-UBC, Deli Wang, Jilin University. ? ? Y?lmaz, Math-UBC and Felix J. Herrmann, EOS-UBCSUMMARYIn this abstract, we present a novel primary-multiple separationscheme which makes use of the sparsity of both primaries andmultiples in a transform domain, such as the curvelet trans-form, to provide estimates of each. The proposed algorithmutilizes seismic data as well as the output of a preliminary stepthat provides (possibly) erroneous predictions of the multiples.The algorithm separates the signal components, i.e., the pri-maries and multiples, by solving an optimization problem thatassumes noisy input data and can be derived from a Bayesianperspective. More precisely, the optimization problem can bearrived at via an assumption of a weighted Laplacian distribu-tion for the primary and multiple coefficients in the transformdomain and of white Gaussian noise contaminating both theseismic data and the preliminary prediction of the multiples,which both serve as input to the algorithm.INTRODUCTIONPredictive multiple suppression methods consist of two mainsteps: a prediction step, during which multiples are predictedfrom seismic data, and a primary-multiple separation step, dur-ing which the predicted multiples are ?matched? with the truemultiples in the data and subsequently removed. This sec-ond separation step, which we will call the estimation step,is crucial in practice. An incorrect separation will cause resid-ual multiple energy in the result or may lead to a distortionof the primaries, or both. To reduce these adverse effects,a new transform-domain method was proposed by Herrmannet al. (2007), where primaries and multiples are separated inthe curvelet domain rather than matched.This separation method is carried out on the basis of differ-ences in the multiscale and multidirectional characteristics ofthe primaries and multiples. The sparsity in the curvelet do-main of these two signal components is exploited by solvinga weighted lscript1-norm minimization problem that is designedto drive the two signal components apart. This separation isaccomplished by using the (SRME) predicted primaries andmultiples as weights. During the optimization, the two signalcomponents are separated by enhancing sparseness (throughweighted lscript1-norms) in the transform domain subject to fittingthe observed data as the sum of the separated components towithin a user-defined tolerance level.Even though good results were obtained with this method, prob-lems with the appropriate balancing of the weights were expe-rienced. Lack of a good balance yielded a ?separation? whereall the signal energy ends up in one of the two signal compo-nents. To tackle this issue, we propose an alternative formu-lation for the primary-multiple separation where an additionalconstraint is imposed on the coefficients for the multiples. Thisformulation follows from a Bayesian approach where sparsityfor the two signal components as well as data misfits are takeninto account.Our approach derives from techniques employed in the field ofsignal processing and geophysics where Bayesian approachessupplemented with sparsity-promoting Laplace/Cauchy prob-ability density functions are frequently used. Our method dif-fers from the earlier results in the literature by (i) using sparsitypromoting curvelets (see, e.g., Candes et al., 2006); (ii) us-ing inaccurate predictions for the to-be-separated signal com-ponents and (iii) providing a new iterative solution techniquebased on descent projections and soft thresholding.The outline of our abstract is as follows. First, we brieflydiscuss curvelets, our transform domain of choice, which hasbeen shown to admit sparse coefficients for seismic data (Can-des et al. (2006)). Next, Bayesian formulation of the primary-multiple separation problem is introduced, leading to the ob-jective function to be minimized. Finally, a separation algo-rithm that provably converges to a minimizer of the objectivefunction is given. Results of our algorithm applied to a syn-thetic shot with 3-D complexity are presented as well.CURVELETSCurvelets are localized ?little plane-waves? (see e.g. Hennen-fent and Herrmann, 2006, and ancillary electronic material)that are oscillatory in one direction and smooth in the other di-rection(s). They are multiscale and multi-directional. Curveletshave an anisotropic shape ? they obey the so-called parabolicscaling relationship, yielding a width proportionallength2 for the supportof curvelets. This anisotropic scaling is optimal for detect-ing wavefronts and explains their high compression rates onseismic data and images (Candes et al., 2006; Hennenfent andHerrmann, 2006; Herrmann et al., 2007).Curvelets represent a specific tiling of the 2-D/3-D frequencyplane into strictly localized multiscale and multi-angular wedges.Fig. 1(b)). Because the directional sampling increases every-other scale-doubling, curvelets become more anisotropic at finerscales. They become ?needle-like? as illustrated in Fig. 1(a).The smoothness in the Fourier domain leads to a rapid decayin the physical domain. Their effective support is given byellipsoids with a width proportional 2j/2 and length proportional 2j and an anglethetajl = 2pil2floorleft j/2floorright with j the scale, l the angular index with floorleftfloorrighta rounding off towards the nearest smaller integer. Curveletsare indexed by the multi-index ? :=( j, l, k) element M with M themulti-index set that runs over the scales, j, angles, l, and posi-tions k (see for details Candes et al., 2006; Ying et al., 2005).The number of wedges doubles every other scale-doubling andsuch that directional selectivity increases for finer scales (seeFig. 1(b)). Curvelets compose an arbitrary column vector f,with the reordered samples, according to f =CT Cf with C andCT , the forward/inverse discrete curvelet transform matrices(defined by the fast discrete curvelet transform, FDCT, withCurvelet-based primary multiple separationwrapping Candes et al., 2006; Ying et al., 2005). CT denotesthe transpose of C, which is equal to the pseudoinverse of Cfor our choise of discrete curvelet transform: the rows of Cconstitute a tight frame with frame bound 1 and with a mod-erate redundancy (a factor of roughly 8 for d = 2 and 24 ford = 3 with d the number of dimensions). Such tight frames(see e.g. Daubechies, 1992) are signal representations thatpreserve energy. Consequently, CT C = I with I the identitymatrix. However, CCT negationslash=I.(a)(b)Figure 1: Curvelet examples. (a)-(b) spatial and frequencyrepresentation of four different curvelets in the spatial domainat three different scales and in the Fourier domain (b). TheFourier spectrum is overlaid with the dyadic partitioning ofcurvelets in the frequency. Each wedge corresponds to thecompact frequency support of a curvelet. The curvelet trans-form itself consists for each wedge of a collection of shiftedcurvelets on a scale-dependent grid. Curvelets are not strictlylocalized in the space but they are of rapid decay. They oscil-late in one direction and are smooth in the other(s).BAYESIAN FORMULATIONHere, we approach the primary-multiple separation problemfrom a probabilistic perspective, where we impose an a pri-ori distribution for the unknown coefficients of the primariesand multiples. Estimates for the curvelet coefficients of theprimaries and multiples are obtained by using a Laplacian-type probability density function for the curvelet coefficientsand weights defined in terms of (SRME) predictions for thetwo signal components. The estimates are calculated withina Bayesian framework which leads to the minimization of aweighted sparsity-promoting functional. The choice of a Lap-lacian-type density function for the curvelet coefficients alsoserves as a sparsity promoting prior (see e.g. Zibulevsky andPearlmutter, 2001; Li et al., 2004) and is consistent with thehigh compression rates that curvelets display for seismic data(Candes et al., 2006; Hennenfent and Herrmann, 2006; Her-rmann et al., 2007).More precisely, letb =s1 +s2 +n (1)be the recorded total data set comprised of a noisy sum of pri-maries s1 and multiples s2 and letb2 =s2 +n2 (2)be the predicted multiples, assumed to be erroneous versions ofthe true underlying multiples. Moreover, let the multiples andprimaries admit (curvelet) transform coefficients satisfyingsi =Axi, {i =1, 2}, (3)where A is the inverse curvelet transform operator (defined bythe fast discrete curvelet transform (FDCT by wrapping, seee.g. Candes et al., 2006; Ying et al., 2005)), which is linear.We can now rewrite the above equations asb2 = Ax2 +n2b1 = Ax1 +n-n2. (4)Thus, our objective is to estimate the curvelet coefficients ofthe primaries x1 and multiples x2, given the predictions b1and 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 objec-tive is to find x1 and x2 such that the conditional probabil-ity, P(x1, x2|b1, b2), is maximized. In other words, we aim tomaximizeP(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 curveletcoefficients for the primaries and multiples that maximize theposterior probability in Eq. 5, i.e.,argmaxx1,x2P(x1, x2|b1, b2) =argmaxx1,x2P(x1, x2)P(n)P(n2). (6)If we assume independent identically distributed white Gaus-sian noise distributions for n and n2, with possibly differentCurvelet-based primary multiple separationvariances, and independent weighted Laplacian priors for x1and x2, then it can be shown that equation (6) reduces toargmaxx1,x2 P(x1, x2|b1, b2) =argminx1,x2 f(x1, x2)withf(x1, x2) = bardblx1bardbl1,w1 +bardblx2bardbl1,w2 +bardblAx2 -b2bardbl22+etabardblA(x1 +x2) -(b1 +b2)bardbl22, (7)where bardblxibardbl1,wi =P? |wi,?xi,?|, ? elementM is a weighted lscript1-norm.The weights are given by the predictions for the two signalcomponents, i.e., w1 =lambda1|AT b2| and w2 =lambda2|AT b1| with ATthe forward curvelet transform, and lambda1 and lambda2 parameters thatcontrol the sparsity of the coefficient vectors. The parametereta controls the tradeoff between fitting the total data and fittingthe predicted multiples.Hence, the maximization of the posterior probability,P(x1, x2|b1, b2), leads to a minimization problem involvingthe weighted lscript1-norms of the coefficients and the followingdata misfits: the misfit between the predicted and estimatedmultiples and between the sum of the estimated primaries andmultiples and the observed total data. As previously men-tioned, the above formalism is similar to approaches in thefield of signal processing where Laplacian priors have beenused extensively (see e.g. Zibulevsky and Pearlmutter, 2001;Li et al., 2004). Long-tailed Laplacian or Cauchy-distributionsalso have extensively been used in (exploration) seismology(see e.g. Sacchi and Ulrych, 1996). What is new in our ap-proach, is the use of the sparsity promoting curvelet transformand the use of weights to formulate the primary-multiple sep-aration problem. Moreover, the above approach differs fromearlier work (see e.g. Herrmann and Verschuur, 2005; Her-rmann et al., 2007) since it adds additional constraints thatguarantee the consistency of the estimated multiples with thepredicted multiples which are the input to our algorithm. Thisadditional constraint prevents all the signal?s energy from end-ing up in one of the two unknown coefficients vectors.SEPARATION ALGORITHMThe primary-multiple separation thus entails the minimizationof the objective function f(x1, x2) for appropriately chosenlambda1, lambda2, eta. To minimize f(x1, x2) in equation 7 above, we de-vised an iterative thresholding algorithm in the spirit of theone developed by Daubechies et al. (2003). The proposedalgorithm provably converges to the minimizer of f(x1, x2),provided that all the weights in the vectors wi, i element {1, 2} arestrictly positive. Starting from any initial estimates x01 and x02of x1 and x2, the iterations of the algorithm proceed asxn+11 = S w12eta`AT b2 -AT Axn2 +AT b1 -AT Axn1 +xn1?xn+12 = S w22(1+eta)hAT b2 -AT Axn2 +xn2 + etaeta+1 `AT b1 -AT Axn1?i (8)where Salpha is the soft-thresholding operator acting elementwiseasSalpha? (v?) =sgn(v?) ? max(0, |v?| -|alpha?|). (9)EXPERIMENT AND RESULTSTo test the performance of our algorithm, it is applied on oneshot of a three-dimensional data volume generated by a subsur-face velocity model with two-dimensional inhomogeneities.This velocity model consists of a high-velocity layer, whichrepresents salt, surrounded by sedimentary layers and a wa-ter bottom that is not completely flat. Applying an acousticfinite-difference modeling algorithm, 361 shots with 361 re-ceivers are simulated on a fixed receiver spread with receiverslocated from 0 to 5400 m with steps of 15 m. The completeprestack dataset can be represented as a three-dimensional vol-ume along the shot, receiver and time coordinates from whichmultiples are predicted and removed with the iterative matched-filter method introduced by Berkhout and Verschuur (1997)(see figure 2). The results of the above method are compared tothose obtained by using the method of Herrmann et al. (2007)as well as to those obtained by the proposed algorithm (ap-plied here to one shot of the data). The results can be seenin figure 3. The comparison shows an improved removal oflate arriving multiple energy by using the curvelet-based al-gorithms. However, the thresholded result of Herrmann et al.(2007) shows a large remaining multiple around 0.6s and anoffset of 2400m which is no longer visible in the result ob-tained with our iterative method. Moreover, the artifacts thatthe thresholded method generates appear much less severe inthe proposed method.DISCUSSION AND CONCLUSIONWe have presented a primary-multiple separation algorithmwhich uses seismic data as well as predictions of the primariesand multiples as input and exploits the sparsity of the compo-nents in the curvelet domain to perform the separation. Thealgorithm, which minimizes an objective function involvingweighted lscript1-norms of the curvelet coefficients as well as mis-fits with the total data and the predicted multiples was derivedfrom a Bayesian formulation of the separation problem. More-over, the algorithm was applied to a synthetic data set and theresults were compared to those of SRME and of the first iter-ation of curvelet-based method proposed by Herrmann et al.(2007). Our algorithm was shown to produce favorable re-sults, by removing late arriving multiple energy when com-pared to SRME results. It was also shown to remove large mul-tiples that remained upon applying the single-threshold algo-rithm and also to produce less artifacts. The improved perfor-mance relative to the SRME method can be partly attributed tothe ability of the curvelet transform to represent seismic data,while the improved results of the proposed algorithm relativeto the thresholded results can be attributed to the use of thenew cost function f(x1, x2) which mitigates the effect of ar-tifacts and is simultaneously relatively robust to the choice ofcontrol 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 asin Herrmann et al. (2007). (c) Estimated primaries after 10 iterations of our separation algorithm (equation 8) with lambda1 = 0.8,lambda2 =1.2 and eta =1.2. Comparing the SRME predicted primaries with the curvelet-based estimates show an improved removal oflate arriving multiple energy. In addition, the thresholded result gives rise to a large remaining multiple around 0.6s and an offsetof 2400m 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 separationREFERENCESBerkhout, 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 sparsityconstraint.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: Pre-sented 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? 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: NeuralComputation, 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