UBC Faculty Research and Publications

Curvelet-based primary-multiple separation from a Bayesian perspective Saab, Rayan 2008

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

Item Metadata


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

Full Text

Curvelet-based primary-multiple separation from a Bayesian perspective Rayan Saab, ECE-UBC, Deli Wang, Jilin University. Özgür Yılmaz, Math-UBC and Felix J. Herrmann, EOS-UBC 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 trans- form, 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 pri- maries 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 distribu- tion 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, dur- ing which the predicted multiples are ’matched’ with the true multiples 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 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 differ- ences in the multiscale and multidirectional characteristics of the primaries and multiples. The sparsity in the curvelet do- main 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, prob- lems with the appropriate balancing of the weights were expe- rienced. Lack of a good balance yielded a ’separation’ where all 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 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 prob- ability density functions are frequently used. Our method dif- fers from the earlier results in the literature by (i) using sparsity promoting 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 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 (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 objective function is given. Results of our algorithm applied to a syn- thetic shot with 3-D complexity are presented as well. CURVELETS Curvelets 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. 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 detect- ing 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 every- other scale-doubling, curvelets becomemore 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 = 2pil2b j/2c with j the scale, l the angular index with bc a rounding off towards the nearest smaller integer. Curvelets are indexed by the multi-index µ := ( j, l, k) ∈M withM the multi-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 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=CTCf 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 mod- erate 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, CTC = I with I the identity matrix. However, CCT 6= I . k1 0-0.25-0.50 0.25 0.50 -0.50 k 2 -0.25 0 0.25 0.500 100 200 300 400 500 0 100 200 300 400 500 Samples S a m p le s (a) k1 0-0.25-0.50 0.25 0.50 -0.50 k 2 -0.25 0 0.25 0.500 100 200 300 400 500 0 100 200 300 400 500 Samples S a m p le s (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 trans- form 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 oscil- late in one direction and are smooth in the other(s). BAYESIAN FORMULATION Here, we approach the primary-multiple separation problem from a probabilistic perspective, where we impose an a pri- ori 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 Laplacian- type 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 Lap- lacian-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; Her- rmann et al., 2007). More precisely, let b= s1+ s2+n (1) be the recorded total data set comprised of a noisy sum of pri- maries s1 and multiples s2 and let b2 = s2+n2 (2) 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 si = Axi, {i= 1,2}, (3) 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 b2 = Ax2+n2 b1 = 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 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 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., argmax x1,x2 P(x1,x2|b1,b2) = argmaxx1,x2 P(x1,x2)P(n)P(n2). (6) If we assume independent identically distributed white Gaus- sian 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 argmaxx1,x2 P(x1,x2|b1,b2) = argminx1,x2 f (x1,x2) with f (x1,x2) = ‖x1‖1,w1+‖x2‖1,w2 +‖Ax2−b2‖22 +η‖A(x1+x2)− (b1+b2)‖22, (7) where ‖xi‖1,wi = P µ |wi,µxi,µ |, µ ∈M is a weighted `1-norm. The weights are given by the predictions for the two signal components, i.e., w1 = λ1|ATb2| and w2 = λ2|ATb1| 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 men- tioned, 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 ap- proach, is the use of the sparsity promoting curvelet transform and the use of weights to formulate the primary-multiple sep- aration problem. Moreover, the above approach differs from earlier work (see e.g. Herrmann and Verschuur, 2005; Her- rmann 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 end- ing 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 de- vised 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 x 0 2 of x1 and x2, the iterations of the algorithm proceed as xn+11 = S w1 2η ` AT b2−ATAxn2 +AT b1−ATAxn1 +xn1 ´ xn+12 = S w2 2(1+η) h AT b2−ATAxn2 +xn2 + η η+1 ` AT b1−ATAxn1 ´i (8) where Sα is the soft-thresholding operator acting elementwise as Sαµ (vµ ) = sgn(vµ ) ·max(0, |vµ |− |αµ |). (9) 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 subsur- face velocity model with two-dimensional inhomogeneities. This velocity model consists of a high-velocity layer, which represents salt, surrounded by sedimentary layers and a wa- ter bottom that is not completely flat. Applying an acoustic finite-difference modeling algorithm, 361 shots with 361 re- ceivers 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 vol- ume along the shot, receiver and time coordinates from which multiples 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 to those 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 seen in figure 3. The comparison shows an improved removal of late 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 an offset of 2400m which is no longer visible in the result ob- tained 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 compo- nents 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 mis- fits with the total data and the predicted multiples was derived from a Bayesian formulation of the separation problem. More- over, the algorithm was applied to a synthetic data set and the results 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 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 ar- tifacts 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.6s and an offset of 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 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: 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és, 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:


Usage Statistics

Country Views Downloads
United States 1 2
Japan 1 0
City Views Downloads
Redmond 1 2
Tokyo 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}


Share to:


Related Items