UBC Faculty Research and Publications

Curvelet-domain multiple elimination with sparseness constraints. Herrmann, Felix J.; Verschuur, Eric 2004-02-21

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

Item Metadata


52383-SEGM2004.pdf [ 1.12MB ]
JSON: 52383-1.0107388.json
JSON-LD: 52383-1.0107388-ld.json
RDF/XML (Pretty): 52383-1.0107388-rdf.xml
RDF/JSON: 52383-1.0107388-rdf.json
Turtle: 52383-1.0107388-turtle.txt
N-Triples: 52383-1.0107388-rdf-ntriples.txt
Original Record: 52383-1.0107388-source.json
Full Text

Full Text

Curvelet-domain multiple elimination with sparseness constraintsFelix J. Herrmann1 and Eric Verschuur21 Department of Earth and Ocean Sciences, University of British Columbia, Canada2 Faculty of Applied Sciences, Delft University of Technology, The NetherlandsAbstractPredictive multiple suppression methods consist of two main steps: aprediction step, in which multiples are predicted from the seismic data,and a subtraction step, in which the predicted multiples are matchedwith the true multiples in the data. The last step appears crucial inpractice: an incorrect adaptive subtraction method will cause multiplesto be sub-optimally subtracted or primaries being distorted, or both.Therefore, we propose a new domain for separation of primaries andmultiples via the Curvelet transform. This transform maps the datainto almost orthogonal localized events with a directional and spatial-temporal component. The multiples are suppressed by thresholding theinput data at those Curvelet components where the predicted multipleshave large amplitudes. In this way the more traditional filtering ofpredicted multiples to fit the input data is avoided. An initial field dataexample shows a considerable improvement in multiple suppression.IntroductionIn complex areas move-out filtering multiple suppression techniquesmay fail because underlying assumptions are not met. Several attemptshave been made to address this problem by either extending move-outdiscrimination methods towards 3D complexities (e.g. by introducingapex-shifted hyperbolic transforms [9]) or by coming up with matchingtechniques in the wave-equation based predictive methods [see e.g.19, 1]. Least-squares matching the predicted multiples in time andspace overlapping windows, [18] provides a straightforward subtractionmethod, where the predicted multiples are matched to the true multiplesfor 2-D input data. Unfortunately, this matching procedure fails whenthe underlying 2D assumption are severely violated. There have beenseveral attempts to address this issue and the proposed solutions rangefrom including surrounding shot positions [13] to methods based onmodel- [16], data-driven [12] time delays and separation of predictedmultiples into (in)-coherent parts [14]. Even though these recentadvances in adaptive subtraction and other techniques have improvedthe attenuation of multiples, these methods continue to suffer from (i) arelative strong sensitivity to the accuracy of the predicted multiples; (ii)creation of spurious artifacts or worse (iii) a possible distortions of theprimary energy. For these situations, subtraction techniques based on adifferent concept are needed to complement the processor?s tool box.The method we are proposing here holds the middle between two com-plementary approaches common in multiple elimination: prediction incombination with subtraction and filtering [17, 9]. Whereas the firstapproach aims to predict the multiples and then subtract, the secondapproach tries to find a domain in which the primaries and multiplesseparate, followed by some filtering operation and reconstruction. Ourmethod is not distant from either since it uses the predicted multiples tonon-linearly filter data in a domain spanned by almost orthogonal andlocal basis functions. We use the recently developed Curvelet trans-form [see e.g. 3], that decomposes data into basis functions that not onlyobtain optimal sparseness on the coefficients and hence reduce the di-mensionality of the problem but which are also local in both locationand angle/dip, facilitating the definition of non-linear estimators basedon thresholding. Main assumption of this proposal is that multiples andprimaries have locally a different temporal, spatial and dip behavior, andtherefore map into different areas in the Curvelet domain. Multiples giverise to large Curvelet coefficients in the input and these coefficients canbe muted by our estimation procedure when the threshold is set accord-ing to the Curvelet transform of the predicted multiples. As such, oursuppression technique has at each location in the transformed domainone parameter, namely the threshold yielded by the predicted multiple,beyond which the input data is suppressed. In that sense, our proce-dure is similar to the ones proposed by [20] and [17], although the latteruse the non-localized FK/Radon domains for their separation while weuse localized basis functions and non-linear estimation by thresholding.Non-locality and non-optimality in their approximation renders the firstfiltering techniques less effective because primaries and multiples willstill have a considerable overlap. The Curvelet transform is able to makea local discrimination between interfering events with different temporaland spatial characteristics.The inverse problem underlying adaptive subtractionVirtually any linear problems in seismic processing and imaging can beseen as a special case of the following generic problem: how to obtainm form data d in the presence of (coherent) noise n:d = Km + n, (1)in which, for the adaptive subtraction problem, K and Kasteriskmath, represent therespective time convolution and correlation with an effective wavelet Phithat minimizes the following functional [see e.g. 8]minPhi = bardbld -Phi t mbardblp, (2)where d represents data with multiples; m the predicted multiples and? = minPhi bardbld -Phi asteriskmath mbardblp the primaries only data, obtained by mini-mizing the Lp-norm. For p = 2, Eq. 2 corresponds to a standard least-squares problem where the energy of the primaries is minimized.Seismic data and images are notoriously non-stationary and non-Gaussian, which may give rise to a loss of coherent events in the de-noised component when employing the L2-norm. The non-stationarityand non-Gaussianity have with some success been addressed by solvingEq. 2 within sliding windows and for p = 1 [8]. In this paper, we followa different strategy replacing the above variational problem by? : minm 12bardblC-1/2n (d -m) bardbl22 + ?J(m). (3)In this formulation, the necessity of estimating a wavelet, Phi has beenomitted. Data with multiples is again represented by d but now m isthe primaries only model on which an additional penalty function is im-posed, such as a L1-norm, i.e. J(m) = bardblmbardbl1. The noise term is nowgiven by the predicted multiples and this explains the emergence of thecovariance operator, whose kernel is given byCn = E{nnT }. (4)Question now is can we find a basis-function representation which is(i) sparse and local on the model (primaries only data), m, data, dand predicted multiples n and (ii) almost diagonalizes the covarianceof the model and the predicted multiples. Answer is affirmative, eventhough the condition of locality is non-trivial. For instance, decomposi-tions in terms of principle components/Karhuhnen-Lo? e-basis (KL) orindependent components may diagonalize the covariance operators butthese decomposition are generally not local and not necessary the samefor both multiples and primaries.By selecting a basis-function decomposition that is local, sparse and al-most diagonalizing (i.e. C?n approxequal diag{C?n} = diag{G2?n}), we find asolution for the above denoising problem that minimizes (mini) the max-imal (max) mean-square-error given the worst possible Bayesian priorfor the model [15, 7]. This minimax estimator minimizes?0 : minm 12bardblG-1?nparenleftBig?d - ?parenrightBigbardbl22 (5)by a simple diagonal non-linear thresholding operator? = B?GThetalambda parenleftbigG-1Bdparenrightbig = B?ThetalambdaGparenleftBig?dparenrightBig. (6)The threshold is set to lambdadiag{G} (we dropped the subscript ?) with lambda anadditional control parameter, which sets the confidence interval (e.g. theconfidence interval is 95 % for lambda = 3) (de)-emphasizing the threshold-ing. This thresholding operation removes the bulk of the multiples byshrinking the coefficients that are smaller then the noise to zero. Thesymbol ? indicates (pseudo)-inverse, allowing us to use Frames ratherthen orthonormal bases. Before applying this thresholding to actualdata, let us first focus on the selection of the appropriate basis-functiondecomposition that accomplishes the above task of replacing the adap-tive subtraction problem given in Eq. 2 by its diagonalized counterpart(cf. Eq. 6).The basis functionsCurvelets as proposed by [3], constitute a relatively new family ofnon-separable wavelet bases that are designed to effectively representseismic data with reflectors that generally tend to lie on piece-wisesmooth curves. This property makes Curvelets suitable to representevents in seismic whether these are located in shot records or timeslices. For these type of signals, Curvelets obtain nearly optimalsparseness, because of (i) the rapid decay for the reconstruction error asa function of the largest coefficients; (ii) the ability to concentrate thesignal?s energy in a limited number of coefficients; (iii) the ability tomap noise and signal to different areas in the Curvelet domain. So howdo Curvelets obtain such a high non-linear approximation rate? Withoutbeing all inclusive [see for details 4, 2, 3, 6], the answer to this questionlies in the fact that Curvelets are? multi-scale, i.e. they live in different dyadic corona (see for moredetail [3] or the other contributions of the first author to the pro-ceedings of this conference) in the FK-domain.? multi-directional, i.e. they live on wedges within these corona. SeeFig. 1.? anisotropic, i.e. they obey the following scaling law width proportionallength2.? directional selective with # orientations proportional 1scale .? local both in (x,t)) and KF.? almost orthogonal, they are tight frames with a moderate redun-dancy. Contourlets implement the pseudo-inverse in closed-formwhile Curvelets provide the transform and its adjoint, yielding apseudo-inverse computed by iterative Conjugate Gradients.Fig. 1: Top left: Curvelet partitioning of the frequency plane [modifiedfrom [5]]. Bottom right: Comparison of non-linear approximation ratesCurvelets and Wavelets [modified from [6]].Curvelets live in a wedges of the 2-D Fourier plane and become moredirectional selective and anisotropic for the higher frequencies. Theyare localized in both the space (or (x,t)) and spatial KF-domains andhave, as consequence of their partitioning, the tendency to align them-selves with curves/wavefronts. As such they are more flexible then a rep-resentation yielded by high-resolution Radon [as described by e.g. 17]because they are local and able to follow any piece-wise smooth curve.In Figure 1, a number of Curvelet properties are detailed [adapted from[5, 3, 6]].As the examples in the next section clearly demonstrate, the optimaldenoising capabilities for incoherent noise carry over to coherent noiseremoval provided we have reasonable accurate predictions for the noise,the multiples in this case. By choosing a threshold defined by the pre-dicted multiples, i.e.G = lambda|Bn|, (7)we are able to awe are able to adaptively decide whether a certain eventbelongs to primary or multiple energy. The n contains the predictedmultiples and lambda represents an additional control parameter which setsthe confidence interval (de)-emphasizing the thresholding.Results of SRME with thresholdingThe above methodology is tested on a field dataset from offshore Scot-land. This is a 2D line, which is known to suffer from strong surface-related multiples. Because the geology of the shallow sub-bottom layersis laterally complex, multiples have been observed to exhibit 3D char-acteristics. Applying surface-related multiple elimination to this datagives only an attenuation of surface multiples, but not a complete sup-pression. This lack of suppression can be observed for a shot recordshown in Fig. 2. When comparing the predicted multiples with the inputdata a good resemblance is observed in a global sense. However, whena matching filter is calculated, it appears that the predicted multiples donot coincide well enough with the true multiples, mainly due to 3D ef-fects. Locally the temporal and spatial shape of the predicted multiples- especially for the higher order multiples - differs from the true events.Estimating shaping filters in spatially and temporally varying, overlap-ping windows, as described by [18] can not resolve this mismatch com-pletely. The adaptive subtraction result is displayed in Fig. 2 in the thirdpanel. To improve this result, more freedom can be incorporated in thesubtraction process. However, note that the more the predicted multiplescan adapt to the input data, the higher the chance that primary energy willbe distorted as well, as the primary and multiple energy are not orthog-onal to each other in the least-squares sense. Multi-gather subtractionusing 9 surrounding multiple panels gives better suppression of the mul-tiples but leaves a lot of multiple remnants to be observed, especially forthe large offset at large travel times. The output gather of the Curveletmethod (see fourth panel in Fig. 2) looks much cleaner. Furthermore,clear events have been restored from interference with the multiples inthe lower left area. Also note the good preservation of a primary eventaround -1000 meter offset and 1.1 seconds. The amplitudes seem to beidentical to the ones in the input data, whereas the adaptive subtractionhas distorted these amplitudes to some extend. For the Curvelet-domainprocedure, the threshold value can be set by the user, thus creating moreor less suppression. Too weak thresholding leaves too much multiples inthe data and that a too strong thresholding procedure seems to removetoo much energy.The filtering via the Curvelet domain can also be applied within each 2Dcross-section through the seismic data volume. Time slices at t = 1.48through the prestack shot-offset data and predicted-multiple volumes aredisplayed in Fig. 3 (top row). Results after application of the adaptivesubtraction per shot record and the time-slice only Curvelet subtractionare included in Fig. 3 (bottom row). Again, the Curvelet results lookmuch cleaner than the least-squares subtraction result. Much of the noisyremnants at small offsets in the least-squares subtraction result have beenremoved. Notice also the improved suppression of the higher order peg-legs in the area around shot 900 and offset 1500 meter.ConclusionsIn this paper a new concept related to multiple subtraction has been de-scribed, based on the Curvelet transform. The Curvelet transform is analmost orthogonal transformation into local basis functions parameter-ized by their relate temporal- and spatial-frequency content. Because oftheir anisotropic shape, Curvelets are directional selective, i.e. they havelocal angle-discrimination capabilities. Our method uses the predictedmultiples from the surface-related multiple prediction method, as a guideto suppress the Curvelet coefficients related to the multiple events di-rectly in the original data. Therefore, the multiples are not actually sub-tracted, but areas in the Curvelet domain related to multiple energy aremuted. The success of this method depends on the assumption that pri-maries and multiples map into different areas in the Curvelet domain.Based on a field data example, we can conclude that the Curvelet-basedmultiple filtering is effective and is able to suppress multiples, while pre-serving primary energy. Especially in situations with clear 3D effects inthe 2D seismic data, it appears to perform better than the more tradi-tional least-squares adaptive subtraction methods. This success does notreally come as a surprise given the successful application of these tech-niques to the removal of noise colored by migration [10, 11] and to thecomputation of 4D difference cubes (see for details in both other contri-butions of the first to proceedings of this conference). Leaves us to hopefor future higher dimensional implementation of the Curvelet transformpossibly supplemented by constraint optimization imposing sparsenessconstraints (see the migration paper in these proceedings).AcknowledgmentsThe authors thank Elf Caledonia Ltd (now part of Total) for providingthe field data and Emmanuel Cand? and David Donoho for making anearly version of their Curvelet code (Digital Curvelet Transforms viaUnequispaced Fourier Transforms, presented at the ONR Meeting, Uni-versity of Minnesota, May, 2003) available for evaluation. This workwas in part financially supported by a NSERC Discovery Grant.References[1] A. J. Berkhout and D. J. Verschuur. Estimation of multiple scatter-ing by iterative inversion, part I: theoretical considerations. Geo-physics, 62(5):1586?1595, 1997.[2] E. J. Cand` and David Donoho. New tight frames of curveletsand optimal representations of objects with smooth singularities.Technical report, Stanford, 2002. URL http://www.acm.caltech.edu/?emmanuel/publications.html. sub-mitted.[3] E. J. Cand` and F. Guo. New multiscale transforms,minimum total variation synthesis: Applications to edge-preserving image reconstruction. Signal Processing, pages1519?1543, 2002. URL http://www.acm.caltech.edu/?emmanuel/publications.html.[4] Emmanual J. Cand` and David L. Donoho. Curvelets ? a surpris-ingly effective nonadaptive representation for objects with edges.Curves and Surfaces. Vanderbilt University Press, 2000.[5] Emmanual J. Cand` and David L. Donoho. Recovering Edges inIll-posed Problems: Optimality of Curvelet Frames. Ann. Statist.,30:784?842, 2000.[6] M. Do and M. Vetterli. Beyond wavelets, chapter Contourlets. Aca-demic Press, 2002.[7] David L. Donoho and Iain M. Johnstone. Minimax estimationvia wavelet shrinkage. Annals of Statistics, 26(3):879?921, 1998.URL citeseer.nj.nec.com/donoho92minimax.html.[8] A. Guitton. Adaptive subtraction of multiples using thel1-norm. Geophys Prospect, 52(1):27?27, 2004. URLhttp://www.blackwell-synergy.com/links/doi/10.1046/j.1365-2478.2004.00401.x/abs.[9] N. Hargreaves, B. verWest, R. Wombell, and D. Trad. Multipleattenuation using apex-shifted radon transform. pages 1929?1932,Dallas, 2003. SEG, Soc. Expl. Geophys., Expanded abstracts.[10] Felix J. Herrmann. Multifractional splines: application to seismicimaging. In Andrew F. Laine; Eds. Michael A. Unser, Akram Al-droubi, editor, Proceedings of SPIE Technical Conference onWavelets: Applications in Signal and Image Processing X, volume5207, pages 240?258. SPIE, 2003. URL http://www.eos.ubc.ca/?felix/Preprint/SPIE03DEF.pdf.[11] Felix J. Herrmann. Optimal seismic imaging with curvelets. In Ex-panded Abstracts, Tulsa, 2003. Soc. Expl. Geophys. URL http://www.eos.ubc.ca/?felix/Preprint/SEG03.pdf.[12] L. T. Ikelle and S. Yoo. An analysis of 2D and 3D inverse scatteringmultiple attenuation. pages 1973?1976, Calgary, 2000. SEG, Soc.Expl. Geophys., Expanded abstracts.[13] H. Jakubowicz. Extended subtraction of multiples. Private com-munication, Delft, 1999.[14] M. M. N. Kabir. Weighted subtraction for diffracted multiple at-tenuation. pages 1941?1944, Dallas, 2003. SEG, Soc. Expl. Geo-phys., Expanded abstracts.[15] S. G. Mallat. A wavelet tour of signal processing. Academic Press,1997.[16] W. S. Ross. Multiple suppression: beyond 2-D. part I: theory. DelftUniv. Tech., pages 1387?1390, Dallas, 1997. Soc. Expl. Geophys.,Expanded abstracts.[17] D. O. Trad. Interpolation and multiple attenuation with migrationoperators. Geophysics, 68(6):2043?2054, 2003.[18] D. J. Verschuur and A. J. Berkhout. Estimation of multiple scatter-ing by iterative inversion, part II: practical aspects and examples.Geophysics, 62(5):1596?1611, 1997.[19] D. J. Verschuur, A. J. Berkhout, and C. P. A. Wapenaar. Adap-tive surface-related multiple elimination. Geophysics, 57(9):1166?1177, 1992.[20] B. Zhou and S. A. Greenhalgh. Wave-equation extrapolation-basedmultiple attenuation: 2-d filtering in the f-k domain. Geophysics,59(10):1377?1391, 1994.Fig. 2: Comparison between adaptive least-squares (using time andspace varying least-squares filters)and Curvelet subtraction of predictedmultiples for one shot record. Input data (left), predicted multiples (sec-ond), least-squares subtraction (third) and the results by thresholdingCurvelet (right). Fig. 3: Time slice through the prestack data volume of all shot records.Top: Time slice through the input data with all multiples. Second: Slicethrough the 2D predicted surface-related multiples. Third: Time slicethrough the shot records after adaptive subtraction per shot record of thepredicted multiples. Bottom: The Curvelet equivalent of the adaptivesubtraction computed for a single time horizon.


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