UBC Faculty Research and Publications

Seismic data processing with curvelets: a multiscale and nonlinear approach. Herrmann, Felix J.; Wang, Deli; Hennenfent, Gilles; Moghaddam, Peyman P. 2007-03-10

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

Item Metadata


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

Full Text

Seismic data processing with curvelets: a multiscale and nonlinear approachFelix J. Herrmannasteriskmath, EOS-UBC, Deli Wang?, Jilin University and Gilles Hennenfentasteriskmath and Peyman MoghaddamasteriskmathSUMMARYIn this abstract, we present a nonlinear curvelet-based sparsity-promoting formulation of a seismic processing flow, consist-ing of the following steps: seismic data regularization andthe restoration of migration amplitudes. We show that thecurvelet?s wavefront detection capability and invariance underthe migration-demigration operator lead to a formulation thatis stable under noise and missing data.INTRODUCTIONIn this abstract, recent applications of the discrete curvelettransform (see e.g. Candes et al., 2006; Hennenfent and Herr-mann, 2006b) are presented that range from data recovery fromacquisitions with large percentages of traces missing to therestoration of migration amplitudes. Our approach derives fromtwo properties of curvelets, namely the detection of wave-fronts, without prior information on the positions and localdips (see e.g. Candes et al., 2006; Hennenfent and Herrmann,2006b) and the relative invariance of curvelets under wavepropagation (see e.g. Cand` and Demanet, 2005). Theseproperties render this transform suitable for a robust formula-tion of data regularization (see e.g. Hennenfent and Herrmann,2006a; Herrmann and Hennenfent, 2007); primary-multipleseparation (Herrmann et al., 2007); of migration-amplitude re-covery (Herrmann et al., 2006) and of wavefield extrapolation(Lin and Herrmann, 2007). All these methods derive fromsparsity in the curvelet domain that is a consequence of theabove properties. This sparsity corresponds to a rapid decayfor the magnitude-sorted curvelet coefficients and facilitates aseparation of (coherent) ?noise? and ?signal?. This separationunderlies the successful applications of this transform to explo-ration seismology (see e.g. Hennenfent and Herrmann, 2006b;Herrmann et al., 2007).CURVELETSCurvelets are localized ?little plane-waves? (see e.g. Hennen-fent and Herrmann, 2006b) that are oscillatory in one directionand smooth in the other direction(s). They are multiscale andmulti-directional. Curvelets have an anisotropic shape ? theyobey the so-called parabolic scaling relationship, yielding awidth proportional length2 for the support of curvelets. This anisotropicscaling is optimal for detecting wavefronts and explains theirhigh compression rates on seismic data and images (Candeset al., 2006; Hennenfent and Herrmann, 2006b; Herrmann et al.,2007).Curvelets represent a specific tiling of the 2-D/3-D frequencyplane into strictly localized multiscale and multi-angular wed-ges. Because the directional sampling increases every-otherscale doubling, curvelets become more anisotropic at finer scales.Curvelets compose an arbitrary column vector f, with the re-ordered samples, according to f =CTCf with C and CT , theforward/inverse discrete curvelet transform matrices (definedby the fast discrete curvelet transform, FDCT, with wrappingCandes et al., 2006; Ying et al., 2005). The symbol T repre-sents the transpose, which is equivalent to the pseudoinversefor our choise of discrete curvelet transform, which is a tightframe with a moderate redundancy (a factor of roughly 8 ford = 2 and 24 for d = 3 with d the number of dimensions).Tight frames (see e.g. Daubechies, 1992) are signal represen-tations that preserve energy. Consequently, CTC = I with Ithe identity matrix. Because of the redundancy, the converse isnot the identity, i.e., CCT negationslash=I .SPARSITY-PROMOTING INVERSIONTo exploit curvelets, (in)complete and noisy measurements arerelated to a sparse curvelet coefficient vector, x0, according toy =Ax0 +n (1)with y a vector with noisy and possibly incomplete measure-ments; A the synthesis matrix that includes the inverse curvelettransform (CT ); and n, a zero-centered white Gaussian noise.The matrix A is a wide rectangular matrix, so the vector x0can not readily be calculated from the measurements, becausethere exist infinitely many vectors that match y.Recent work in ?compressive sampling? (Cand` et al., 2006;Donoho, 2006) has shown that rectangular matrices can stablybe inverted by solving a nonlinear sparsity promoting program(Elad et al., 2005). These inversions require a fast decay forthe magnitude-sorted curvelet coefficients. Following these re-sults, the vector x0 can be recovered from noise-corrupted andincomplete data. Sparsity-promoting norm-one penalty func-tionals are not new to the geosciences (see for instance theseminal work of Claerbout and Muir (1973), followed by manyothers). New are (i) the curvelet transform that obtains nearoptimal theoretical and empirical (Candes et al., 2006; Hen-nenfent and Herrmann, 2006b) compression rates on seismicdata and images and (ii) the theoretical understanding of theconditions for a successful recovery. In this work, problemsin seismic imaging and processing are solved by the norm-onenonlinear program:Pepsilon :(e =argminx bardblxbardbl1 s.t. bardblAx-ybardbl2 <=<epsilonef =SSST ex (2)in which epsilon is a noise-dependent tolerance level. Eq. 2 is gen-eral and the (curvelet-based) synthesis matrix, A, and the in-verse sparsity transform, ST , are defined in accordance withthe application. The vectorerepresents the estimated solution(denoted by the symbol e). The above nonlinear program issolved with a threshold-based cooling method following ideasfrom Figueiredo and Nowak (2003) and Elad et al. (2005).SEISMIC DATA RECOVERYSeismic data acquisition is often based on equally-spaced sam-pling. Indeed, when Nyquist?s sampling theorem is met, equidis-tant sampling allows for a perfect reconstruction of bandwidth-limited signals. Unfortunately, sampling rates are often inad-equate, which leads to difficult to remove aliasing. For theCurvelet-based seismic processingsame average sampling interval, random subsampling leads toan easily denoiseable spectrum, as long as the to-be-recoveredsignal is compressible (Hennenfent and Herrmann, 2007). Ran-dom (compressive) sampling has been a topic of intense scien-tific debate (see e.g. Sun et al., 1997) and the response byVermeer (1998). What differs in our approach is the nonlinearrecovery from incomplete data. Recovery from compressivesampling depends on sparsity and mixing. Mixing turns harm-ful aliasing into relatively harmless noise and depends on therandomness of the acquisition and the incoherence (max corre-lation) between the measurement basis (Diracs) and curvelets.The better the mixing the better the recovery (Donoho et al.,2006; Herrmann and Hennenfent, 2007).Curvelet-based recovery: In our formulation, seismic dataregularization involves the solution of Pepsilon with A :=RCT , S :=C, given incomplete data, y = Rf, with f the fully sampleddata and R the picking matrix. In recent years, the authorsrepeatedly reported on succesful curvelet-based recovery ofseismic data (see e.g. Herrmann, 2005; Hennenfent and Herr-mann, 2006a, 2007). Compared to other methods, such assparse Fourier recovery (Sacchi and Ulrych, 1996; Zwartjesand Gisolf, 2006) and plane-wave destruction (Fomel and Gui-tton, 2006), curvelet-based methods work for data with con-flicting dips. Fig. 1, shows an unfavorable recovery example,where aliased groundroll is recovered by interpolation from a5m grid to a grid of 2.5m, where the groundroll is no longeraliased. Large amplitude contrast and regular subsamplingmake this example difficult.Focused recovery: Combining the non-adaptive curvelet trans-form with the data-adaptive focal transform (Berkhout and Ver-schuur, 2006), leads to a powerful formulation where data isfocused by inverting the primary operator (= a multidimen-sional ?convolution? with an estimate of the major primaries).During this curvelet-regularized inversion of the primary op-erator, deltaP, propagation paths that include the surface are re-moved, yielding a more focused wavefield and hence a morecompressed curvelet vector. This improved focussing is achievedby Pepsilon with the synthesis matrix A :=RdeltaPCT and inverse spar-sity transform ST := deltaPCT . The solution of Pepsilon now entailsthe inversion of deltaP, yielding the sparsest set of curvelet co-efficients that matches the incomplete data when ?convolved?with the primaries. By virtue of the improved compression,focusing improves the recovery (compare Fig 2(c) and 2(d)).MIGRATION-AMPLITUDE RECOVERYAfter processing seismic data can be modeled by the linearizedBorn scattering/demigration operator, d = Km +n, with theremaining nonlinear signal components and measurement er-rors modeled by Gaussian noise, n. A seismic image is cre-ated by applying the adjoint of the modeling operator. Thisimage, y = KT d or y = Psim+e, contains the imprint of clut-ter, e = KT n, and the Gramm (demigration-migration) opera-tor, Psi :=KT K, and serves as input to our amplitude recoveryscheme which estimates the reflectivity, m. The following ap-proximate identity is usedAAT r similarequalPsir (3)with r a reference vector. This expression follows from acurvelet decomposition of the Gramm operator, PsirsimilarequalCT DPsiCr.Following work by Guitton (2004) and more recently on op-timal scaling of reverse-time migration (W. W. Symes, per-sonal communication, 2007), this diagonal can be estimatedby solving a least-squares problem that uses the reference anddemigrated-migrated reference vector as input. The factoriza-tion in Eq. 3 is based on a synthesis matrix given by A :=CT Gwith G :=radicalDPsi. This factorization leads to y similarequalAx0 +e as anapproximate image representation that is amenable to Pepsilon . Af-ter solving for x0, the reflectivity is obtained by applying thesynthesis operator, ST := `AT ??, with ? the pseudoinverse, toe.The diagonal approximation serves two purposes. It approxi-mately corrects the amplitudes and it whitens the colored clut-ter. As the results in Fig. 3(b) indicate, a recovery with Pepsilon ,leads to a stable (under noise) recovery of the imaged reflec-tivity. These improvements (cf. Fig. 3(a) and 3(b)) are obtainedby a penalty functional that jointly promotes the curvelet spar-sity and the continuity along the reflectors. Results for theSEG AA? dataset (O?Brien and Gray, 1997; Aminzadeh et al.,1997) are summarized in Fig. 3. These results are obtained forlinearized Born data, modeled with the adjoint of reverse-timemigration with optimal checkpointing (W. W. Symes, personalcommunication, 2007). The estimated images show a nice am-plitude recovery and clutter removal for reflection data with asignal-to-noise ratio (SNR) of only 3dB.DISCUSSION AND CONCLUSIONSThe seismic data regularization derived its performance fromcurvelet compression, random subsampling and the inclusionof an additional focusing operator. The approximate invarianceunder demigration-migration allowed for a computationally-efficient migration amplitude recovery scheme. Since our schemeis based on a phase-space approximation, our scaling is lessrestrictive with respect to conflicting dips. The successful ap-plication of curvelets, juxtaposed by sparsity-promoting inver-sion, opens a range of new perspectives on seismic data pro-cessing, wavefield extrapolation and imaging. Because of theirsingular wavefront detection capability, curvelets represent inour vision an ideal domain for future developments in explo-ration seismology.ACKNOWLEDGMENTSThe authors would like to thank Eric Verschuur and Chris Stolkfor their input. We also would like to thank the authors ofCurveLab for making their codes available and William Symesfor his reverse-time migration code. The examples presentedwere prepared with Madagascar (rsf.sourceforge.net/),supplemented by SLIMPy operator overloading, developed bySean Ross Ross. ExxonMobil and Norsk Hydro are thankedfor making the field datasets available. Eric Verschuur andM. O?Brien, S. Gray and J. Dellinger are thanked for providingthe synthetic multiple dataset for the SEG AA? available. Thiswork was in part financially supported by the NSERC Dis-covery (22R81254) and CRD Grants DNOISE (334810-05) ofF.J.H. and was carried out as part of the SINBAD project withsupport, secured through ITF, from BG Group, BP, Chevron,ExxonMobil and Shell.Curvelet-based seismic processing(a) (b)Figure 1: Illustrations of curvelet-based seismic data recovery. (a) Original shot record with aliased groundroll sampled with aninterval of 5m. (b) Interpolated result, yielding a sample interval of 2.5m for which the groundroll is no longer aliased.(a) (b)(c) (d)Figure 2: Comparison between curvelet-based recovery by sparsity-promoting inversion with and without focusing. (a) Full realSAGA data volume. (b) Randomly subsampled data with 80% of the traces missing. (c) Curvelet-based recovery. (d) Curvelet-based recovery with focusing. Notice the significant improvement from the focusing with the primary operator.Curvelet-based seismic processing(a)(b)Figure 3: Image amplitude recovery for a migrated image calculated from noisy data (SNR 3dB). (a) Image with clutter. (b) Imageafter nonlinear recovery. The clearly visible non-stationary noise in (a) is mostly removed during the recovery while the amplitudesare also restored. Steeply dipping reflectors under the salt are also well recovered.Curvelet-based seismic processingREFERENCESAminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D Salt and Overthrust Model: Society of Exploration Geophysicists.Berkhout, A. J. and D. J. Verschuur, 2006, Focal transformation, an imaging concept for signal restoration and noise removal:Geophysics, 71.Cand` E., J. Romberg, and T. Tao, 2006, Stable signal recovery from incomplete and inaccurate measurements: Comm. PureAppl. Math., 59, 1207?1223.Cand` E. J. and L. Demanet, 2005, The curvelet representation of wave propagators is optimally sparse: Comm. Pure Appl. Math,58, 1472?1528.Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul.,5, 861?899.Claerbout, J. and F. Muir, 1973, Robust modeling with erratic data: Geophysics, 38, 826?844.Daubechies, I., 1992, Ten lectures on wavelets: SIAM.Donoho, D., Y. Tsaig, I. Drori, and J.-L. Starck, 2006, Sparse solution of underdetermined linear equations by stagewise orthonor-mal matching pursuit. Submitted for publication.Donoho, D. L., 2006, Compressed sensing: IEEE Trans. Inform. Theory, 52, 1289?1306.Elad, M., J. L. Starck, P. Querre, and D. L. Donoho, 2005, Simulataneous Cartoon and Texture Image Inpainting using Morpholog-ical Component Analysis (MCA): Appl. Comput. Harmon. Anal., 19, 340?358.Figueiredo, M. and R. Nowak, 2003, An EM algorithm for wavelet-based image restoration: IEEE Trans. Image Processing, 12,906?916.Fomel, S. and A. Guitton, 2006, Regularizing seismic inverse problems by model reparameterization using plane-wave construction:Geophysics, 71, A43?A47.Guitton, A., 2004, Amplitude and kinematic corrections of migrated images for nonunitary imaging operators: Geophysics, 69,1017?1024.Hennenfent, G. and F. Herrmann, 2006a, Application of stable signal recovery to seismic interpolation: Presented at the SEGInternational Exposition and 76th Annual Meeting.???, 2007, Irregular sampling: from aliasing to noise: Presented at the EAGE 69th Conference & Exhibition.Hennenfent, G. and F. J. Herrmann, 2006b, Seismic denoising with non-uniformly sampled curvelets: IEEE Comp. in Sci. andEng., 8, 16?25.Herrmann, F. J., 2005, Robust curvelet-domain data continuation with sparseness constraints: Presented at the EAGE 67th Confer-ence & Exhibition Proceedings.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 G. Hennenfent, 2007, Non-parametric seismic data recovery with curvelet frames. Submitted for publication.Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2006, Sparsety- and continuity-promoting seismic imaging with curvelet frames.In revision.Lin, T. and F. J. Herrmann, 2007, Compressed wavefield extrapolation. in revision.O?Brien, M. and S. Gray, 1997, Can we image beneath salt?: Leading Edge, 15, 17?22.Sacchi, M. and T. Ulrych, 1996, Estimation of the discrete fourier transform, a linear inversion approach: Geophysics, 61, 1128?1136.Sun, Y., G. T. Schuster, and K. Sikorski, 1997, A quasi-Monte Carlo approach to 3-D migration: Theory: Geophysics, 62, 918?928.Vermeer, G., 1998, On: ?A quasi-monte carlo approach to 3-D migration: Theory,? (Y. Sun, G. T. Schuster, and K. Sikorski,GEOPHYSICS, 62): Geophysics, 63, 1475.Ying, L., L. Demanet, and E. J. Cand? 2005, 3D discrete curvelet transform: Wavelets XI, Expanded Abstracts, 591413, SPIE.Zwartjes, P. and A. Gisolf, 2006, Fourier reconstruction of marine-streamer data in four spatial coordinates: Geophysics, 71,V171?V186.


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