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

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

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 approach Felix J. Herrmann∗, EOS-UBC, Deli Wang†, Jilin University and Gilles Hennenfent∗ and Peyman Moghaddam∗ SUMMARY In 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 and the restoration of migration amplitudes. We show that the curvelet’s wavefront detection capability and invariance under the migration-demigration operator lead to a formulation that is stable under noise and missing data. INTRODUCTION In this abstract, recent applications of the discrete curvelet transform (see e.g. Candes et al., 2006; Hennenfent and Herr- mann, 2006b) are presented that range from data recovery from acquisitions with large percentages of traces missing to the restoration of migration amplitudes. Our approach derives from two properties of curvelets, namely the detection of wave- fronts, without prior information on the positions and local dips (see e.g. Candes et al., 2006; Hennenfent and Herrmann, 2006b) and the relative invariance of curvelets under wave propagation (see e.g. Candès and Demanet, 2005). These properties render this transform suitable for a robust formula- tion of data regularization (see e.g. Hennenfent and Herrmann, 2006a; Herrmann and Hennenfent, 2007); primary-multiple separation (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 from sparsity in the curvelet domain that is a consequence of the above properties. This sparsity corresponds to a rapid decay for the magnitude-sorted curvelet coefficients and facilitates a separation of (coherent) ’noise’ and ’signal’. This separation underlies the successful applications of this transform to explo- ration seismology (see e.g. Hennenfent and Herrmann, 2006b; Herrmann et al., 2007). CURVELETS Curvelets are localized ’little plane-waves’ (see e.g. Hennen- fent and Herrmann, 2006b) 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, 2006b; Herrmann et al., 2007). Curvelets represent a specific tiling of the 2-D/3-D frequency plane into strictly localized multiscale and multi-angular wed- ges. Because the directional sampling increases every-other scale doubling, curvelets becomemore anisotropic at finer scales. Curvelets compose an arbitrary column vector f, with the re- ordered 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 wrapping Candes et al., 2006; Ying et al., 2005). The symbol T repre- sents the transpose, which is equivalent to the pseudoinverse for our choise of discrete curvelet transform, which is a tight frame with a moderate redundancy (a factor of roughly 8 for d = 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 I the identity matrix. Because of the redundancy, the converse is not the identity, i.e.,CCT 6= I . SPARSITY-PROMOTING INVERSION To exploit curvelets, (in)complete and noisy measurements are related to a sparse curvelet coefficient vector, x0, according to y= Ax0+n (1) with y a vector with noisy and possibly incomplete measure- ments; A the synthesis matrix that includes the inverse curvelet transform (CT ); and n, a zero-centered white Gaussian noise. The matrix A is a wide rectangular matrix, so the vector x0 can not readily be calculated from the measurements, because there exist infinitely many vectors that match y. Recent work in ’compressive sampling’ (Candès et al., 2006; Donoho, 2006) has shown that rectangular matrices can stably be inverted by solving a nonlinear sparsity promoting program (Elad et al., 2005). These inversions require a fast decay for the magnitude-sorted curvelet coefficients. Following these re- sults, the vector x0 can be recovered from noise-corrupted and incomplete data. Sparsity-promoting norm-one penalty func- tionals are not new to the geosciences (see for instance the seminal work of Claerbout andMuir (1973), followed by many others). New are (i) the curvelet transform that obtains near optimal theoretical and empirical (Candes et al., 2006; Hen- nenfent and Herrmann, 2006b) compression rates on seismic data and images and (ii) the theoretical understanding of the conditions for a successful recovery. In this work, problems in seismic imaging and processing are solved by the norm-one nonlinear program: Pε : (ex= argminx ‖x‖1 s.t. ‖Ax−y‖2 ≤ εef= STex (2) in which ε 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 with the application. The vectoref represents the estimated solution (denoted by the symbol e). The above nonlinear program is solved with a threshold-based cooling method following ideas from Figueiredo and Nowak (2003) and Elad et al. (2005). SEISMIC DATA RECOVERY Seismic 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 the Curvelet-based seismic processing same average sampling interval, random subsampling leads to an easily denoiseable spectrum, as long as the to-be-recovered signal 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 by Vermeer (1998). What differs in our approach is the nonlinear recovery from incomplete data. Recovery from compressive sampling depends on sparsity and mixing. Mixing turns harm- ful aliasing into relatively harmless noise and depends on the randomness 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 data regularization involves the solution of Pε with A := RCT , S := C, given incomplete data, y = Rf, with f the fully sampled data and R the picking matrix. In recent years, the authors repeatedly reported on succesful curvelet-based recovery of seismic data (see e.g. Herrmann, 2005; Hennenfent and Herr- mann, 2006a, 2007). Compared to other methods, such as sparse Fourier recovery (Sacchi and Ulrych, 1996; Zwartjes and 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 a 5m grid to a grid of 2.5m, where the groundroll is no longer aliased. Large amplitude contrast and regular subsampling make this example difficult. Focused recovery: Combining the non-adaptive curvelet trans- formwith the data-adaptive focal transform (Berkhout and Ver- schuur, 2006), leads to a powerful formulation where data is focused 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, ∆P, propagation paths that include the surface are re- moved, yielding a more focused wavefield and hence a more compressed curvelet vector. This improved focussing is achieved by Pε with the synthesis matrix A := R∆PCT and inverse spar- sity transform ST := ∆PCT . The solution of Pε now entails the inversion of ∆P, 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 RECOVERY After processing seismic data can be modeled by the linearized Born scattering/demigration operator, d = Km+ n, with the remaining 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. This image, y = KTd or y = Ψm+ e, contains the imprint of clut- ter, e = KTn, and the Gramm (demigration-migration) opera- tor, Ψ := KTK , and serves as input to our amplitude recovery scheme which estimates the reflectivity, m. The following ap- proximate identity is used AAT r'Ψr (3) with r a reference vector. This expression follows from a curvelet decomposition of the Gramm operator,Ψr'CTDΨCr. 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 estimated by solving a least-squares problem that uses the reference and demigrated-migrated reference vector as input. The factoriza- tion in Eq. 3 is based on a synthesis matrix given by A :=CTΓ with Γ := √ DΨ. This factorization leads to y' Ax0+ e as an approximate image representation that is amenable to Pε . Af- ter solving for x0, the reflectivity is obtained by applying the synthesis operator, ST := ` AT ´†, with † the pseudoinverse, toex. 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 Pε , leads to a stable (under noise) recovery of the imaged reflec- tivity. These improvements (cf. Fig. 3(a) and 3(b)) are obtained by a penalty functional that jointly promotes the curvelet spar- sity and the continuity along the reflectors. Results for the SEG AA’ dataset (O’Brien and Gray, 1997; Aminzadeh et al., 1997) are summarized in Fig. 3. These results are obtained for linearized Born data, modeled with the adjoint of reverse-time migration with optimal checkpointing (W. W. Symes, personal communication, 2007). The estimated images show a nice am- plitude recovery and clutter removal for reflection data with a signal-to-noise ratio (SNR) of only 3dB. DISCUSSION AND CONCLUSIONS The seismic data regularization derived its performance from curvelet compression, random subsampling and the inclusion of an additional focusing operator. The approximate invariance under demigration-migration allowed for a computationally- efficient migration amplitude recovery scheme. Since our scheme is based on a phase-space approximation, our scaling is less restrictive 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 their singular wavefront detection capability, curvelets represent in our vision an ideal domain for future developments in explo- ration seismology. ACKNOWLEDGMENTS The authors would like to thank Eric Verschuur and Chris Stolk for their input. We also would like to thank the authors of CurveLab for making their codes available andWilliam Symes for his reverse-time migration code. The examples presented were prepared withMadagascar (rsf.sourceforge.net/), supplemented by SLIMPy operator overloading, developed by Sean Ross Ross. ExxonMobil and Norsk Hydro are thanked for making the field datasets available. Eric Verschuur and M. O’Brien, S. Gray and J. Dellinger are thanked for providing the synthetic multiple dataset for the SEG AA’ available. This work was in part financially supported by the NSERC Dis- covery (22R81254) and CRD Grants DNOISE (334810-05) of F.J.H. and was carried out as part of the SINBAD project with support, 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 an interval 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 real SAGA 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) Image after nonlinear recovery. The clearly visible non-stationary noise in (a) is mostly removed during the recovery while the amplitudes are also restored. Steeply dipping reflectors under the salt are also well recovered. Curvelet-based seismic processing REFERENCES Aminzadeh, 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ès, E., J. Romberg, and T. Tao, 2006, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Candès, 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 SEG International 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. and Eng., 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és, 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