Curvelet-Based Migration AmplitudeRecoverybyPeyman P. MoghaddamB.Sc. in Electrical Engineering, Amirkabir University of Technology, Tehran, Iran,1995M.Sc. in Electrical Engineering, Amirkabir University of Technology, Tehran, Iran,1998A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate Studies(Geophysics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2010c©Peyman P. Moghaddam 2010AbstractMigration can accurately locate reflectors in the earth but in most casesfailstocorrectlyresolvetheiramplitude. Thismightleadtomis-interpretationof the nature of reflector.In this thesis, I introduced a method to accurately recover the ampli-tude of the seismic reflector. This method relies on a new transform-basedrecovery that exploits the expression of seismic images by the recently devel-oped curvelet transform. The elements of this transform, called curvelets,are multi-dimensional, multi-scale, and multi-directional. They also remainapproximately invariant under the imaging operator.I exploit these properties of the curvelets to introduce a method calledCurvelet Match Filtering (CMF) for recovering the seismic amplitude inpresence of noise in both migrated image and data.I detail the method and illustrate its performance on synthetic dataset. Ialso extend CMF formulation to other geophysical applications and presentresults on multiple removal. In addition of that, I investigate preconditioningof the migration which results to rapid convergence rate of the iterativemethod using migration.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiiStatement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xxiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Theme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7iiiTable of Contents2 Sparsity and continuity promoting seismic image recoverywith curvelet frames . . . . . . . . . . . . . . . . . . . . . . . . 92.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.1.1 Existing approaches . . . . . . . . . . . . . . . . . . . 112.1.2 Our approach . . . . . . . . . . . . . . . . . . . . . . 122.1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Approximation of the normal operator . . . . . . . . . . . . 182.2.1 Zero-order imaging operators . . . . . . . . . . . . . . 192.2.2 Curvelet frames . . . . . . . . . . . . . . . . . . . . . 202.2.3 Diagonal approximation of PsDO’s . . . . . . . . . . 222.2.4 Decomposition of the normal operator . . . . . . . . 252.2.5 Estimation of the diagonal . . . . . . . . . . . . . . . 262.3 Stable seismic image recovery . . . . . . . . . . . . . . . . . . 312.3.1 Curvelet frames for seismic images . . . . . . . . . . . 312.3.2 Stable recovery . . . . . . . . . . . . . . . . . . . . . 352.3.3 Solution by iterative thresholding . . . . . . . . . . . 372.3.4 Stable seismic recovery . . . . . . . . . . . . . . . . . 382.4 Sparsity- and continuity-enhanced seismic image recovery . . 402.5 Practical considerations . . . . . . . . . . . . . . . . . . . . . 442.6 Numerical experiments . . . . . . . . . . . . . . . . . . . . . 452.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 562.7.1 Initial findings . . . . . . . . . . . . . . . . . . . . . . 562.7.2 Extension to prestack imaging . . . . . . . . . . . . . 592.7.3 Recent related work . . . . . . . . . . . . . . . . . . . 602.7.4 Choice of the reference vector . . . . . . . . . . . . . 60ivTable of ContentsBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623 Curvelet-based migration preconditioning and scaling . . 673.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 673.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . 693.3 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . 723.3.1 Left preconditioning by fractional differentiation . . . 733.3.2 Right preconditioning by scaling in the physical do-main . . . . . . . . . . . . . . . . . . . . . . . . . . . 733.3.3 Right preconditioning by scaling in the curvelet do-main . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.4 Convergence of least-squares migration . . . . . . . . . . . . 763.5 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 783.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 803.7 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 81Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864 Curvelet-based seismic data processing . . . . . . . . . . . . 894.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 894.2 Curvelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 904.3 Common problem formulation by sparsity-promoting inver-sion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.4 Seismic data recovery . . . . . . . . . . . . . . . . . . . . . . 934.4.1 Curvelet-based recovery . . . . . . . . . . . . . . . . . 934.4.2 Focused recovery . . . . . . . . . . . . . . . . . . . . . 94vTable of Contents4.5 Seismic signal separation . . . . . . . . . . . . . . . . . . . . 964.6 Migration-amplitude recovery . . . . . . . . . . . . . . . . . 994.7 Discussion and conclusions . . . . . . . . . . . . . . . . . . . 1004.8 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 102Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045 True amplitude depth migration using curvelets . . . . . . 1085.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1085.2 Outline of the chapter . . . . . . . . . . . . . . . . . . . . . . 1105.3 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . 1105.4 Diagonal Approximation of the Normal Operator . . . . . . 1145.5 Curvelets and their invariance under normal operator . . . . 1155.6 Diagonal approximation . . . . . . . . . . . . . . . . . . . . . 1175.7 Practical Workflow . . . . . . . . . . . . . . . . . . . . . . . 1195.8 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.8.1 Noise-free case . . . . . . . . . . . . . . . . . . . . . . 1215.8.2 Noisy data . . . . . . . . . . . . . . . . . . . . . . . . 1305.9 Tolerance analysis . . . . . . . . . . . . . . . . . . . . . . . . 1365.10 Verification of approximation . . . . . . . . . . . . . . . . . . 1385.11 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 1385.12 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1405.13 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 140Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142viTable of Contents6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1466.1 Main contributions . . . . . . . . . . . . . . . . . . . . . . . 1466.1.1 Sparsity and continuity promoting seismic image re-covery with curvelet frames . . . . . . . . . . . . . . . 1476.1.2 Curvelet-based migration preconditioning and scaling 1486.1.3 Curvelet-based seismic data processing . . . . . . . . 1496.1.4 True amplitude depth migration using curvelets . . . 1496.2 Follow-up work . . . . . . . . . . . . . . . . . . . . . . . . . . 1506.2.1 Full waveform inversion . . . . . . . . . . . . . . . . . 1506.2.2 3D true amplitude migration . . . . . . . . . . . . . . 1506.3 Adding more priori knowledge to solution . . . . . . . . . . . 1506.4 Current limitations . . . . . . . . . . . . . . . . . . . . . . . 1516.4.1 Curvelet code . . . . . . . . . . . . . . . . . . . . . . 151Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153AppendicesA Proofs of lemma 1 and theorem 1 . . . . . . . . . . . . . . . 154B Reverse time wave equation migration . . . . . . . . . . . . 158B.1 Single scattering . . . . . . . . . . . . . . . . . . . . . . . . . 158B.2 Shot-geophone modeling and migration . . . . . . . . . . . . 159Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161viiList of Tables2.1 Algorithm for the estimation of the diagonal via regularizedleast-squares inversion. The Lagrange multiplier, η, is in-creased up to the point that all entries in the vector for thediagonal are positive. . . . . . . . . . . . . . . . . . . . . . . . 302.2 Sparsity-and continuity-enhancing recovery of seismic ampli-tudes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44viiiList of Figures2.1 Example of the recovery from different percentages of curveletcoefficients. The data set concerns a migrated image derivedfrom the Mobil data set, (a) near perfect recovery from p =99% of the coefficients, (b) recovery from only p = 3% of thecoefficients, (c) the difference between near perfect recovery(a) and approximate recovery (b). This difference does notcontain substantial coherent energy. . . . . . . . . . . . . . . 16ixList of Figures2.2 Invariance of curvelets under the discretized normal operatorΨ for a smoothly varying background model (a so-called lensmodel see Fig. 2.4(a)). Three coarse-scale curvelets in thephysical domain before (a) and after application of the nor-mal operator (b) in the physical (a-b) and Fourier domain(e-f). The results for three fine-scale curvelets are plotted in(c-d) for the physical domain and in (g-h) for the Fourierdomain. Remark: The curvelets remain close to invariantunder the normal operator, a statement which becomes moreaccurate for finer scale which is consistent with Theorem 1.The example also shows that this statement only holds forcurvelets that are in the support of the imaging operator ex-cluding steeply dipping curvelets. . . . . . . . . . . . . . . . . 172.3 Spatial and frequency representation of curvelets, (a)four dif-ferent curvelets in the spatial domain at three different scales,(b) dyadic partitioning in the frequency domain, where eachwedge corresponds to the frequency support of a curvelet inthe spatial domain. This figure illustrates the microlocal cor-respondence between curvelets in the physical and Fourier do-main. Curvelets are characterized by rapid decay in the physi-cal space and of compact support in the Fourier space. Noticethe correspondence between the orientation of curvelets in thetwo domains. The 90◦ rotation is a property of the Fouriertransform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22xList of Figures2.4 Stylized example of the diagonal estimation for a smooth ve-locity model and a reflectivity consisting of three reflectionevents including a fault, (a) the smooth background modelwith the low-velocity lens. This velocity model is used to cal-culate the linearized scattering and migration operators, (b)the bandwidth limited reflectivity, (c) the imaged reflectiv-ity from data consisting of 32 shot records with 500 traceseach, (d) the approximated normal operator for diag(tildewideu) es-timated with η = 1. The bandwidth limited reflectivity in(b) and the image in (c) served as input for the referenceand ’data’ vectors for the diagonal estimation procedure out-lined in Table 2.1. Remarks: the main dimming of the normaloperator is captured by the diagonal approximation. The rel-ative lscript2-error between the actual and the approximate normaloperator is 6.1%. . . . . . . . . . . . . . . . . . . . . . . . . . 322.5 Estimates for the diagonal tildewideu are plotted in (a-d) for increas-ing η ={0.01, 0.1, 1, 10}. The diagonal is estimated accord-ing the procedure outlined in Table 2.1 with the reference and’data’ vectors, v and b, plotted in Fig. 2.4(b) and 2.4(c). Asexpected the diagonal becomes more positive for increasing η. 33xiList of Figures2.6 Decay rate for the curvelet coefficients compared to the de-cay for the coefficients of Fourier (dashed line) and discretewavelet transforms (dot-dashed line), (a) decay rate for thesorted curvelet coefficients (solid line) of the real migrated im-age included in Fig. 2.1(a), (b) the same as (a) but now forthe synthetic reflectivity of the SEG/EAGE AAprime salt model(cf. Fig. 2.8(a)). The decay for the curvelet transform clearlycompares favorably to these other two transforms. . . . . . . 352.7 The SEG/EAGE AAprime salt model (Aminzadeh et al., 1997),(a) The original velocity model, (b) the smoothed velocitymodel with a window size of 720×720m. Remark: Imagingthis model is a challenge because of the presence of the high-velocity (bright) salt. . . . . . . . . . . . . . . . . . . . . . . . 472.8 ImagingoftheSEG/EAGEAAprime saltmodel(Aminzadehetal.,1997), (a) reflectivity defined in terms of the band-pass fil-tered difference between the smoothed velocity model and aslightly less smoothed velocity model, (b) the imaged reflec-tivity according to Eq. 2.24. The main reflection events arepresent but suffer from deteriorated amplitudes, especiallyunder the high-velocity salt and for steep reflectors and faults. 48xiiList of Figures2.9 Images that are input to the diagonal approximation schemeoutlined in Table 2.1, (a) the reference vector, r, derivedfrom Fig. 2.8(b), after applying corrections for the spherical-divergenceandreceiver-arraytaper, (b)thedemigrated-migratedreference vector, b = Ψr, that serves as ’data’ for the diag-onal estimation, (c) diagonal approximation on the originalbandwidth-limited reflectivity plotted in (a). This diagonalapproximationbyAATmcapturesthenormaloperator, Ψm,quite well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502.10 Gradient vectors for the reference vector r plotted in Fig.2.9(a). The gradients (white ↑’s) are used to calculate thetangential directions along which the additional anisotropicsmoothing is applied. . . . . . . . . . . . . . . . . . . . . . . . 512.11 Images after nonlinear recovery (Pepsilon1). (a) result with the lscript1-norm recovery only (β = 0); (b) recovery with the combinedlscript1 and continuity recovery for α = β = 1/2. The amplitudesare recovered. The anisotropic diffusion successfully removesthe artifacts. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53xiiiList of Figures2.12 Recovery of the amplitudes after normalization with respectto the first event, (a) seismic traces of the original reflectivity(dot-dashed line), the migrated reflectivity (dashed line) andthe recovered reflectivity (solid line), (b) amplitudes alongthe bottom most horizontal reflector, (c) average amplitudespectra of the depth-dependent reflectivity. The original andrecovered spectra are normalized to match. Remarks: Thenonlinear recovery corrects for the amplitudes and restoresthe spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . . . 542.13 Image amplitude recovery for noisy data (SNR 3dB), (a)noise image according to Eq. 2.24, (b) image after nonlin-ear recovery from noisy data (Pepsilon1). The clearly visible non-stationary noise in (a) is removed during the recovery whilethe amplitudes are also restored. Steeply dipping reflectorsdenoted by the arrows are also well recovered. . . . . . . . . 572.14 ’Denoising’ of a shot record, (a) the noise-free data receivedby the receiver array, (b) noisy data with a SNR of 3dB, (c)forward modeled data after amplitude recovery. Observe thesignificant improvement in the data quality, reflected in anincrease for the SNR of 19.2dB. . . . . . . . . . . . . . . . . 58xivList of Figures3.1 The SEG/EAGE AAprime salt model, (a) reflectivity defined bythe high-pass filtered velocity model, (b) smoothed velocitymodel, (c) the migrated image according to Equation 3.1.This image suffers from deteriorated amplitudes, especiallyunder the high-velocity salt and for steep reflectors and faults. 823.2 Migrated images for different levels of preconditioning, (a)result for left preconditioning (level I, cf. Equation 3.6), (b)result for left-right (including depth-correction) precondition-ing (level II, cf. Equation 3.7), (c)the same but now includingcurvelet-domain scaling (level III, cf. Equation 3.10). . . . . . 833.3 Residual decays for different levels of preconditioning. Thedotted blue lines corresponds to least-squares migration with-out preconditioning, the dash-dotted lines to level I precon-ditioning, the dashed black lines to level II preconditioning,and the red solid lines to level III preconditioning. This isoffset by one iteration to account for the overhead, (a) plotfor the decay of the data-space normalized residues µk as afunction of the number of LSQR iterations, (b) the same butnow for the model-space normalized residuals νk. . . . . . . . 843.4 Least-squares migration without and with preconditioningjuxtaposed with the original reflectivity (in light blue), (a)least-squares migrated image, (b) least-squares image withlevel III preconditioning. Notice the improvement in the re-covered reflectivity. . . . . . . . . . . . . . . . . . . . . . . . . 85xvList of Figures4.1 Comparison between 3-D curvelet-based recovery by sparsity-promoting inversion with and without focusing, (a)fully sam-pled real SAGA data shot gather, (b) randomly subsampledshot gather from a 3-D data volume with 80% of the tracesmissing in the receiver and shot directions, (c) curvelet-basedrecovery, (d) curvelet-based recovery with focusing. Noticethe improvement (denoted by the arrows) from the focusingwith the primary operator. . . . . . . . . . . . . . . . . . . . 954.2 3-DPrimary-multipleseparationwithPepsilon1 fortheSAGAdataset,(a) near-offset section including multiples, (b) the SRME-predicted multiples, (c) the estimated primaries accordingto lscript2-matched filtering, (d) the estimated primaries obtainedwith Pepsilon1. Notice the improvement, in areas with small 3-Deffects (ellipsoid) and residual multiples. . . . . . . . . . . . . 984.3 Image amplitude recovery for a migrated image calculatedfrom noisy data (SNR 3dB), (a) image with clutter, (b) im-ageafternonlinearrecovery. Theclearlyvisiblenon-stationarynoise in (a) is mostly removed during the recovery while theamplitudes are also restored. Steeply dipping reflectors (de-noted by the arrows) under the salt are also well recovered. . 101xviList of Figures5.1 Spatial and frequency representation of curvelets, (a)four dif-ferent curvelets in the spatial domain at three different scales,(b) dyadic partitioning in the frequency domain, where eachwedge corresponds to the frequency support of a curvelet inthe spatial domain. This figure illustrates the microlocal cor-respondence between curvelets in the physical and Fourier do-main. Curvelets are characterized by rapid decay in the phys-ical space and of compact support in the Fourier space. No-tice the correspondence between the orientations of curveletsin the two domains. The 90o rotation is a property of theFourier transform. Courtesy of Herrmann et al. (2008). . . . 1155.2 Conflicting dips velocity model and reflectivity (a) velocitymodel, (b) smooth velocity model used to generate our exam-ple, (c) reflectivity generated by subtracting smooth velocitymodel from original one . . . . . . . . . . . . . . . . . . . . . 1225.3 Diagonal approximation of normal operator, (a)reference im-age (depth corrected migrated image), (b) normal operatoron the reference image (i.e., Ψr) (c) approximate normal op-erator on the reference image (i.e., CTWCr), approximationerror is %3.71 . . . . . . . . . . . . . . . . . . . . . . . . . . 124xviiList of Figures5.4 Curvelet-domainrepresentationofthecurveletvectorobtainedfrom (a) reference vector, (b) scaling coefficients withoutsmoothing constraints, (c), scaling coefficients with smooth-ingconstraint. Thedifferentsub-imagesrepresentthecurveletcoefficients at different scales (coarsest in the center) and dif-ferent angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1265.5 Examples of the solution for optimization problem QP withsoft-thresholding, (a) λ = 0, (b) λ = σb (c) λ = 2σb . . . . 1275.6 Examples of the solution for optimization problem BPDNwith SPGL1, (a) epsilon1 = σb, (b) epsilon1 = 50σb (c) epsilon1 = 100σb . . . . 1295.7 Reflectivity model, the vertical lines are locations of investi-gation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1305.8 Amplitude analysis of the recovered images using proposedmethods, (a) offset=3480 (m) , (b) offset=5160 (m), (c) off-set=5880 (m) . . . . . . . . . . . . . . . . . . . . . . . . . . . 1315.9 Amplitude analysis of the recovered images using proposedmethods, (a) offset= 8640 (m), (b) offset=10080 (m) , (c)offset=12720 (m), note the salt bottom amplitude mismatch 1325.10 Syntheticdata, (a)before, (b)afteraddingthenoise, SNR =10dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1335.11 Depth corrected migrated noisy image . . . . . . . . . . . . . 1345.12 Examples of the solution for optimization problem QP withsoft-thresholding, (a) λ = σb, (b) λ = 2σb (c) λ = 3σb . . . 1355.13 Examples of the solution for optimization problem BPDNwith SPGL1, (a) epsilon1 = 10σb, (b) epsilon1 = 140σb (c) epsilon1 = 200σb . . 137xviiiList of Figures5.14 Comparison between the approximation terms in Equation5.21, (a) left hand term El = CTW−12CCTWx , (b) righthand term Er = CTW12x, relative error is 4.09 percent. . . . 139xixPrefaceThis thesis was prepared with Madagascar, a reproducible research soft-ware package available at rsf.sf.net, in such a way that most of the repro-ducible results are linked to the code that generated them. Reproducibilityfacilitates the dissemination of knowledge not only within the Seismic Lab-oratory for Imaging and Modeling (SLIM) but also between SLIM and itssponsors, and more generally, the entire research community.The programs required to reproduce the reported results are Mada-gascar programs written in C/C++, MatlabR©, or Python. The numeri-cal algorithms and applications are mainly written in Python as part ofSLIMpy(slim.eos.ubc.ca/SLIMpy)withafewexceptionswritteninMatlabR©or Python.xxAcknowledgmentsI am greatly indebted to my advisor Felix Herrmann. A few words ofacknowledgment are just not enough to express my gratitude and appreci-ation for all he has done for me. I also wish to thank him for his kind andsometimes (very) energetic encouragements.I would like to express my gratitude to Reza Shahidi, Deli Wang, CodyBrown, Gilles Hennenfent, and the remainder of the SLIM team for theirhelp and friendships. I was very fortunate to be part of this group.The next persons who have influenced my scientific life are Chris Stolk,William Symes, Uri Ascher and Chen Grief. They helped my research in-cursions in the areas of applied mathematics and optimization. I would liketo thank them for their time and patience.My appreciation goes to Total E&P for the hospitality during my intern-ship. In particular, I wish to thank Eric Dussaud, and Paul Williamson forgreat technical discussions.I would like to express my gratitude to Chen Grief, Michael Bostock,Doug Oldenburg for serving on my supervisory committee.Finally, many thanks to my family for their love, support and encour-agement throughout my studies.xxiTo my familyxxiiStatement of Co-AuthorshipChapter 2 was published jointly with Felix J. Herrmann and Chris Stolkin the mathematical literature where alphabetical ordering of authors iscustomary. Aside from providing the numerical examples for this paper, Ico-authored the text for the examples and I wrote the description of thepart that involved our nonlinear solver.Chapter 3 was published with Felix J. Herrmann, Cody R. Brown, andYogi A. Erlangga. My contribution consisted of the development of thecurvelet matched filter and the design of the imaging experiment. Codyassisted with generating the results and Felix wrote most of the manuscriptwith significant inputs from me and the other co-authors.Chapter 4 is a review paper highlighting three distinctively different ap-plications of the curvelet transform to seismic data processing and imaging.I was responsible for the imaging example and I authored the correspondingtext.Chapter 5 was prepared for submission for publication jointly with Fe-lix J. Herrmann with input from Reza Shahidi. I was responsible for allexamples and I authored the text.xxiiiChapter 1IntroductionExploration seismology is a complex and costly operation. On land,dynamiteorVibroseissources canbe usedtosendenergyintothesubsurface.This energy propagates and partially reflects at interfaces due to a changein rock properties. The reflected wavefield is recorded at the earth’s surfaceby an array of geophones. At sea, the principle remains the same but theseismic source is typically an air gun and the receivers are hydrophones onstreamer lines towed by a seismic vessel. Seismic processing is a techniquewidely used by the oil industry to explore and identify potential oil-richareas into the earth using the recorded data.Seismic migration is vital part of seismic processing. Migration generatesan image of the subsurface that is finally interpreted by geo-scientists.Migration difficulties generally arise from assumptions made in algo-rithms, that are not met by underlying theory. In particular, modern mi-gration algorithms can accurately position reflecting features in the Earth’ssubsurface, however they generally do not correctly resolve the amplitudes.This problem has different causes. First, the migration operator is the ad-joint but not the inverse of the linearized Born modeling or scattering op-erator (see e.g. Guitton, 2004; Claerbout, 1985; Gray, 1997; Symes, 2008).11.1. ThemeSecond, There are also other sources contributing to this inaccuracy in-cluding geometrical spreading of wavefront energy, acquisition limitations,velocity heterogeneity and anomalous focusing.There exists a wide variety of migration amplitude recovery techniques:• iterative-based methods use migration and scattering operatorsand iteratively recover the amplitude of the reflectors through a Krylovsubspace-based method such as LSQR or conjugate-gradient (Kuehland Sacchi, 2003). They typically fairly computationally intensive• match filter-based methods represent another type of amplituderecovery approach that include the migrated and re-migrated image.They require ”only one” time modeling and migration of the imageand computationally cheaper than iterative based methods. (Guitton,2004; Claerbout, 1985; Rickett, 2003; Symes, 2008) are examples ofthis approach.1.1 ThemeThe method proposed in this thesis follows the match filter-based ap-proach. The main theme of this thesis is a practical, robust, and geomet-rical—i.e., transform-based—approach to the seismic migration amplituderecovery problem. The motivation of this approach is that two key featuresof seismic data are, in my opinion, not used to their full extent in existingapproaches, namely• Representation of the seismic normal operator in transform21.1. Themedomain The seismic normal operator typically belongs to a class ofoperators called Pseudo-differential operators. These operators haveinteresting properties to make them suitable to represent in the trans-form domain.• Strong geometrical structureSeismic images are a spatio-temporalsampling of the reflectors in the Earth’s subsurface. These reflectorsare mathematically continuous functions of reduced dimension (i.e.,curves in 2D or sheets in 3D)To make the most of these properties, my approach uses the curvelettransform (Cand`es and Donoho, 2004) which is data-independent, multi-scale, and multidirectional. The elements of this transform, the curvelets,are localized in the frequency domain and of rapid decay in the physicaldomain. They are very efficient at representing curve-like singularities—e.g., reflectors. In other words, only a few curvelets are needed to representthe complexity of real seismic image. I use this piece of information, calledsparsity, to help solve the imaging problem. The idea of sparsity-promotinginversion is in itself not new to geophysics.The other interesting property of the curvelets is their behavior underthe action of acoustic wave propagator (e.g., Cand`es and Demanet, 2003;Smith, 1998, 1997) . A theorem in chapter 2 proved that the curvelets areapproximately invariant under the action of the normal operator . Further-more I adopt the result of this theorem to introduce match-filtering approachfor estimation of this operator.31.2. Objectives1.2 ObjectivesThe objectives of this thesis are two-fold:• develop an in-depth understanding of successful curvelet match filter-ing for the approximation of normal operator and its key ingredients.• formulate a practical sparsity-promoting seismic migration amplituderecovery algorithm.1.3 OutlineIn chapter 2, I observe that curvelets, as a directional frame expansion,lead to sparsity of seismic images and exhibit invariance under the normaloperator of the linearized imaging problem. Based on this observation Iderive a method for stable recovery of the migration amplitudes from noisydata. The method corrects the amplitudes during a post-processing stepafter migration, such that the main additional cost is one application of thenormal operator, i.e. a modeling followed by a migration. Asymptoticallythis normal operator belongs to a class of operators known as pseudodiffer-ential operators, for which a diagonal approximation in the curvelet domainis derived, including a bound for its error and a method for the estimation ofthe diagonal from a compound operator consisting of discrete implementa-tions for the scattering operator and its adjoint the migration operator. Thesolution is formulated as a nonlinear optimization problem where sparsityin the curvelet domain, as well as continuity along the imaged reflectors,are jointly promoted. To enhance sparsity, the lscript1-norm on the curvelet41.3. Outlinecoefficients is minimized, while continuity is promoted by minimizing ananisotropic diffusion norm on the image. The performance of the recoveryscheme is evaluated with a time-reversed ’wave-equation’ migration code onsynthetic datasets, including the complex SEG/EAGE AAprime salt model.In chapter 3, I introduce a preconditioner for the inversion of the lin-earized Born scattering operator. This preconditioner approximately cor-rects for the “square root” of the normal operator. This approach consistsof three parts, namely (i) a left preconditioner, defined by a fractional timeintegration designed to make the migration operator zero order, and tworight preconditioners that apply (ii) a scaling in the physical domain ac-counting for a spherical spreading, and (iii) a curvelet-domain scaling thatcorrects for spatial and reflector-dip dependent amplitude errors. I showthat a combination of these preconditioners lead to a significant improve-ment of the convergence for iterative least-squares solutions to the seismicimaging problem based on reverse-time migration operators.In chapter 4, I show that other geophysical problems—e.g., focused re-covery, seismic signal separation, and migration amplitude recovery—can bere-cast in the formulation used for curvelet sparsity inversion formulation.This puts in a broader perspective the insights gained during the develop-ment of curvelet sparsity inversion formulation.In chapter 5, I focused on the practical approach of amplitude recoverymethod. I revisit the amplitude recovery method and extend it threefold.First, I replace the linear least-squares formulation for the estimation ofthe curvelet-domain coefficients (Herrmann et al. (2008)) by a non-linearleast-squares formulation that imposes positivity on the scaling coefficients.51.3. OutlineThis comes from the fact that the normal operator is positive semi-definite.Second, I consider the noise term both in the data and that generated bymigration due to incomplete data (i.e., few shots). Third, I investigate twomethodsforrecoveryandcomparethemintermsofqualityandperformance,a soft thresholding-correction technique and a sparsity promotion technique.In Chapter 6, I summarize the work done in this thesis, and discuss someof its aspects in a broader context. Conclusions and recommendations forfuture research follow.Appendix A contains further details about the curvelet transform. Itproves the theorem and associated lemma regarding the invariance of thecurvelet elements under the action of normal operator. This appendix pairswith chapter 2.6BibliographyCand`es, E. J. and L. Demanet, 2003, Curvelets and Fourier integral opera-tors: Compte Rendus de l’Acad´emie des Sciences, Paris, 336, no. 1, 395– 398.Cand`es, E. J. and D. L. Donoho, 2004, New tight frames of curvelets andoptimal representations of objects with C2 singularities: Communicationson Pure and Applied Mathematics, 57, no. 2, 219 – 266.Claerbout, J., 1985, Earth soundings analysis, processing versus inversion:Blackwell Scientific Publications.Gray, S. H., 1997, True amplitude seismic migration: A comparison of threeapproaches: Geophysics, 62, 929–936.Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor nonunitary imaging operators: Geophysics, 69, 1017–1024.Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2008, Sparsity- andcontinuity-promoting seismic image recovery with curvelet frames: Ap-plied and Computational Harmonic Analysis, 24, 150–173.Kuehl, H. and M. Sacchi, 2003, Least-squares wave-equation migration forAVP/AVA inversion: Geophysics, 68, 262–273.Rickett, J., 2003, Illumination based normalization for wave-equation depthmigration: Geophysics, 68, 208–221.7Chapter 1. BibliographySmith, H., 1997, A Hardy space for Fourier integral operaotrs: Geom. Anal.,7, 784–842.——–, 1998, A parametrix construction for wave equations with C1,1 coef-ficients: , Ann. Inst. Fourier, Grenoble, 48, 797–835.Symes, W., 2008, Optimal scaling for reverse time migration: submitted toGeophysics.8Chapter 2Sparsity and continuitypromoting seismic imagerecovery with curveletframes2.1 IntroductionThis paper introduces a general-purpose algorithm for the stable recoveryof the amplitudes in seismic images. The method involves the inversion indimension d of a zero-order pseudodifferential operator (PsDO) of the typeparenleftbigΨfparenrightbig(x) = integraldisplayRde−ix·ξa(x,ξ)ˆf(ξ)dξ (2.1)withˆf(ξ) = 1(2pi)dintegraldisplayRdf(x)eix·ξdx (2.2)A version of this chapter has been published. Herrman, F.J., Moghaddam, P.P. andStolk, C. (2008) Sparsity- and continuity-promoting seismic image recovery with curveletframes. Applied and Computational Harmonic Analysis, Vol. 24, No. 2, pp. 150-173,2008.92.1. Introductionthe Fourier transform with x, ξ∈Rd the spatial coordinate and wavenum-ber vectors and a(x,ξ) the symbol of the pseudodifferential operator. Undercertain conditions that include no grazing rays and high-frequency asymp-totics, the above linear operator describes seismic data after imaging (tenKroode et al., 1998; de Hoop and Brandsberg-Dahl, 2000; Stolk and Symes,2003). A seismic image is derived from noisy data given byd(xs,xr,t) = parenleftbigKmparenrightbig(xs,xr,t)+n(xs,xr,t) (2.3)with K the Born scattering operator, m(x) the model with the reflectivityand n zero-mean standard deviation σ white Gaussian noise. The seismicdata volume is a function of the source and receiver positions, xs ∈ Rd−1and xr ∈Rd−1 and of time, t∈R+0 . After applying the migration operatorKT, (the symbol T is reserved for the transpose), to the data volume an“image” is obtained according to(KTd)(x) = parenleftbigKTKmparenrightbig(x)+parenleftbigKTparenrightbign(x)y(x) = parenleftbigΨmparenrightbig(x)+e(x). with x∈Rd (2.4)This equation for the image y corresponds to the normal equation and con-tains contributions from the normal operator, Ψ, as well as from imagednoise e. The operator K preserves the singularities in y, i.e., KTKm≈Id mwith Id the identity matrix and the symbol≈indicating ”approximated by”in the locations of the singularities. With the increased demand for high-quality images, the above approximation of Ψ≈Id is no longer justifiable102.1. Introductionbecause it may lead to errors in the relative amplitudes of the imaged re-flectors. Our main goal is to derive a method that recovers these imagedamplitudes, i.e. estimate models m from “images” y with clutter.2.1.1 Existing approachesWith the increasing demand for hydrocarbons and the increasingly hardto image (subsalt), the recovery of correct seismic amplitudes has becomemore necessary. Approaches in the extensive literature on this topic rangefrom least-squares migration, where the normal or Gramm matrix is invertedwith Krylov subspace methods (see e.g. Nemeth. et al., 1999; Chavent andPlessix, 1999; de Hoop and Brandsberg-Dahl, 2000), to high-frequency (mi-crolocal)methods(seee.g. tenKroodeetal.,1998;deHoopandBrandsberg-Dahl, 2000; Stolk and Symes, 2003), where the normal operator is invertedbased on ray-asymptotic arguments, to methods that apply a diagonal scal-ing (see e.g. (Rickett, 2003; Guitton, 2004; Plessix and Mulder, 2004) andmore recently Symes (2006a)) to approximately invert the Gramm operatorof ’wave-equation’ migration. These scaling methods either calculate theweighting analytically (ten Kroode et al., 1998; Plessix and Mulder, 2004),based on certain assumptions regarding the acquisition and velocity model,or estimate the diagonal weighting from the action of the normal operatoron some reference vector, an idea first suggested by Symes and reportedin (Claerbout and Nichols, 1994). In this paper, a similar approach is fol-lowed, where the normal operator is replaced by a diagonal matrix actingon the curvelet coefficients of the image (Candes et al., 2006a; Hennenfentand Herrmann, 2006).112.1. Introduction2.1.2 Our approachThe computational cost of evaluating the migration operator isO(nd+2) fora d-dimensional image with n samples in each direction. This large costmakes it a major challenge to conduct least-squares migration for d > 2 orlarge n Nemeth. et al. (1999); Chavent and Plessix (1999); de Hoop andBrandsberg-Dahl (2000); Hu et al. (2001); Kuhl and Sacchi (2003); Plessixand Mulder (2004); Hu et al. (2001). We address this issue by exploit-ing recently developed curvelet frames. These frame expansions compressseismic images (see e.g. Candes et al., 2006a; Hennenfent and Herrmann,2006, for the compression of seismic data and Fig. 2.1 for the compressionof a typical migrated seismic image) and consist of a collection of frameelements ’curvelets’ that are approximately invariant under PsDO’s. Thesetwo properties allow for the development of an approach similar to the so-called wavelet-vaguelette method (WVD) (see e.g. Donoho, 1995; Candesand Donoho, 2000b; Lee and Lucier, 2001). In this approach, scale-invarianthomogeneous operators are nonlinearly inverted, using the eigenfunction-likebehavior of wavelets and curvelets.Our main contribution to earlier ideas on stable seismic image recovery(see e.g. Herrmann, 2003) is threefold: First, the WVD method is extendedto expanding PsDO’s with respect to redundant curvelet frames. Second, itis observed that order-zero pseudodifferential operators with homogeneousprincipal symbols act approximately as a multiplication on curvelets. Thisproperty implies that such operators can be approximated by a diagonalmatrix acting on curvelet coefficients of an image, followed by an inverse122.1. Introductioncurvelet transform. We refer to this procedure as a curvelet-domain diag-onal approximation (or weighting) that becomes more accurate for smallerscales (higher frequencies). Third, a formulation for the amplitude recov-ery is presented in terms of a nonlinear sparsity-and continuity-enhancingoptimization problem.After discretization, the nonlinear optimization problem for the seismicamplitude recovery (see e.g. Tsaig and Donoho, 2006; Candes et al., 2006b;Daubechies et al., 2005) has the following formPepsilon1 :tildewidex = argminxJ(x) = Js(x)+Jc(x) subject to bardbly−Axbardbl2≤epsilon1tildewidem = parenleftbigATparenrightbig† tildewidex.(2.5)During the optimization, the vector x is optimized with respect to thepenalty functional J(x) and the lscript2-data misfit. The penalty functional J(x)is combination of sparsity penalty function Js(x) and continuity penaltyfunction Jc(x). The term sparsity vector is used for x to stress the pointthat the magnitude-sorted coefficients of this vector are of rapid decay, byvirtue of compression by curvelet frames.The above optimization problem is solved for the model by finding acoefficient vector tildewidex that minimizes the penalty term subject to fitting thedata to within a noise-level dependent tolerance level epsilon1. The ’tilde’ symbol(tildewide) is reserved for vectors that solve a nonlinear optimization problem. Thepenalty functional J(x) is designed to promote the sparsity of the imagein the curvelet domain (i.e., Js(x)) as well as continuity along the imagedreflectors (i.e., Jc(x)). Bold lowercase symbols refer to discretized vectors132.1. Introductionand bold uppercase symbols to discretized operators. Non-boldface symbolsrefer to continuous infinite dimensional functions and operators.Estimates for the model tildewidem are calculated by applying the pseudo-inverse(denoted by the symbol†) of the analysis matrix to tildewidex. This analysis matrix isgiven by the adjoint of the synthesis matrix (AT). This synthesis is given bythe diagonally-weighted curvelet synthesis matrix with a weighting designedsuch thatAATrsimilarequalΨr. (2.6)In this expression, r represents an appropriately chosen discrete referencevector and Ψ the discrete normal operator, formed by compounding thediscrete scattering and its transpose the migration operator, i.e., Ψ := KTKwith the symbol := denoting ’defined as’. The symbolsimilarequalis used to denotethat this expression is approximate, a statement we will make more precisewith a bound on the error.After forming the normal operator with one’s favorite numerical imple-mentation for the migration operator and its adjoint, our amplitude recoverymethod consists of the following steps:1. Calculate an appropriate reference vector, r, from the data that isclose enough to the unknown image. Typically, this reference vector isdefined by a migrated image to which a simple amplitude correctionis applied;2. Estimate a diagonal curvelet-domain weighting that approximates thenormal operator on the reference vector. This diagonal defines thesynthesis matrix A;142.1. Introduction3. Estimate tildewidex by solving the nonlinear optimization problem Pepsilon1. Thisprogram inverts the synthesis matrix;4. Calculate the model tildewidem from the recovered coefficient vector tildewidex throughthe pseudo-inverse of AT.The proposed recovery method derives from two essential properties ofcurvelet frames, namely the ability of this signal representation to compressseismic images and the invariance of curvelets under the normal operator.Fig. 2.1 and 2.2 are included to stress the importance of these propertiesby confirming that a real migrated image can accurately be recovered fromonly 3% of the curvelet coefficients, while curvelets remain invariant underthe normal operator defined for a smooth background velocity model. Dur-ing this paper, we use the first property to formulate a sparsity promotingseismic image in term of a nonlinear optimization and the second propertyfor diagonal approximation of the normal operator in curvelet domain.2.1.3 OutlineFirst, a diagonal approximation of the normal operator in the curvelet do-main is derived. The error of this approximation is bounded, using proper-ties of the curvelet tiling of phase space. In this derivation, use is made ofthe property that the normal operator can be described as a zero-order pseu-dodifferential operator. Next, a method is proposed to estimate the diago-nal weighting matrix that uses the property that the symbol of the normaloperator is smooth in phase space and positive. This diagonal approxima-tion leads to an approximate seismic image representation by a diagonally-152.1. Introduction(a)(b)(c)Figure 2.1: Example of the recovery from different percentages of curveletcoefficients. The data set concerns a migrated image derived from the Mobildata set, (a) near perfect recovery from p = 99% of the coefficients, (b)recovery from only p = 3% of the coefficients, (c) the difference betweennear perfect recovery (a) and approximate recovery (b). This differencedoes not contain substantial coherent energy.162.1. Introduction(a) (b) (c) (d)(e) (f) (g) (h)Figure 2.2: Invariance of curvelets under the discretized normal operatorΨ for a smoothly varying background model (a so-called lens model seeFig. 2.4(a)). Three coarse-scale curvelets in the physical domain before(a) and after application of the normal operator (b) in the physical (a-b) and Fourier domain (e-f). The results for three fine-scale curvelets areplotted in(c-d)for the physical domain and in(g-h)for the Fourier domain.Remark: The curvelets remain close to invariant under the normal operator,a statement which becomes more accurate for finer scale which is consistentwith Theorem 1. The example also shows that this statement only holds forcurvelets that are in the support of the imaging operator excluding steeplydipping curvelets.172.2. Approximation of the normal operatorweighted curvelet transform. With this representation, the seismic imageamplitudes can be recovered by solving a nonlinear optimization problem.During this optimization problem, the image is recovered by minimizing thedata mismatch while promoting sparsity in the curvelet domain and con-tinuity along imaged reflectors. This paper is concluded by applying thealgorithm to a synthetic data set derived from the SEG AAprime salt model.2.2 Approximation of the normal operatorThe primary focus of seismic imaging is to locate singularities in the earth’selastic properties from seismic data recorded at the surface (Beylkin, 1984;de Hoop and Bleistein, 1997; ten Kroode et al., 1998; Stolk, 2000; Bleis-tein et al., 2001; Brandsberg-Dahl and de Hoop, 2003). A seismic surveyconsists of multiple seismic experiments in which both the location of thesources and receivers are varied along the surface. The acquired data is sub-sequently used to create images of the singularities in the subsurface. Themain purpose of the recovery method presented in this paper is to recoverthe relative amplitudes of seismic images from data that is possibly contam-inated with noise. For this purpose, a diagonal approximation of the normaloperator in the curvelet domain is presented. This approximation derivesfrom the invariance properties of curvelets under the normal operator (seeFig. 2.2). The approximation leads to a SVD-like decomposition for thenormal operator and makes the large-scale seismic image recovery problemamenable to a solution by nonlinear optimization.182.2. Approximation of the normal operator2.2.1 Zero-order imaging operatorsIn the high-frequency limit, the scattering operator and the normal oper-ator can, under certain conditions on the medium and ray-geometry, beconsidered as Fourier Integral Operators (FIO’s) (ten Kroode et al., 1998).In dimension two (d = 2), the scattering operator, K, and its adjoint themigration operator, KT, can both be considered as FIO’s of order 14, whilethe leading behavior for their composition, the normal operator Ψ, corre-sponds to an order-one invertible elliptic PsDO(ten Kroode et al., 1998;de Hoop and Brandsberg-Dahl, 2000; Stolk and Symes, 2003). To make thisPsDOamenable to an approximation by curvelets, the following substitu-tions are made for the scattering operator and the model: Kmapsto→K (−∆)−1/2and mmapsto→(−∆)1/2 m with ((−∆)αf)∧(ξ) =|ξ|2α· ˆf(ξ), with ∆ is the Lapla-cian operator. Alternatively, these operators can be made zero-order bycomposing the data side with a 1/2-order fractional integration along thetime coordinate, i.e., K mapsto→ ∂−1/2t K. After these substitutions, the normaloperator Ψ becomes zero-order. Similar substitutions are made in the WVDmethods where vaguelettes are introduced according to the same mappings.Our seismic image amplitude problem is different from the recovery ofimages in ill-posed problems such as inverting the Radon transform (Candesand Donoho, 2000b). In our case, the ill-posedness comes from small valuesof the symbol at certain positions and wave numbers and is not related tothe ill-posedness stemming from the inversion of a normal operator thatacts as a (fractional) inverse Laplacian. Our normal operator, without theabove substitutions, acts as a (fractional) Laplacian. This behavior makes192.2. Approximation of the normal operatorthe problem ill-posed for the low-frequencies. Since seismic data and theapproximations of the operators are both high-frequency, this ill-posednessis negated. The ill-posedness that is a concern are the entries in the symbolof the normal operator that correspond to regions in the model space thatare badly insonified.2.2.2 Curvelet framesCurveletsaredirectionalframesthatrepresentaspecifictilingofthetwo/three-dimensional frequency plane into multiscale and multi-angular wedges (seeFig. 5.1). Because the directional sampling increases every-other scale dou-bling, curvelets become more and more anisotropic at finer and finer scales.They become ’needle-like’ as illustrated in Fig. 5.1. Curvelets are localizedin Fourier space and their smoothness in this domain leads to a rapid decayin the physical domain. Their effective support is given by ellipsoids witha width ∝ 2j/2 and length ∝ 2j and an angle θjl = 2pil2floorleftj/2floorright with j thescale and l the angular index. Curvelets are indexed by the multi-indexµ := (j, l, k) ∈M with M the multi-index set running over all scales, j,angles, l, and positions k (see for details Candes et al., 2006a; Ying et al.,2005). Following Candes et al. (2006a), define the continuous curvelet familyfor x∈R2 asϕµ(x) = 23j/4ϕ(DjRθjlx−k), (2.7)where• ϕ(x) is a smooth bell-shaped function in the horizontal direction andoscillatory in the vertical direction;202.2. Approximation of the normal operator• Dj = diag(2j,2floorleftj/2floorright) is the parabolic scaling matrix;• Rθjl is the rotation matrix over angles θjl = 2pil2−floorleftj/2floorright with 0≤θjl <2pi;• k = (k1,k2)∈Z2 the translation parameters.In the frequency domain, curvelets are compactly supported and eachelement hatwideϕµ(ξ) is localized near the symmetric wedgeWjl ={±ξ,2j ≤|ξ|≤2j+1,|θ−θjl|≤pi.2−floorleftj/2floorright},The number of wedges doubles every other scale doubling and hence thedirectional selectivity increases for finer scales (see Fig. 5.1). For each µ∈M, the curvelet coefficients are given by the inner productcµ =〈f,ϕµ〉:=integraldisplayR2f(x)ϕµ(x)dx (2.8)of a function f ∈L2(R2) with a curvelet ϕµ ∈L2(R2). Curvelets are tightframes, so we have the following reconstruction formulaf =summationdisplayµ∈Mcµϕµ = CTCf, (2.9)by virtue of the energy conservationsummationdisplayµ∈M|cµ|2 =bardblfbardblL2(R2), ∀f ∈L2(R2), (2.10)212.2. Approximation of the normal operatorwith C and CT the curvelet decomposition and composition, respectively.Refer to Candes et al. (2006a), for details on the construction of the fastdiscrete curvelet transform (FDCT) in dimension two.k10-0.25-0.50 0.25 0.50-0.50k2-0.2500.250.5001002003004005000 100 200 300 400 500SamplesSamplesFigure 2.3: Spatial and frequency representation of curvelets, (a) four dif-ferent curvelets in the spatial domain at three different scales, (b) dyadicpartitioning in the frequency domain, where each wedge corresponds to thefrequency support of a curvelet in the spatial domain. This figure illustratesthe microlocal correspondence between curvelets in the physical and Fourierdomain. Curvelets are characterized by rapid decay in the physical spaceand of compact support in the Fourier space. Notice the correspondencebetween the orientation of curvelets in the two domains. The 90◦ rotationis a property of the Fourier transform.2.2.3 Diagonal approximation of PsDO’sData, and the scattering operator K can be so defined that the normaloperator is a PsDOof order 0, which is real and self-adjoint, and has ahomogeneous principal symbol a(x,ξ). We will show that we can approxi-mate Ψf in the curvelet domain by application of a diagonal matrix, if the222.2. Approximation of the normal operatorwavenumbers in f are sufficiently high, relative to the smoothness of thesymbol of Ψ.So let Ψ = Ψ(x,D) be a pseudodifferential operator of order 0, withhomogeneous principal symbol a(x,ξ). Assuming the operator is self-adjointand real, the principal symbol a is also real and has an even symmetry underthe transformation ξ↔−ξ. In addition, we make the technical assumptionthat Ψ is compact, i.e., if f has compact support, then Ψf also has compactsupport.Curvelets in R2 are denoted by ϕµ and we define |µ| = j. The scalej ranges from some positive number jmin, to infinity (or jmax in an imple-mentation). The central position and wave vector for the curvelet will bedenoted by (xµ,ξµ), we let θµ = ξµ/bardblξµbardblbe the angle of the curvelet.Next, we give some consequences of the localization of the curvelet inthe space and Fourier domains. The support of a curvelet in the Fourierdomain is contained in one of the domains ±ξ2 > epsilon1|ξ1|, or ±ξ1 > epsilon1|ξ2|,for some small ε, say we are in the first case, then ξ2 can be used as theradial coordinate, with ξ1/ξ2 the angular coordinate in the Fourier domain.From the localization properties w.r.t. the angular coordinate, it follows that(ξ1/ξ2−ξµ,1/ξµ,2) is bounded by C12−floorleftj/2floorright on the support of the curveletand with C1 some finite positive constant (from here on we assume that anyconstant Ci is finite and positive), so thatbardbl(ξ1/ξ2−ξµ,1/ξµ,2)hatwideϕµbardblL2(R2)≤C12−floorleftj/2floorright. (2.11)Localization in the space domain follows from the smoothness of the defining232.2. Approximation of the normal operatorfunctions in the Fourier domain, in the radial direction we havebardblθµ·(x−xµ)ϕµbardblL2(R2)≤C22−jwith C2 another finite constant. For the other directions we have a weakerestimate, as the curvelet is in the Fourier domain narrower, and therefore inthe space domain wider in the directions normal to θµbardbl(x−xµ)ϕµbardblL2(R2)≤C32−floorleftj/2floorright.In fact parallel results hold for curvelets in Rd.We now compare the application of a PsDOto a curvelet with simplymultiplying by the value of the symbol at (xµ,ξµ). The result will be provedfor x in Rd.Lemma 1 Suppose a is in the symbol class S01,0 (Hormander, 1987), then,with Cprime some constant, the following holdsbardbl(Ψ(x,D)−a(xν,ξν))ϕνbardblL2(Rn)≤Cprime2−|ν|/2. (2.12)To approximate Ψ, we define the sequence u := (uµ)µ∈M = a(xµ,ξµ) withMthe set of all the curvelet indices. Let DΨ be the diagonal matrix withentries given by u. Next we state our result on the approximation of Ψ byCTDΨC.Theorem 1 The following estimate for the error holdsbardbl(Ψ(x,D)−CTDΨC)ϕµbardblL2(Rn)≤Cprimeprime2−|µ|/2, (2.13)242.2. Approximation of the normal operatorwhere Cprimeprime is a constant depending on Ψ.This main result is proved in Appendix A and it shows that the approxi-mation error for the diagonal approximation decreases for increasingly finerscales, j. The approximation derives from the property that the symbol isslowly varying over the support of a curvelet when the velocity model issufficiently smooth.2.2.4 Decomposition of the normal operatorBy virtue of Theorem 1, the normal operator can approximately be factor-ized intoparenleftbigΨϕµparenrightbig(x) similarequal parenleftbigCTDΨCϕµparenrightbig(x) (2.14)= parenleftbigAATϕµparenrightbig(x)with A := CT√DΨ and AT :=√DΨC. Because the seismic reflectivity canbe written as a superposition of curvelets, the ϕµ can be replaced with themodel m. The normal equation (cf. Eq. 2.4) can now be rewritten in thefollowing approximate formy(x) = parenleftbigΨmparenrightbig(x)+e(x)similarequal parenleftbigAATmparenrightbig(x)+e(x)= Ax0 +e. (2.15)The symbolsimilarequalin these expressions is used to indicate that the approxima-tion is only valid for models that are sufficiently close to a reference vector.252.2. Approximation of the normal operatorMoreover, the identity only holds to within some constant. This approxi-mation for the forward model forms the point of departure for our seismicamplitude recovery algorithm.To illustrate the above results, a numerical approximation of a PsDOisapplied to a number of curvelets at different location, angles and scales. ThisPsDOcorresponds to the normal operator associated with a seismic experi-ment with a velocity model containing a low-velocity lens (see Fig. 2.4(a)).TheresultsoftheseexperimentsaresummarizedinFig.2.2. Fig.2.2(b),2.2(d)display the results of applying the PsDOto three different coarse- and fine-scale curvelets (Fig. 2.2(a),2.2(c)). The Fourier-domain images are plottedin Fig. 2.2(f),2.2(h). Comparing the imaged curvelets with their originals,shows that curvelets remain invariant as long as they lie in the support ofthe normal operator. Steep dipping curvelets are outside the support of theoperator and this explains the deterioration of their amplitudes for steepdips. As predicted by Theorem 1, the normal operator becomes more diag-onal for increasingly finer scales, an observation reflected in the behavior ofthe coarse- versus the fine-scale curvelets. Not only is this observation con-sistent with the microlocal correspondence of curvelets reported by Candesand Donoho (2000b) but it is also consistent with the fact that the PsDO’scorrespond to high-frequency approximations of the normal operator.2.2.5 Estimation of the diagonalThe discrete normal operator: Before proceeding with the formulationof the seismic image recovery problem, a method is introduced to estimatethe diagonal entries of DΨ. These entries are calculated from applying a262.2. Approximation of the normal operatordiscrete implementation of the normal operator to some appropriately cho-sen reference vector. This discrete normal operator is given by the com-position of the discretized scattering operator that relates the discretizedreflectivity to the discretized data and its adjoint the migration operator.This compound operator corresponds to a matrix-free implementation forΨ := KTK∈RM×M with K∈RL×M the scattering matrix and L and Mthe length of the data and image vectors.The discrete curvelet transform: For the numerical implementationof the curvelet transform, the fast discrete curvelet transform via wrap-ping (Candes et al., 2006a) is used, which corresponds to a matrix-free im-plementation of the tall curvelet decomposition matrix C ∈ RN×M withN = #{M}greatermuch M with the symbol # denoting ’number of’. This trans-form yields a redundant coefficient vector c = Cr with c ∈ RN. For thischoice of curvelet transform, the pseudo-inverse equals the transpose, i.e.,CTc = C†c. The transform is a numerical isometry that preserves energy,i.e.,bardblrbardbl=bardblCrbardbl, so we have CTCr = Ir in the lscript2 sense. Since the discretecurvelet representation is overcomplete, with a moderate redundancy (a fac-tor of roughly 8 for d = 2), the converse is not the identity, i.e., CCTr isa projection. Because CCT is a projection, not every curvelet vector is theforward transform of some discrete vector f.Non-uniqueness: The estimation of the diagonal DΨ involves the cal-culation of a diagonal curvelet-domain weighting vector that approximatesthe action of the normal operator on an appropriately chosen reference vec-272.2. Approximation of the normal operatortor r, i.e., CTDΨCrsimilarequalΨr (cf. Eq. 2.13). This estimation problem can beformulated in terms of a least-squares minimization problem. Define the’data’ vector as b = Ψr, i.e., the migrated demigrated reference vector. LetDΨ = diag(u), with u ∈ RN, be the diagonal curvelet-domain weightingmatrix with coefficients (uµ)µ∈M on its diagonal and calculate the curvelettransform of the reference vector v = Cr. The diagonal estimation problemcan now be formulated as finding the vector u that minimizes bardblb−Pubardbl2with P := CT diag(v) and b∈RM the known demigrated-migrated refer-ence vector. Since the curvelet transform is redundant (M lessmuchN), there aremany possible solutions to this minimization problem, which gives rise tonon-uniqueness.Regularization via smoothness in phase space: One of the conse-quences of the non-uniqueness is the existence of negative entries in theestimated diagonal. These entries are inconsistent with the fact that Ψ ispositive definite. Therefore, a regularization is needed that leads to an esti-mate for the diagonal with positive entries only. We argue that this can beaccomplished by imposing smoothness in phase space. This smoothness isa property of the symbol that describes the behavior of the normal opera-tor. In the curvelet domain, this smoothness translates to a smoothness inthe amplitudes amongst neighboring points on the curvelet grid. For eachentry in the diagonal, this smoothness can be imposed by including an ad-ditional penalty term in the least-squares formulation for the estimation of282.2. Approximation of the normal operatorthe diagonal. This damped least-squares problem can be formulated astildewideu = argminu 12bardblb−Pubardbl22 +η2bardblLubardbl22, (2.16)where L = bracketleftbigDT1 DT2 DTθ bracketrightbigT is a so-called sharpening operator, penalizingfluctuations amongst neighboring coefficients in u. D1,2 contain first-orderdifferences at each scale in the x1,2 directions, i.e., D1,2 =bracketleftBigDjmin1,2 ···Djmax1,2bracketrightBigwith Dj1,2 the differences at the jth scale in the spatial directions, and Dθ thefirst-order differences in the θ direction, i.e., Dθ =bracketleftBigDjminθ ···DjmaxθbracketrightBigwithDjθ the difference in the angle direction at the jth scale. These differencesare scale dependent because the curvelet grid changes for each scale. Theemphasis on the penalty term with respect to the misfit is controlled by η.Estimation with positivity constraints: To ensure a correct balancebetween the misfit and the penalty functional, Eq. 2.16 is solved for a seriesof increasing Lagrange multipliers η. For each η, the following system ofequations is inverted PηLu≈b0 (2.17)with ≈ the approximation in least squares sense and 0 a column vectorwith N zero entries. For a given reference vector, the diagonal estimationprocedure consists of the following steps. Apply the normal operator tothe reference vector and calculate its curvelet coefficients. Set the Lagrangemultiplier to ηmin and solve Eq. 2.16 by inverting the system of equations inEq. 2.17. Subsequently, increase η by ∆η and repeat this procedure until all292.2. Approximation of the normal operatordiagonal elements of tildewideu are nonnegative. Even though we do not have a proofof convergence for this method, increasing the η leads to positive entries forlarge enough η. The steps of this estimation procedure are summarized inTable. 2.1.Calculate: b = Ψr and v = Cr.Set: η = ηmin;while∃(˜uµ)µ∈M < 0 doSolvetildewideu = argminu 12bardblb−Pubardbl22 +η2bardblLubardbl22Increase the Lagrange multiplierλ = η +∆ηend whileTable 2.1: Algorithm for the estimation of the diagonal via regularized least-squares inversion. The Lagrange multiplier, η, is increased up to the pointthat all entries in the vector for the diagonal are positive.Example: The procedure outlined above is tested on a stylized exam-ple with a velocity model given by a low-velocity lens (Fig. 2.4(a)) and abandwidth-limited reflectivity (Fig. 2.4(b)), consisting of three events. Thelow-velocity zone is representative of a gas lens in the overburden. Data isgenerated from this reflectivity with a zero-order linearized Born modelingand consists of 32 shot records with 500 receivers each. For more detailson the migration code, refer to Symes (2006b) and to the numerical ex-periment section at the end of this paper. The data is subsequently imaged(Fig. 2.4(c)). This image and the bandwidth-limited reflectivity (Fig. 2.4(b))define the ’data’ vector b (cf. Eq. 2.17) and the reference vector r used toestimate tildewideu according to the algorithm presented in Table 2.1. Fig. 2.5 dis-302.3. Stable seismic image recoveryplays the curvelet vectors for increasing η ={0.01, 0.1, 1, 10}. These plotsclearly illustrate that tildewideu becomes positive for η large enough (η≥1 in thiscase). Fig. 2.4(d) contains the diagonal approximation of the normal opera-tor for the diagonal estimated with η = 1. Comparison between the resultsof applying the normal operator (Fig. 2.4(c)) and its diagonal approximation(Fig. 2.4(d)), shows a nice correspondence. The relative lscript2-error betweenthe actual image its approximation is 6.1%. This value corresponds to amoderate size for the constant Cprimeprime in the bound for the approximation error(Theorem 1).2.3 Stable seismic image recoveryOur formulation for the seismic amplitude recovery problem banks on twofundamental properties of curvelets, namely their ability to detect wave-fronts and their invariance under the normal operator. The combinationof these two properties allow for a robust solution of the recovery problemin terms of a nonlinear optimization problem. In this section, the differentsteps that lead to this formulation are discussed starting with an argumen-tation why curvelets compress seismic images, followed by a quick review oflscript1-promoted image recovery, its extension to the seismic situation and itssolution by a cooling method.2.3.1 Curvelet frames for seismic imagesThe sedimentary crust consists of sheet-like layers that correspond to func-tions with singularities along piece-wise smooth curves. These singularities312.3. Stable seismic image recovery(a) (b)(c) (d)Figure 2.4: Stylized example of the diagonal estimation for a smooth velocitymodel and a reflectivity consisting of three reflection events including a fault,(a) the smooth background model with the low-velocity lens. This velocitymodel is used to calculate the linearized scattering and migration operators,(b) the bandwidth limited reflectivity, (c) the imaged reflectivity from dataconsisting of 32 shot records with 500 traces each, (d) the approximatednormal operator for diag(tildewideu) estimated with η = 1. The bandwidth limitedreflectivity in (b) and the image in (c) served as input for the referenceand ’data’ vectors for the diagonal estimation procedure outlined in Table2.1. Remarks: the main dimming of the normal operator is captured by thediagonal approximation. The relative lscript2-error between the actual and theapproximate normal operator is 6.1%.322.3. Stable seismic image recovery(a) (b)(c) (d)Figure 2.5: Estimates for the diagonal tildewideu are plotted in (a-d) for increasingη = {0.01, 0.1, 1, 10}. The diagonal is estimated according the procedureoutlined in Table 2.1 with the reference and ’data’ vectors, v and b, plottedin Fig. 2.4(b) and 2.4(c). As expected the diagonal becomes more positivefor increasing η.332.3. Stable seismic image recoveryare associated with rapid variations in the medium fluctuations at whichseismic waves reflect. During seismic imaging, these singularities are recov-ered. Examples of singularities imaged from real seismic data can be foundin Fig. 2.1(a). This image can be considered as a bandwidth limited versionof a function with singularities on piece-wise C2 curves that may includefaults and pinch outs (see e.g., Fig. 2.4(b) and 2.8(a) that are examplesof synthetic reflectivity models). Fourier (and also wavelet) transforms areknown to perform poorly on this type of functions, (Candes and Donoho,2000a; Candes and Guo, 2002). This poor performance can be explained bythe inability of the Fourier transform to localize and the lack of directivity ofwavelets. Only when oriented perpendicular to the interfaces, wavelets decayrapidly. Since wavelets lack directionality, they can not resolve these direc-tions, known as the wavefront set. This inability explains why wavelets donot significantly improve the compression rate for seismic images. Curveletframes, on the other hand, are designed to detect the wavefront set (Candesand Donoho, 2000a, 2004, 2005a,b) and lead, as shown in Fig. 2.6, to betterempirical compression rates for real and synthetic seismic images. This im-proved compression rate explains the accurate reconstruction in Fig. 2.1(a)of a real seismic image from only 3% of the curvelet coefficients and is con-sistent with the theoretical nonlinear approximations rates reported in theliterature (see e.g. Candes and Donoho, 2000a). These rates correspond toa rate of O(P−2) for curvelets, with P the number of largest entries usedin the approximation, while Fourier and wavelets only attainO(P−1/2) andO(P−1), respectively. The empirical rates plotted in Fig. 2.6 seem to confirmthese theoretical findings.342.3. Stable seismic image recovery(a) (b)Figure 2.6: Decay rate for the curvelet coefficients compared to the decayfor the coefficients of Fourier (dashed line) and discrete wavelet transforms(dot-dashed line), (a) decay rate for the sorted curvelet coefficients (solidline) of the real migrated image included in Fig. 2.1(a), (b) the same as(a) but now for the synthetic reflectivity of the SEG/EAGE AAprime salt model(cf. Fig. 2.8(a)). The decay for the curvelet transform clearly comparesfavorably to these other two transforms.2.3.2 Stable recoveryBecause of the redundancy of the curvelet transform, CCT is a projectionand this means that not every curvelet vector is a forward transform of a vec-tor f. Therefore, the vector x0 can’t readily be calculated from f = CTx0,because there exist infinitely many coefficient vectors whose inverse trans-form equals f. Recent work on in the field of compressive sensing has shownthat rectangular matrices such as the curvelet matrix, C ∈ RN×M withN greatermuch M, can stably be inverted through a nonlinear sparsity promotingoptimization program. These nonlinear recovery methods require fast decayfor the magnitude sorted curvelet coefficients (for details refer to the ex-tensive literature on stable recovery also known as compressed sensing. Seee.g. Starck et al., 2004; Elad et al., 2005; Candes et al., 2006b).Following these results, the vector x0 can be approximately recovered352.3. Stable seismic image recoveryfrom noise-corrupted datay = Ax0 +n (2.18)with A := CT. As long as the signal is sufficiently compressible x0 cansuccessfully be recovered to within the noise level. This recovery involvesthe solution of a constrained convex optimization problemP1 :tildewidex = argminxbardblxbardbl1 := summationtextNj=1|xj| subject to bardbly−Axbardbl2≤epsilon1tildewidem = Ax.(2.19)For certain matricesA, the solution of this unconstrained optimization prob-lem lies to within the noise level (see e.g. Candes et al., 2006b; Elad et al.,2005). The optimization problem P1 is known as the constrained variationof the LASSO (Tibshirani, 1997) and the basis-pursuit denoising (BPDN)(Chen et al., 2001) algorithms.As part of the optimization, the sparsity vector is fitted within the tol-erance epsilon1. This tolerance depends on the noise level given by the standarddeviation of the noise vector. Since n1,n2,···nM ∈N(0,σ2), the probabil-ity of bardblnbardbl22 exceeding its mean by plus or minus two standard deviationsis small. The bardblnbardbl22 is distributed according the χ2-distribution with meanM·σ2 and standard deviation√2M·σ2. By choosing epsilon12 = σ2(M +ν√2M)with ν = 2, we remain within the mean plus or minus two standard devi-ations. Following (Elad et al., 2005), the above constrained optimizationproblem (P1), is replaced by a series of simpler unconstrained optimization362.3. Stable seismic image recoveryproblemsPλ :argminxbardbly−Axbardbl22 +λbardblxbardbl1tildewidem = Atildewidex.(2.20)These optimization problems depend on the Lagrange multiplier λ. A cool-ing method is used, where Pλ is solved for a Lagrange multiplier λ that isslowly decreased from a large starting value λ1 <bardblATybardbl∞. The optimal tildewidexis found for the largest λ for whichbardbly−Atildewidexbardbl2≤epsilon1. During the optimization,the underdetermined frame matrix A is inverted by imposing the sparsitypromoting lscript1-norm. This norm regularizes the inverse problem of findingthe unknown coefficient vector (see also Daubechies et al., 2005). We referto (Donoho et al., 2006; Tropp, 2006) for the recovery conditions for Eq.’s2.19 and 2.20.2.3.3 Solution by iterative thresholdingFollowing(Daubechiesetal.,2005;Eladetal.,2005)andideas datingbackto(Figueiredo and Nowak, 2003), Eq. 2.27 is solved by an iterative thresholdingtechnique that derives from the Landweber descent method. The methodconsist of an outer loop during which the Lagrange multiplier is lowered andan inner loop during which Pλ is solved. The inner loop is initialized by thesolution obtained from the previous outer loop starting with zero vector forthe first iteration. Aftermiterations of the outer cooling loop, the estimatedcoefficient vector is computed for fixed λ by the following inner loopxm+1 = Tλparenleftbigxm +AT(y−Axm)parenrightbig (2.21)372.3. Stable seismic image recoverywith λ = λm andTλ(x) := sgn(x)·max(0,|x|−|λ|). (2.22)As shown by (Daubechies et al., 2005), this iteration for fixed λ convergesto the solution of Eq. 2.27 for m large enough and bardblAbardbl< 1. The cost foreach iteration is a curvelet synthesis and subsequent analysis.2.3.4 Stable seismic recoveryThe main purpose is to estimate the relative amplitudes of seismic imagesfrom data that are possibly contaminated with noise. This data is repre-sented by the discretized linearized forward modeld = Km+n (2.23)withKthe discrete linearized Born modeling operator. Applying the adjointof this operator to the data vector creates a noisy image (cf. Eq. 2.4)y = Ψm+e. (2.24)With the curvelet decomposition of the normal operator (cf. Eq. 3.9) indiscrete form, i.e.,ΨmsimilarequalAATm, (2.25)382.3. Stable seismic image recoverythis expression can be rewritten into the following approximate formysimilarequalAx0 +e (2.26)with A := CTΓ and Γ := √DΨ. To avoid instabilities related to smallentries in the estimated diagonal u, we set DΨ = diag(tildewideu + δ)/δ with δsome small parameter. Comparing this image representation with the onein Eq. 2.18 shows that these expressions are equivalent, aside from the noiseterm. While the noise in Eq. 2.18 is white and Gaussian, the noise termin Eq. 2.26 is colored by the migration operator. The expectation for thecovariance of this colored noise term equals E{eeT} = σ2Ψ similarequal σ2AAT.Because the nonlinear recovery of the vector x0 involves the inversion ofthe matrix A in which the covariance is factored, the noise term is approxi-mately whitened. To understand this whitening, consider the pseudo-inverseA† = Γ−1C which, when applied to e, approximately whitens the noise inthe curvelet domain. This whitening is a well-known property of WVDtechniques.The nonlinear recovery of seismic images now involves solvingPprime1 :tildewidex = argminxbardblxbardbl1 := summationtextNj=1|xj| subject to bardbly−Axbardbl2≤epsilon1tildewidem = parenleftbigATparenrightbig† tildewidex.(2.27)The main assumption underlying this formulation is that the curvelet vectorof the model remains (approximately) sparse after application of the normaloperator. This assumption is valid when the curvelets remain sufficiently392.4. Sparsity- and continuity-enhanced seismic image recoveryinvariant under this operator. In that case, the nonlinear program (Pprime1)approximately inverts the normal operator in two stages, namely duringthe lscript1-norm regularized inversion of the synthesis matrix, yielding a sparseestimate for the curvelet coefficient vector, and during the calculation ofthe model vector via the pseudo-inverse of AT. During each stage, the’square-root’ of the normal operator is approximately inverted.2.4 Sparsity- and continuity-enhanced seismicimage recoveryThe above formulation provides a stable recovery for the amplitudes of seis-mic images, the relative strengths of the fluctuations in m, without thenecessity of multiple evaluations of the normal operator. This recovery isonly accurate when the proposed diagonal approximation (cf. Eq. 3.3) pro-vides an adequate approximation to the normal operator. The accuracyof this approximation depends on the available scales in the seismic image(cf. Theorem. 1), on the complexity of the background velocity model, theacquisition geometry and the closeness of the reference r to the actual un-known model m. For the remainder of this paper, we assume the acquisitionto be ideal and the reference vector to be close to the actual model.Anisotropic diffusion: Aside from the aforementioned factors that enterinto the seismic recovery problem, there is the issue of how to limit spuriouscurvelet artifacts. These artifacts are either related to the so-called “pseudoGibb’s” phenomenon (or better side-band effects (see Candes and Guo,402.4. Sparsity- and continuity-enhanced seismic image recovery2002)), inherent to curvelet or other harmonic expansions. To reduce thesespurious artifacts, the sparsity enhancing penalty functional in Eq. 2.27 iscomplemented with an additional continuity-enhancing penalty functional.This functional enhances the continuity along the imaged reflectors, whichtend to be smooth in the tangential directions. By applying an anisotropicsmoothing technique, the ’wavefront’ set of the imaged reflectivity is pre-served. This technique differs from commonly used edge-preserving penaltyfunctionals such as total variation (TV) (see e.g. Claerbout and Muir, 1973;Schertzer, 2003) that tend to remove the oscillatory behavior of the reflectorsin the normal direction.The anisotropic-diffusion penalty term (see e.g. Fehmers and Hocker,2003) is given byJc(m) =bardblΛ1/2∇mbardbl22, (2.28)with ∇ the discretized gradient matrix defined as ∇ = bracketleftbigDT1 DT2bracketrightbigT. Theblock-diagonal matrix Λ is location dependent (see Fig. 2.10, which plotsthe gradients) and rotates the gradient towards the tangents of the reflectingsurfaces. This rotation matrix is given byΛ[r] = 1bardbl∇rbardbl22 +2υ+D2r−D1rparenleftbigg+D2r −D1rparenrightbigg+υI, (2.29)with Di the discretized derivative in the ith coordinate direction and υ aparameter that controls the fluctuations for regions where the gradient issmall. Following Black et al. (1998), this control parameter is set propor-tional to the median of|∇r|with||the length of each gradient vector (white412.4. Sparsity- and continuity-enhanced seismic image recoveryarrows in Fig. 2.10). Similar to the diagonal approximation, a reference vec-tor derived from the migrated image (cf. Eq. 2.24) is used to calculate thetangential directions of the reflecting surfaces.By combining the two different penalty terms that promote sparsity andcontinuity, we finally arrive at our formulation for the seismic-amplituderecovery problemPepsilon1 :argminxJ(x) subject to bardbly−Axbardbl2≤epsilon1tildewidem = parenleftbigATparenrightbig†x,(2.30)in which the composite penalty term J(x) is given byJ(x) = αJs(x)+βJc(x), (2.31)with α, β≥0 and α+β = 1. The Js(x) =bardblxbardbl1 is the lscript1-norm. The secondterm in the penalty term is given by Jc(x) =bardblΛ1/2∇parenleftbigATparenrightbig†xbardbl22. Becausethe optimization is carried out over x and not over the model vector m, thisexpression includes a pseudo-inverse that is calculated with a few iterationsof the LSQR algorithm (Paige and Saunders, 1982).The cooling method: The above nonlinear optimization problem (Pepsilon1)is solved with a cooling method (see e.g. Starck et al., 2004). This methodconsists of a series of thresholded Landweber iterations that solve a series ofunconstrained subproblems for decreasing λ. Since this method only requiresknowledge of the Jacobians at each iteration, it is relatively straightforwardto include the Jacobian of the additional continuity-enhancing penalty func-422.4. Sparsity- and continuity-enhanced seismic image recoverytional Jc(x). For data given by Eq. 2.18, the iterations of the cooling methodfor a particular cooling parameter λ consist of the following three main steps:Step 1: update of the Jacobian of 12bardbly−Axbardbl22:x←x+AT (y−Ax); (2.32)Step 2: project onto the lscript1 ball S ={bardblxbardbl1≤bardblx0bardbl1}by soft thresholdingx←Tλ(x); (2.33)Step 3: project onto the anisotropic diffusion ball C ={x : J(x)≤J(x0)}byx←x−β∇xJc(x) (2.34)with∇xJc(x) = 2A†∇·parenleftBigΛ∇parenleftBigparenleftbigATparenrightbig†xparenrightBigparenrightBig. (2.35)Remember that steps 1&2 withbardblAbardbl≤1 converge to the solution of Eq. 2.27for a fixed λ. During step 3, the coefficients are updated according to thegradient of the anisotropic diffusion norm designed to reduce spurious arti-facts. The different steps of our algorithm are summarized in Table 2.2. Toensure sparsity, the algorithm is started by setting x0 = 0.432.5. Practical considerationsInitialize:m = 0;x0 = 0;y = KTd;Choose:M and LbardblATybardbl∞ > λ1 > λ2 >···whilebardbly−Atildewidexbardbl2 > epsilon1 dom = m+1;xm = xm−1;for l = 1 to L doxm = Tλmparenleftbigxm +AT (y−xm)parenrightbig{Iterative thresholding}end forAnisotropic descent update;xm = xm−β∇xm Jc(xm);end whiletildewidex = xm; tildewidem = parenleftbigATparenrightbig† tildewidex.Table2.2: Sparsity-andcontinuity-enhancingrecoveryofseismicamplitudes.2.5 Practical considerationsThe presented method depends on the availability of a reference vector withan image that is sufficiently close to the unknown model. As in most Krylov-subspace solvers, we derive the reference vector from the migrated image(matched filter) as an initial guess. Since waves are subject to sphericalspreading, a depth-dependent correction is applied to the migrated image(cf. Eq. 2.24). This corrected image serves as input for the diagonal estima-tion procedure and the calculation of the anisotropic difussion norm. Theimage recovery procedure consists of the following steps1. Obtain the reference vector, r, by imaging the data with Eq. 2.24,442.6. Numerical experimentsfollowed by applying the depth correction.2. Use the reference vector, r, to estimate the diagonal matrix, diag(tildewideu),with the procedure outlined in Table 2.1;3. Define the synthesis and analysis matrices A and AT;4. Calculate the rotation matrix, Λ, (cf. Eq. 2.29) for the continuitypromoting penalty functional;5. Solve the optimization problem Pepsilon1 (Eq. 2.30) with the algorithm out-lined in Table 2.2. This yields the estimate for the amplitude-correctedimage, tildewidem;6. Optionally goto 1 with r := tildewidem;2.6 Numerical experimentsThe recovery method is applied to a synthetic imaging example for data ofthe SEG/EAGE AAprime salt model (O’Brien and Gray, 1997; Aminzadeh et al.,1997). The original velocity model and its smoothed version that defines thebackground velocity model, used by the migration operators, are includedin Fig. 2.7. To avoid multiple reflections, the background velocity modelis averaged over a window of size 720×720m. The reflectivity (velocityperturbation) is defined as the band-pass filtered difference between thesmoothed velocity model and a slightly less smoothed velocity model. Thedensity is assumed to be constant. Fig. 2.8(a) shows the reflectivity modelthat is 3.6km by 15.3km with a grid spacing of 24m in the vertical and452.6. Numerical experimentshorizontal coordinate directions. This reflectivity was used to generate datawith the linearized Born modeling operator.Themigratedimage: Atwo-wayreverse-timemigrationcodewithcheck-pointing was used to conduct the experiments. This two-dimensional (d = 2)time-domain code, developed by Symes (2006b), has absorbing boundaryconditions and is based on the adjoint-state method. The migration oper-ator is derived from the gradient of a descent algorithm for least-squaresinversion (see Symes, 2006b, and the references therein). The data is gen-erated with the adjoint of the migration operator, which corresponds tothe linearized Born approximation. This Born modeling operator generatesdata without nonlinear events such as multiple reflections. The migrationoperator and its adjoint are made zero-order by applying a half-order frac-tional integration on the source wavelet. We choose a source wavelet witha relatively broad temporal frequency band [5−60]Hz. The data consistsof 324 shots with 176 receivers each with a shot at every 48m, yielding amaximum offset of 4.224km. Each geophone records 626 time samples witha sample interval of 8ms, which amounts to 5s of data.The 8000 time-step linearized Born modeling takes about 64min with68 CPU’s on a MPI cluster, while the migration takes about 294min. Thesecomputation times explain why calculating the pseudo-inverse of the normaloperator becomes computationally prohibitive in real applications where theimages are three dimensional (d = 3).Fig. 2.8 shows the original reflectivity and the migrated image obtainedwith Eq. 2.24. The image suffers from amplitude deterioration due to spher-462.6. Numerical experiments(a)(b)Figure 2.7: The SEG/EAGE AAprime salt model (Aminzadeh et al., 1997), (a)The original velocity model, (b) the smoothed velocity model with a windowsize of 720×720m. Remark: Imaging this model is a challenge because ofthe presence of the high-velocity (bright) salt.472.6. Numerical experimentsical spreading and poor illumination. Despite the overall amplitude deteri-oration for steep reflectors and events at increasing depths, the migrationresolves most of the singularities in the image (cf. Fig 2.8(a) and 2.8(b)).(a)(b)Figure 2.8: Imaging of the SEG/EAGE AAprime salt model (Aminzadeh et al.,1997), (a) reflectivity defined in terms of the band-pass filtered differencebetween the smoothed velocity model and a slightly less smoothed velocitymodel, (b) the imaged reflectivity according to Eq. 2.24. The main reflec-tion events are present but suffer from deteriorated amplitudes, especiallyunder the high-velocity salt and for steep reflectors and faults.482.6. Numerical experimentsDiagonalestimation: Afterapplyingcorrectionsforthesphericalspread-ing, the migrated image is considered as reference vector. This referencevector is used to compute the diagonal approximation for the normal op-erator and for the calculation of the gradients. The image deteriorationbetween the reference vector (Fig. 2.9(a)) and the demigrated-migrated ref-erence vector (Fig. 2.9(b)) is similar to the amplitude deterioration betweenthe ’unknown’ true reflectivity and migrated image (cf. Fig. 2.8), there-for the images in Fig’s. 2.9(a) and 2.9(b) should be adequate as input forthe diagonal estimation procedure outlined in Table 2.1. To avoid insta-bilities related to small entries in the estimated diagonal diagtildewideu, we setΓ = diag(tildewideu + δ)/δ with δ = 0.2. The result of applying the estimateddiagonal to the bandwidth-limited reflectivity of Fig. 2.8(a) is included inFig. 2.9(c) and shows that most of the imprint of the normal operator iscaptured by the diagonal approximation.Seismic amplitude recovery for the SEG AAprime model: With the di-agonal approximation in place, the amplitude corrections are applied withthe procedure presented in Table 2.2. The anisotropic diffusion penaltyterm is calculated from the gradients plotted in Fig. 2.10. The results of theamplitude recovery are summarized in Fig. 2.11. Fig. 2.11(a) depicts theresult obtained by only applying a lscript1-norm penalty (β = 0). The image iscalculated by lowering the Lagrange multiplier λ in 20 steps from a valuethat corresponds to a threshold that removes all but 5% of the coefficientsto a Lagrange multiplier that only removes 1% of the curvelet coefficients.For each intermediate λ, the inner loop is repeated 5 times (L=5). The492.6. Numerical experiments(a)(b)(c)Figure 2.9: Images that are input to the diagonal approximation schemeoutlined in Table 2.1, (a) the reference vector, r, derived from Fig. 2.8(b),after applying corrections for the spherical-divergence and receiver-array ta-per, (b) the demigrated-migrated reference vector, b = Ψr, that servesas ’data’ for the diagonal estimation, (c) diagonal approximation on theoriginal bandwidth-limited reflectivity plotted in (a). This diagonal approx-imation by AATm captures the normal operator, Ψm, quite well.502.6. Numerical experimentsGradient of the reference vectorlateral (m)depth (m)2000400060008000100001200014000500100015002000250030003500Figure 2.10: Gradient vectors for the reference vector r plotted in Fig.2.9(a). The gradients (white ↑’s) are used to calculate the tangential di-rections along which the additional anisotropic smoothing is applied.512.6. Numerical experimentsresult of this norm-one optimization are plotted in Fig. 2.11(a). The im-age shows a clear improvement over the image plotted in Fig. 2.9(a). Notonly is the spherical divergence corrected in a stable manner, but the ampli-tude of the steep reflections events are recovered as well. There are, however,some small remaining artifacts related to side-band effects (Candes and Guo,2002). By including the continuity-enhancing penalty functional most of theaforementioned artifacts can be removed as can be seen from Fig. 2.11(b).In this example, the sparsity- and continuity-enhancing norms are weightedequally, i.e. α = β = 1/2 in Eq. 2.31. To further illustrate the migrationamplitude recovery, detailed plots for different traces of the recovered theamplitudes are included in Fig. 2.12. Fig. 2.12(a) compares vertical tracesof the original, migrated and amplitude recovered images. The traces arenormalized with respect to the first reflector and the plots show that theamplitudes are nicely recovered. Similarly, the amplitudes along the hori-zontal reflector at z = 3432m, are nicely recovered. The improvement inthe recovered amplitudes is also apparent from Fig. 2.12(c) that contains theaverage vertical-wavenumber amplitude spectra of the original, the imagedand the recovered reflectivity.Seismic amplitude recovery from noisy data: So far, the exampleswere noise free. In practice, seismic data volumes contain several sources ofclutter, ranging from coherent nonlinear events, such as multiple reflections,to incoherent measurement errors. In this paper, only the recovery fromincoherent contamination is considered. Sources of coherent clutter canbe removed by using curvelet-based signal separation techniques reported522.6. Numerical experiments(a)(b)Figure 2.11: Images after nonlinear recovery (Pepsilon1). (a) result with the lscript1-norm recovery only (β = 0); (b) recovery with the combined lscript1 and con-tinuity recovery for α = β = 1/2. The amplitudes are recovered. Theanisotropic diffusion successfully removes the artifacts.532.6. Numerical experiments(a) (b)(c)Figure 2.12: Recovery of the amplitudes after normalization with respectto the first event, (a) seismic traces of the original reflectivity (dot-dashedline), the migrated reflectivity (dashed line) and the recovered reflectivity(solid line), (b) amplitudes along the bottom most horizontal reflector, (c)average amplitude spectra of the depth-dependent reflectivity. The originaland recovered spectra are normalized to match. Remarks: The nonlinearrecovery corrects for the amplitudes and restores the spectrum.542.6. Numerical experimentselsewhere (Herrmann et al., 2007). An image computed from data with asignal-to-noise ratio (SNR) of 3dB is included in Fig. 4.3(a). As shownin Eq. 2.24, the image now contains an additional contribution due to thenoise. This noise contribution is clearly visible throughout the image as anon-stationary clutter and leads to a deterioration of the image quality. Theclutter in the image space, however, is significantly smaller than in the dataspace, as can be observed from the noisy shot record plotted in Fig. 2.14(b).The difference in noise levels between the data and image spaces can beexplained by the 154-fold redundancy of the data space compared to theimage space (L = 154×M).The result for the amplitude recovery of the noisy data are included inFig. 4.3(b) and this figure shows that the migration amplitudes can stablybe recovered, leading to a significant improvement in the SNR for the image(9.2dB for the image in Fig. 4.3(b) as opposed to (1.5dB for the noisyimage which includes the error due to the imprint of the normal operator).This improved image was obtained using the algorithm of table 2.2, withthe same settings as in the noise-free example except for the lowest valueof the Lagrange multiplier, which is now set to a larger value dependingcalculated from the noise level. From Fig. 4.3(b), one can see that the noisecontamination has largely been removed. Moreover, the amplitudes havebeen restored in particular for the steep events at depth. Data generatedfrom the estimated image, tildewided = Ktildewidem shows a significant removal of thenoise (cf. Fig. 2.14(b)-2.14(c)), with reflection events that match the noise-free data plotted in Fig. 2.14(a). This visual improvements leads to animprovement of SNR for the data of 19.2dB.552.7. Conclusions2.7 Conclusions2.7.1 Initial findingsThe method presented in this paper combines the compression of images bycurvelets with the invariance of this transform under the normal operator.This combination allows for a formulation of a stable recovery method forseismic amplitudes. During the recovery, the normal operator is approxi-mately inverted. Compared to other approaches for migration scaling, thepresented method (i) includes a theoretical bound on the L2 −error forthe diagonal approximation in the curvelet domain; (ii) prescribes a proce-dure for the estimation of the diagonal from numerical implementations ofthe imaging operators; (iii) formulates the amplitude recovery problem asa nonlinear optimization problem, where the inversion of the diagonalizednormal operator is regularized by imposing sparsity in the curvelet domainand continuity along the imaged reflectors.As long as the background velocity model is sufficiently smooth andthere is a reference vector available close enough to the actual reflectivity,the normal operator can be approximated by a diagonal weighting in thecurvelet domain. The theoretical approximation error is scale-dependentand decreases for finer scales, an observation consistent with the microlocalproperties of the normal operator and curvelets in the high-frequency limit.Estimation of the diagonal exploits smoothness of the symbol by penalizingneighboring entries in the diagonal that fluctuate. These fluctuations areinconsistent with the smoothness of the symbol. For a sufficiently large La-grange multiplier, the regularized least-squares estimation for the diagonal562.7. Conclusions(a)(b)Figure 2.13: Image amplitude recovery for noisy data (SNR 3dB), (a) noiseimage according to Eq. 2.24, (b) image after nonlinear recovery from noisydata (Pepsilon1). The clearly visible non-stationary noise in (a) is removed dur-ing the recovery while the amplitudes are also restored. Steeply dippingreflectors denoted by the arrows are also well recovered.572.7. Conclusions(a) (b)(c)Figure 2.14: ’Denoising’ of a shot record, (a) the noise-free data received bythe receiver array, (b) noisy data with a SNR of 3dB, (c) forward modeleddata after amplitude recovery. Observe the significant improvement in thedata quality, reflected in an increase for the SNR of 19.2dB. 582.7. Conclusionsleads to positive entries. Numerical experiments with the diagonal approx-imation showed a moderate-sized constant for the bound of the theoreticalerror.The diagonal approximation is used to formulate the seismic amplituderecovery in terms of a constrained optimization problem. The amplitude-correctedimageisobtainedbysolvingthissparsity-andcontinuity-promotingoptimization problem. The invariance of curvelets under the normal opera-tor preserves the sparsity. The cost of computing the diagonal approxima-tion is one demigration-migration per reference vector, which is much lesscompared to the cost of Krylov-based least-squares inversion. The recov-ery results show an overall improvement of the image quality. The joinedsparsity- and continuity-enhanced image has diminished artifacts, improvedresolution and recovered amplitudes.2.7.2 Extension to prestack imagingThe discussion has been on ’post-stack’ images obtained after applying thezero-offset imaging condition. When the velocity model is correct, imagesare smooth along the redundant coordinate in pre-stack angle gathers. Thissmoothness property has been used in velocity analyses methods such asdifferential semblance (Symes and Carazzone, 1991) and is suitable for theframework laid down in this paper. Besides the lscript1- and continuity norms,the amplitude recovery can also exploit the smoothness along the redundantangle coordinate. The fact that smoothness translates to sparsity is espe-cially exciting because it allows to further exploit the results from the fieldof compressed sensing.592.7. Conclusions2.7.3 Recent related workDuring the revision, we were informed of recent work by William Symeson optimal scaling of reverse-time migration Symes (2006a). This workand recent work by Guitton (2004) attempt to invert the normal operatorby estimating a diagonal approximation from a reference vector, typicallygiven by the migrated image, and the demigrated-migrated reference vec-tor. Smoothness of the diagonal approximation is also imposed by bothauthors by parameterizing the diagonal as a smooth function in the spatialdomain. Our approach goes several steps further by imposing smoothnessin phase space and by making use of the curvelet transform that allows forthe inversion of the diagonal using the methods of stable image recovery.2.7.4 Choice of the reference vectorHaving access to a proper reference vector is a prerequisite for the success ofthe method presented in this paper. Selection of the appropriate referenceis somewhat reminiscent of finding a good initial guess for Krylov subspacemethods. In this paper, we took a data-driven approach by using the mi-grated image (matched-filter result) as a first guess. Our formulation allowsfor multiple reference vectors and we hope to report on appropriate selectionstrategies for reference vectors in the future. This investigation will includean analysis of the sensitivity of our method to the accuracy of the velocitymodel.602.7. ConclusionsAcknowledgmentsThe authors would like to thank the reviewers for their constructive com-ments and suggestions. We also wish to thank the authors of the Fast Dis-creteCurveletTransformformakingtheircodeavailableatwww.curvelet.organd Dr. William Symes for making his reverse-time migration code avail-able. M. O’Brien, S. Gray and J. Dellinger from BP Amoco are thankedfor making the SEG AA’ dataset available. The examples presented inthis paper were prepared with CWP’s Seismic Un*x and with Madagascar(rsf.sourceforge.net). This work was in part financially supported by theNatural Sciences and Engineering Research Council of Canada DiscoveryGrant (22R81254) and the Collaborative Research and Development GrantDNOISE (334810-05) of F.J.H. This research was carried out as part of theSINBAD project with support, secured through ITF (the Industry Technol-ogy Facilitator), from the following organizations: BG Group, BP, Chevron,ExxonMobil and Shell. The first two authors would also like to thank theInstitute of Pure and Applied Mathematics at UCLA, supported by theNSF under grant DMS-9810282, for their hospitality. The research of thethird author was supported by the Netherlands Organization for ScientificResearch, through a VIDI grant no. 639.032.509.61BibliographyAminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D salt and overthrust model.SEG/EAGE 3-D Modeling Series, No. 1: Society of Exploration Geo-physicists, Tulsa.Beylkin, G., 1984, The inverse problem and applications of the generalizedRadon transform: Comm. Pure Appl. Math., 37, 579–599.Black, M., G. Sapiro, and D. Marimont, 1998, Robust anisotropic diffusion:IEEE Trans. Signal Process., 7, 421–432.Bleistein, N., J. Cohen, and J. Stockwell, 2001, Mathematics of multidimen-sional seismic imaging, migration and inversion: Springer.Brandsberg-Dahl, S. and M. V. de Hoop, 2003, Focusing in dip and AVAcompensation on scattering-angle/azimuth common image gathers: Geo-physics, 68, 232–254.Candes, E. and D. Donoho, 2005a, Continuous curvelet transform i: Reso-lution of the wavefront set: Appl. Comput. Harmon. Anal., 19, 162–197.Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006a, Fast discretecurvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899.Candes, E. J. and D. Donoho, 2004, New tight frames of curvelets andoptimal representations of objects with piecewise C2 singularities: Comm.Pure Appl. Math., 57, 219–266.62Chapter 2. Bibliography——–, 2005b, Continuous curvelet transform ii: Discretization and frames:Appl. Comput. Harmon. Anal., 19, 198–222.Candes, E. J. and D. L. Donoho, 2000a, Curvelets – a surprisingly effec-tive nonadaptive representation for objects with edges: Presented at theCurves and Surfaces, Vanderbilt University Press.——–, 2000b, Recovering edges in ill-posed problems: Optimality of curveletframes: Ann. Statist., 30, 784.Candes, E. J. and F. Guo, 2002, New multiscale transforms, minimum totalvariation synthesis: Applications to edge-preserving image reconstruction:Signal Processing, 1519–1543.Candes, E. J., J. Romberg, and T. Tao, 2006b, Stable signal recovery fromincomplete and inaccurate measurements: Comm. Pure Appl. Math., 59,1207–1223.Chavent, G. and R. Plessix, 1999, An optimal true-amplitude least-squaresprestack depth-migration operator: Geophysics, 64, 508–515.Chen, S. S., D. L. Donoho, and M. A. Saunders, 2001, Atomic decompositionby basis pursuit: SIAM Journal on Scientific Computing, 43, 129–159.Claerbout, J. and F. Muir, 1973, Robust modeling with erratic data: Geo-physics, 38, 826–844.Claerbout, J. and D. Nichols, 1994, Spectral preconditioning: Technicalreport, Stanford Exploration Project.Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterative thresholdingalgorithm for linear inverse problems with a sparsity constraints: Comm.Pure Appl. Math., 57, 1413–1457.de Hoop, M. V. and N. Bleistein, 1997, Generalized Radon transform inver-63Chapter 2. Bibliographysions for reflectivity in anisotropic elastic media: Inverse Problems, 13,669–690.de Hoop, M. V. and S. Brandsberg-Dahl, 2000, Maslov asymptotic extensionof generalized Radon transform inversion in anisotropic elastic media: aleast-squares approach.: Inverse problems, 16, 519–562.Donoho, D. L., 1995, Nonlinear solution of linear inverse problems bywavelet-vaguelette decomposition: App. and Comp. Harmonic Analysis,2, 101–126.Donoho, D. L., M. Elad, and V. Temlyakov, 2006, Stable recovery of sparseovercomplete representations in the presence of noise: IEEE Trans. In-form. Theory, 52, 6–18.Elad, M., J. L. Starck, P. Querre, and D. L. Donoho, 2005, Simultane-ous cartoon and texture image inpainting using morphological componentanalysis (mca): Appl. Comput. Harmon. Anal., 19, 340–358.Fehmers, G. C. and C. F. Hocker, 2003, Fast structural interpretation withstructure-oriented filtering: Geophysics, 68, 1286–1293.Figueiredo, M. and R. Nowak, 2003, An em algorithm for wavelet-basedimage restoration: IEEE Trans. Image Processing.Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor nonunitary imaging operators: Geophysics, 69, 1017–1024.Hennenfent, G. and F. J. Herrmann, 2006, Seismic denoising with non-uniformly sampled curvelets: Comp. in Sci. and Eng., 8, 16–25.Herrmann, F. J., 2003, Multifractional splines: application to seismic imag-ing: Proceedings of SPIE Technical Conference on Wavelets: Applicationsin Signal and Image Processing X, 240–258.64Chapter 2. BibliographyHerrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Nonlinear primary-multiple separation with directional curvelet frames: Geoph. J. Int., 170,781–799.Hormander, L., 1987, The analysis of linear partial differential operators iii:Pseudo-differential operators.Hu, J., G. T. Schuster, and P. Valasek, 2001, Poststack migration deconvo-lution: Geophysics, 66, 939–952.Kuhl, H. and M. Sacchi, 2003, Least-squares wave-equation migration forAVP/AVA inversion: Geophysics, 68, 262.Lee, N. Y. and B. J. Lucier, 2001, Wavelet methods for inverting the Radontransform with noisy data: IEEE Trans. Image Processing, 10, 79–94.Nemeth., T., C. Wu, and G. T. Schuster, 1999, Least-squares migration ofincomplete reflection data: Geophysics, 64, 208–221.O’Brien, M. and S. Gray, 1997, Can we image beneath salt?: Leading Edge,15, 17–22.Paige, C. C. and M. A. Saunders, 1982, Lsqr: An algorithm for sparse linearequations and sparse least squares: ACM TOMS, 8, 43–71.Plessix, R. and W. Mulder, 2004, Frequency-domain finite differenceamplitude-preserving migration: Geoph. J. Int., 157, 975–987.Rickett, J. E., 2003, Illumination-based normalization for wave-equationdepth migration: Geophysics, 68, no. 4.Schertzer, O., 2003, Scale-space methods and regularization for denosingand inverse problems: Advances in Imaging and Electron Physics, 128,445–530.Starck, J. L., M. Elad, and D. L. Donoho, 2004, Redundant multiscale65Chapter 2. Bibliographytransforms and their application to morphological component separation:Advances in Imaging and Electron Physics, 132.Stolk, C. C., 2000, Microlocal analysis of a seismic linearized inverse prob-lem: Wave Motion, 32, 267–290.Stolk, C. C. and W. W. Symes, 2003, Smooth objective functionals forseismic velocity inversion: Inverse Problems, 19, 73–89.Symes, W. W., 2006a, Optimal scaling for reverse time migration: TechnicalReport TR 06-19, Department of Computational and Applied Mathemat-ics, Rice University, Houston, Texas, USA.——–, 2006b, Reverse time migration with optimal checkpointing: TechnicalReport 06-18, Department of Computational and Applied Mathematics,Rice University, Houston, Texas, USA.Symes, W. W. and J. Carazzone, 1991, Velocity inversion by differentialsemblance optimization: Geophysics, 56, 2061–2073.ten Kroode, A., D.-J. Smit, and A. Verdel, 1998, A microlocal analysis ofmigration: Wave Motion, 28.Tibshirani, R., 1997, Regression shrinkage and selection via the LASSO: J.Royal. Statist. Soc B., 58, 267–288.Tropp, J. A., 2006, Just relax: convex programming methods for identifyingsparse signals in noise: IEEE Trans. Inform. Theory, 52, 1030–1051.Tsaig, Y. and D. L. Donoho, 2006, Extensions of compressed sensing: Signalprocessing, 86, 549–571.Ying, L., L. Demanet, and E. J. Candes, 2005, 3D discrete curvelet trans-form: Presented at the Wavelets XI, SPIE.66Chapter 3Curvelet-based migrationpreconditioning and scaling3.1 IntroductionOver the years, extensive research has been done to reduce the computa-tional costs of (least-squares) seismic imaging. Improvements in this areaare particularly important during iterative least-squares migration, wherethe linear Born-scattering operator is inverted with iterative Lanczos meth-ods, such as LSQR (Paige and Saunders, 1982; De Roeck, 2002). Examplesof these methods can be found in the literature (see e.g. Nemeth et al.,1999; Chavent and Plessix, 1999; Hu et al., 2001; Kuhl and Sacchi, 2003; Yuet al., 2006).The most successful methods to reduce the cost of migration are the so-calledscalingmethodswheretheactionofthecompoundlinearizedmodeling-migration operator—known as the Hessian or normal operator—is replacedby a diagonal scaling in some domain, see e.g. contributions by ClaerboutA version of this chapter has been published. Herrmann, F.J., Brown, C.R., Erlangga,Y.A. and Moghaddam, P.P. (2009) Curvelet-based migration preconditioning and scaling.Geophysics, 74, pp. A41, 2009673.1. Introductionand Nichols (1994); Rickett (2003); Guitton (2004); Plessix and Mulder(2004), and more recently by Symes (2008) and Herrmann et al. (2008a).These methods vary in degree of sophistication with regard to the estimationof the diagonal, e.g. through migrated-image to remigrated-image match-ing. They also differ in the way the scaling is applied—i.e., by ’division’ inthe physical domain or via sparsity promotion in the curvelet domain, asreported recently by Herrmann et al. (2008a). During all these methods, im-aged amplitudes are restored by applying a scaling as a post-processing stepafter migration. Even though our approach and the one of Symes (2008)are similar, the curvelet-domain scaling has the advantage that it is able tohandle conflicting dips.In this paper, we take this line of research a step further by using theabove scaling argument to apply the proper preconditioning to the systemof equations involved in linearized Born scattering. By “proper” we mean apreconditioner with a computational overhead that justifies its improvementby the increase in convergence rate. Because the system is solved iteratively,the preconditioner does not need to be very accurate, avoiding unnecessaryextra computational overhead. Note that we use the term precondition-ing somewhat loosely compared to its formal definition in numerical linearalgebra where the solution is not changed with preconditioning (Calvetti,2007). Therefore, we also use this term to denote changes in the forwardmodel that favor least-squares inversion. To illustrate the improvements inthe images and in the convergence of least-squares migration, we considerthree levels of preconditioning. First, we correct the order of the normaloperator by introducing a left preconditioning, consisting of a fractional683.2. Problem formulationtime integration that corresponds to a scaling in the Fourier domain. Thisfirst level of preconditioning follows directly from earlier work on migration-amplitude recovery reported by Herrmann et al. (2008a) and Symes (2008).The next level of preconditioning consists of diagonal scaling in the physi-cal domain that compensates for spherical spreading of seismic waves. Asa final step, we also include a curvelet-domain scaling as part of the rightpreconditioning. This step corrects for the remaining amplitude errors thatvary spatially as a function of the reflector dip. We conclude by study-ing the performance of these different levels of preconditioning on the SEGAA’ salt model (O’Brien and Gray, 1996; Aminzadeh et al., 1997), usinga reverse-time ‘wave-equation’ migration code with optimal checkpointing(Symes, 2007). Because our main interest is to study the performance of ourpreconditioner, we work on data generated by the linearized Born-scatteringoperator instead of data generated by the solution of the full nonlinear for-ward model.3.2 Problem formulationDuring seismic imaging, the following system of equations needs to be solvedAx≈b, (3.1)where b = Ax0 is the multiple-free data with x0 the true reflector model, Athe linearized Born-scattering operator, and x the unknown model vector.Here, the symbol ≈ refers to “equality” in the least-squares sense. Con-693.2. Problem formulationtrary to many inverse problems, the matrix A, albeit extremely large, isreasonably well behaved and a migrated image can be obtained by applyingthe adjoint of A to the data vector—i.e., tildewidex = A∗b with A∗ the migrationoperator. The symbol ∗ denotes the adjoint, and tildewidex is the estimate (denotedby thetildewide) of the image.Unfortunately, the output of the above procedure, called migration, pro-duces erroneous results for the amplitudes of the imaged reflectors. Torestore these amplitudes, the least-squares solution to Equation 3.1 can beobtained as the solution of the linear systemA∗Ax≈A∗b, (3.2)with A∗A the normal or Hessian operator. Solutions to this system are notunique and correspond to solutions of the least-squares method—i.e.,tildewidexLS = argminx12bardblb−Axbardbl22, (3.3)which finds image vectors, x, that after modeling fit the data vector, b.However, the matrix A∗A is not invertible—i.e., it has zero or small singu-lar values that correspond to shadow zones. Because the system is large,we employ an iterative matrix-free solution method, such as LSQR (Paigeand Saunders, 1982). By limiting the number of iterations for this method,we control the energy of the solution vector tildewidex and we obtain a regularizedsolution that is equivalent to the solution of a damped least-squares prob-lem (Vogel, 2002). Even though these iterative methods converge relatively703.2. Problem formulationquickly (see e.g. Nemeth et al., 1999; Chavent and Plessix, 1999; Hu et al.,2001; Kuhl and Sacchi, 2003; Yu et al., 2006), the sheer size of the imagingproblem calls for a further reduction in the number of iterations.In a perfect world, with infinite computational resources, the ideal pre-conditioning for the system in Equation 3.1 corresponds toAM−1R u≈b, x := M−1R u, (3.4)with the right preconditioning matrix, MR := parenleftbigA∗Aparenrightbig1/2, given by the‘square-root’ of the normal operator. Here, the symbol := refers to ‘definedas’. In this ideal case, migration recovers the image vector x, exactly. Un-fortunately, in practice (albeit some recent exciting progress has been madeby Demanet and Ying, 2008, using discrete symbol calculus for smooth sym-bols, a development on which we intend to report in the future) the quantityparenleftbigA∗Aparenrightbig1/2 cannot be computed and we have to resort to appropriate approx-imations.Inthispaper, weproposeacombinationofleftandrightpreconditioning—i.e., we replace Equation 3.1 bybAbracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownrightM−1L AM−1R u≈bbbracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownrightM−1L b, x := M−1R u, (3.5)with M−1L the left preconditioning matrix. The migrated and least-squaresmigrated images are given by tildewidex = M−1R tildewideu, with tildewideu = hatwideA∗hatwideb, and by tildewidexLS =M−1R tildewideuLS, with tildewideuLS = argminubardblhatwideb− hatwideAubardbl2, respectively. Our precondi-tioners are derived from the following three observations: (i) under certain713.3. Preconditioningconditions—such as the high-frequency limit, smooth background velocitymodels, and the absence of turning waves (Stolk, 2000)—the normal oper-ator is in d dimensions a (d−1)-order pseudodifferential operator (ΨDO,see e.g. recent work by Herrmann et al., 2008a; Symes, 2008, and the refer-ences therein), (ii) migration amplitudes decay with depth due to sphericalspreading of seismic body waves, and (iii) zero-order ΨDO’s can be approx-imated by a diagonal scaling in the curvelet domain (see e.g. Herrmannet al., 2008a). These observations allow us to define a series of increasinglymore accurate approximations to the ‘square-root’ of the normal operator,leading to better and better preconditioners. Finally, we also argue thatusing curvelets will add a certain robustness to Gaussian noise and mod-eling errors, an observation substantiated by successful applications of thistransform in seismic data processing (Herrmann et al., 2008b; Wang et al.,2008).3.3 PreconditioningIn this section, we introduce different types of preconditioners based on theaforementioned observations. For each preconditioned system, we study themigrated images. Later, we study the convergence of the iterative solver.The examples are computed for the reflectivity and smooth velocity back-ground models plotted in Figure 3.1. To test our preconditioner, data with324 shots is generated using Equation 3.1. Each shot consists of 176 tracesof 6.4s and with a trace interval of 24m. The maximum offset of the data is4224m.723.3. Preconditioning3.3.1 Left preconditioning by fractional differentiationAs stated before, under aforementioned conditions the normal operator cor-responds to a (d−1)-order ΨDO. In 2-D image space, this operator corre-sponds to the leading-order behavior of the square root of Laplacian—i.e.,the action of parenleftbig∆parenrightbig1/2·in the physical domain, or to|ξ|·with ξ the wave vec-tor in the spatial Fourier domain. In data space, this action corresponds toa multiplication by |ω| in the temporal Fourier domain (Herrmann et al.,2008a). We compensate for this action by defining the following left precon-ditioning:M−1L := ∂−1/2|t| , (3.6)where ∂−1/2|t| ·:= F∗|ω|−1/2F·with F the Fourier transform and F∗ = F−1 itsinverse. We define the level I preconditioner with M−1L as above and M−1R =I(whichimpliesx = u). Comparisonofthemigratedimagesbeforeandafterleft preconditioning (cf. Figures 3.1(c) and 3.2(a)) shows that the imprint ofthe Laplacian is removed—i.e., some of the low-frequency content is restored.However, the migrated image still contains dimming of the amplitudes. Notethat this left preconditioning corresponds to the solution of the scaled least-squares problem: argminx 12bardblhatwideb−hatwideAxbardbl22 = argminx 12bardblM−1L (b−Ax)bardbl22.3.3.2 Right preconditioning by scaling in the physicaldomainTo further correct the amplitudes, we propose to apply a scaling to com-pensate for the leading-order amplitude decay. This decay is linear becausethe reflected waves travel from the source down to the reflector, experi-733.3. Preconditioningencing an amplitude decay proportional to the square-root of the reflectordepth (in 2-D). We correct this linear amplitude decay by defining the rightpreconditioning matrix:M−1R = Dz := diagparenleftbigzparenrightbig12, (3.7)where zi = i∆z, i = 1···nz, with ∆z the depth sample interval and nz thenumber of samples. Combined with the left preconditioner M−1L , we callthis the level II preconditioner. The results for the migrated image in Fig-ure 3.2(b) now show further improvement. However, amplitude variationsremain, e.g., along the major horizontal reflector just above 3500m.3.3.3 Right preconditioning by scaling in the curveletdomainAfter applying the left preconditioning, the Hessian can be modeled by azero-order ΨDO whose action corresponds to that of a nonstationary dipfilter—i.e., we haveparenleftbigΨfparenrightbig(x)similarequalintegraldisplayξ∈Rdejξ·xa(x,ξ)ˆf(ξ)dξ, (3.8)with Ψ the Hessian of the preconditioned modeling operator and a(x,ξ)a space- and spatial-frequency dependent filter known as the symbol. Weuse the symbolsimilarequalto indicate high-frequency approximation and absence ofturning waves. Following ideas that go back to Taylor (1981), a pseudodiffer-ential operator can approximately be diagonalized by linear combinations of743.3. Preconditioninglocalized oscillatory functions, such as curvelets (Cand`es et al., 2006), waveatoms (Demanet and Ying, 2007), and local Fourier bases (Meyer, 1992).In this paper, we use recent results by Herrmann et al. (2008a) who usecurvelets because they have the additional advantage of being sparse on themodel. The action of the ΨDO can, after discretization, be approximatedby a scaling in the curvelet domain—i.e., we have the following approximateidentityΨr≈C∗D2ΨCr, D2Ψ := diagparenleftbigd2parenrightbig, (3.9)which is accurate for a reference vector r close enough to the actual im-age. In this expression, the matrices C and C∗ represent the 2-D discretecurvelet transform (see e.g. Cand`es et al., 2006) for which the adjoint equalsthe pseudoinverse—i.e., we have C∗C = I with I the identity matrix. Thereciprocal of the curvelet-domain scaling coefficients, d−2, is found by a rem-igrated image-to-image matched-filtering procedure that involves the refer-ence vector, typically derived from a conventional migrated image, and theremigrated reference vector. Hence, the cost of forming Equation 3.9 is ap-proximately one modeling and one migration (for further details refer toHerrmann et al., 2008a,c). As we will show below, these additional costs arewell offset by the increase in convergence.By including the above approximation, we define the right precondition-ing matrix asM−1R = DzC∗D−1Ψ , (3.10)which compensates for the remaining amplitude errors. Because there isno need for high accuracy, Equation 3.9 is a good approximation whose ac-753.4. Convergence of least-squares migrationcuracy increases with frequency (Herrmann et al., 2008a). Note, however,that this preconditioner by virtue of the underlying assumptions will nothelp in the removal of nonlocal imaging artifacts. We call the combinationof Equations 3.6 and 3.10, the level III preconditioner. The image obtainedwith this system is plotted in Figure 3.2(c) and shows, as expected, fur-ther improvement in the amplitudes of the migrated image. Despite thisimprovement, image artifacts and amplitude errors are still present and canbe attributed to the approximation in Equation 3.9 and to the fact thatmigration does not correspond to inversion—i.e., A∗Anegationslash= I.3.4 Convergence of least-squares migrationEventhoughthedifferentpreconditioningoperatorsdefinedsofarleadtoim-provements, problems remain in achieving high-fidelity images. As reportedin the literature, the image quality can be further improved by replacingmigration with least-squares migration where the (preconditioned) scatter-ing matrix is inverted iteratively. After a limited number of iterations, thismethod inverts the (preconditioned) system of equations approximately.The performance of iterative solvers depends on certain properties ofthe (preconditioned) matrix A, which include its condition number (ratioof the largest to smallest singular value) and the clustering of the singu-lar values (see e.g. De Roeck, 2002, where these quantities are discussedfor a relatively small Kirchhoff-based imaging problem). Because LSQRminimizes the residual during each iteration (cf. Equation 3.3), studying itsprogress as a function of the number of iterations gives us some way to gauge763.4. Convergence of least-squares migrationthe performance. For this purpose, we introduce the normalized log-basedleast-squares residual µk = 20logbardblhatwideAuk −hatwidebbardbl2/bardblhatwidebbardbl2, with uk the solutionof the (preconditioned) system after k iterations and with u0 = 0. Follow-ing De Roeck (2002), we also track the log-based least-squares model-spaceresidual—i.e., νk = 20logbardblhatwideA∗parenleftbighatwideAuk−hatwidebparenrightbigbardbl2/bardblhatwideA∗hatwidebbardbl2. Note that in practice,these latter residuals are never computed because they are not a by-productof LSQR, and require an additional matrix-vector multiply per iteration.Contrary to data-space residuals, which possibly contain unmodeled com-ponents that may not be in the range of the modeling operator, model-spaceresiduals typically converge to zero. Both quantities are used to empiricallyestablish the performance of the different levels of preconditioning. The re-sults of this exercise are summarized in Figure 3.3, which plots the decay ofµk and νk as a function of the number of iterations k. To account for theoverhead, plots for the level III preconditioner are offset by one iteration.This is justified because the cost per iteration of evaluating the curvelettransform, and its inverse, is negligible compared to the cost of migrationand demigration.As we move from a single left preconditioner, towards left and right pre-conditioners, the data residuals decay faster with a significant improvementobtained by the curvelet-domain scaling. For instance, the residual after10 iterations for level II preconditioning (fractional integration and depthweighting) is attained by only 5 iterations of level III preconditioning (in-cluding curvelet-domain scaling), whereas the result after 10 iterations isapproximately 2dB better. Including the computational overhead, level IIIpreconditioning gains 8dB in four iterations. Even though the picture for773.5. Extensionsthe model-space residual is less clear, there is a similar trend for the level IIIpreconditioning. For instance, after 10 iterations, we have an improvementof approximately 4.5dB.The improvements in convergence for the preconditioned system are alsoreflected in the least-squares migrated results included in Figure 3.4. Com-paring these images shows a clear enhancement for the preconditioned sys-tem, plotted in Figure 3.4(b), over the least-squares result obtained withoutpreconditioning. Moreover, juxtaposingthe preconditioned least-squares im-age with the solution for the preconditioned migrated image after one itera-tion depicted in Figure 3.2(c) shows a significant enhancement of the overallamplitudes and frequency content (cf. Figure 3.1(a) and 3.4(b)).The bottom line is that each iteration takes 40 and 180 minutes for themodeling and migration on an IBM eServer with 52 processors at 2.2GHz.Aside from one additional modeling and migration, the computation of thediagonal estimation takes approximately 90 minutes on a single CPU. Thecost of the 2-D curvelet-transforms part of the preconditioning is less thanone minute and is negligible compared to the modeling and migration costs.3.5 ExtensionsThe preconditioning methodology presented in this letter constitutes a firststep towards a concerted effort to formulate seismic imaging and inversionas an optimization problem. This type of formulation allows us to createhigh-fidelity images and to make progress towards full-waveform inversion.Having access to appropriate preconditioners is instrumental for this purpose783.5. Extensionsbecause it makes the computations tractable. We envisage the followingextensions of our work:• generalization to 3-D, which entails a different power for the fractionalintegration, a different spherical-spreading correction, and a 3-D for-mulation of our curvelet-domain matched filter. Moreover, in 3-D thetheory of approximating the Hessian by a ΨDO is less well developed.• replacement of the level II physical-space preconditioning by moreaccurate source Green’s function illumination corrections (see e.g.,Plessix and Mulder, 2004).• inclusion of density variations, or in case of elastic wave propagation,the inclusion of variations in the elastic moduli. This extension re-quires a multi-parameter formulation.• regularization by curvelet-domain sparsity promotion, replacing Equa-tion 3.3 bytildewideulscript1 = argminubardblubardbl1 subject to bardblhatwideAu−hatwidebbardbl2≤epsilon1 (3.11)with epsilon1 a noise-dependent parameter, and xlscript1 := M−1R tildewideulscript1. To ben-efit from our preconditioning, Equation 3.11 requires a solver withNewton-type steps, as opposed to most lscript1-norm solvers that are basedon projected gradients. The advantage of this formulation is that ituses curvelet-domain sparsity, which has proven to be a particularlypowerful prior (Wang and Sacchi, 2007; Herrmann et al., 2008b,a; Hen-nenfent et al., 2008). We will report on the solution of Equation 3.11793.6. Conclusionselsewhere.• incorporationofourpreconditionerinsparsity-promotingfull-waveforminversion. This problem requires the solution of the unconstrained op-timization problemminz 12bardblhatwideb−hatwideF[z]bardbl22 +λbardblzbardbl1, (3.12)where z is the curvelet representation of the model, λ a Lagrangemultiplierbalancingtheresidualenergyandlscript1-normpenaltyterm, andhatwideF[z] the nonlinear forward map that links the curvelet-domain modelto data. Again, the solution of this optimization problem requires asophisticated solver that can benefit from our preconditioner.3.6 ConclusionsBecause of the size of the seismic imaging problem, preconditioning of least-squares migration is an elusive topic, where traditional approaches fromnumerical linear algebra have not yet found their way. Lack of direct accessto the matrices involved and the cost of evaluating matrix-free implementa-tions of the operators are both to blame. Nonetheless, the first few iterationsof the LSQR algorithm of least-squares migration are known to make sig-nificant progress towards the solution. Unfortunately, even for this limitednumber of iterations, the computational costs are often still prohibitivelylarge for practical problems. The method presented in this paper partlyresolves this issue through a combination of left and right preconditioning803.7. Acknowledgmentstogether with a curvelet-domain scaling. Inclusion of the latter proved par-ticularly important because it restores the amplitudes and leads to fasterconvergence, at a relatively small computational overhead. Aside from thistangible reduction in computational costs of roughly 50%, the use of curveletframes opens the enticing perspective to use lscript1-norm regularization to im-prove the quality of images. Preconditioning also plays a pivotal role inmaking this approach numerically feasible and may extend to a solutionof the full-waveform inversion problem with sparsity promotion. Both ap-proaches are justified by ample evidence that curvelets are sparse on themodel.3.7 AcknowledgmentsThe authors would like to thank C. C. Stolk and W. W. Symes for their inputon migration-amplitude recovery and for the use of the reverse-time migra-tion code. We also would like to thank Tamas Nemeth and Michael Friedlan-derforinsightfuldiscussionsandtheauthorsofCurveLab(www.curvelet.org).The examples were prepared with Madagascar (rsf.sf.net), supplemented bySLIMpy operator overloading, developed by S. Ross-Ross. This work was inpart financially supported by the NSERC Discovery (22R81254) and CRDGrants DNOISE (334810-05) of F.J.H. and was in part carried out withinthe SINBAD project with support, secured through ITF, from BG Group,BP, Chevron, ExxonMobil and Shell. We finally thank the anonymous re-viewer and W. W. Symes for their constructive comments and suggestionsthat greatly improved our paper.813.7. Acknowledgments(a)(b)(c)Figure 3.1: The SEG/EAGE AAprime salt model, (a) reflectivity defined bythe high-pass filtered velocity model, (b) smoothed velocity model, (c) themigrated image according to Equation 3.1. This image suffers from dete-riorated amplitudes, especially under the high-velocity salt and for steepreflectors and faults.823.7. Acknowledgments(a)(b)(c)Figure 3.2: Migrated images for different levels of preconditioning, (a) resultfor left preconditioning (level I, cf. Equation 3.6), (b) result for left-right(including depth-correction) preconditioning (level II, cf. Equation 3.7), (c)the same but now including curvelet-domain scaling (level III, cf. Equation3.10).833.7. Acknowledgments0 2 4 6 8 10Iterations121086420µk(dB)NoneLevel ILevel IILevel III(a)0 2 4 6 8 10Iterations181614121086420νk(dB)NoneLevel ILevel IILevel III(b)Figure 3.3: Residual decays for different levels of preconditioning. Thedotted blue lines corresponds to least-squares migration without precondi-tioning, the dash-dotted lines to level I preconditioning, the dashed blacklines to level II preconditioning, and the red solid lines to level III precon-ditioning. This is offset by one iteration to account for the overhead, (a)plot for the decay of the data-space normalized residues µk as a function ofthe number of LSQR iterations, (b) the same but now for the model-spacenormalized residuals νk.843.7. Acknowledgments(a)(b)Figure 3.4: Least-squares migration without and with preconditioning jux-taposed with the original reflectivity (in light blue), (a) least-squares mi-grated image, (b) least-squares image with level III preconditioning. Noticethe improvement in the recovered reflectivity.85BibliographyAminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D Salt and Overthrust Model.SEG/EAGE 3-D Modeling Series, No. 1: Society of Exploration Geophysi-cists, Tulsa.Calvetti, D., 2007, Preconditioned iterative methods for linear discrete ill-posed problems from a Bayesian inversion perspective: Journal of Com-putational and Applied Mathematics, 198, 378–395.Cand`es, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006, Fast discretecurvelet transforms: SIAM Multiscale Modeling and Simulation, 5, 861–899.Chavent, G. and R. E. Plessix, 1999, An optimal true-amplitude least-squares prestack depth-migration operator: Geophysics, 64, 508–515.Claerbout, J. and D. Nichols, 1994, Spectral preconditioning: TechnicalReport SEP-82, Stanford Exploration Project.De Roeck, Y., 2002, Sparse linear algebra and geophysical migration: Areview of direct and iterative methods: Numerical Algorithms, 29, 283 –322.Demanet, L. and L. Ying, 2007, Wave atoms and sparsity of oscillatorypatterns: Journal of Applied and Computational Harmonic Analysis, 23,368–387.86Chapter 3. Bibliography——–, 2008, Discrete symbol calculus. (Submitted for publication. Availableat http://math.stanford.edu/ laurent/papers/DSC.pdf).Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor nonunitary imaging operators: Geophysics, 69, 1017–1024.Hennenfent, G., E. van den Berg, M. P. Friedlander, and F. J. Herrmann,2008, New insights into one-norm solvers from the Pareto curve: Geo-physics, 73, no. 4, A23–A26.Herrmann, F. J., P. P. Moghaddam, and C. C. Stolk, 2008a, Sparsity- andcontinuity-promoting seismic imaging with curvelet frames: Journal ofApplied and Computational Harmonic Analysis, 24, 150–173.Herrmann, F. J., D. Wang, G. Hennenfent, and P. P. Moghaddam, 2008b,Curvelet-based seismic data processing: a multiscale and nonlinear ap-proach: Geophysics, 73, no. 1, A1–A5.Herrmann, F. J., D. Wang, and D. J. E. Verschuur, 2008c, Adaptive curvelet-domain primary-multiple separation: Geophysics, 73, A17–A21.Hu, J., G. T. Schuster, and P. A. Valasek, 2001, Poststack migration decon-volution: Geophysics, 66, 939–952.Kuhl, H. and M. D. Sacchi, 2003, Least-squares wave-equation migration forAVP/AVA inversion: Geophysics, 68, 262–273.Meyer, Y., 1992, Wavelets and operators: Cambridge University Press.Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration ofincomplete reflection data: 64, 208–221.O’Brien, M. and S. Gray, 1996, Can we image beneath salt?: The LeadingEdge, 15, 17–22.Paige, C. C. and M. A. Saunders, 1982, LSQR: An algorithm for sparse87Chapter 3. Bibliographylinear equations and sparse least squares: ACM TOMS, 8, 43–71.Plessix, R. and W. Mulder, 2004, Frequency-domain finite differenceamplitude-preserving migration: Geophysical Journal International, 157,975–987.Rickett, J. E., 2003, Illumination-based normalization for wave-equationdepth migration: Geophysics, 68, 1371–1379.Stolk, C. C., 2000, Microlocal analysis of a seismic linearized inverse prob-lem: Wave Motion, 32, 267–290.Symes, W. W., 2007, Reverse time migration with optimal checkpointing:Geophysics, 72, SM213–SM221.——–, 2008, Approximate linearized inversion by optimal scaling of prestackdepth migration: Geophysics, 73, R23–R35.Taylor, M., 1981, Pseudodifferential operators: Princeton University Press.Vogel, C., 2002, Computational Methods for Inverse Problems: SIAM.Wang, D., R. Saab, O. Yilmaz, and F. J. Herrmann, 2008, Bayesian wave-field separation by transform-domain sparsity promotion: Geophysics, 73,A33–A38.Wang, J. and M. D. Sacchi, 2007, High-resolution wave-equation amplitude-variation-with-ray-parameter (AVP) imaging with sparseness constraints:Geophysics, 72, S11–S18.Yu, J., J. Hu, G. T. Schuster, and R. Estill, 2006, Prestack migration de-convolution: Geophysics, 71, S53–S62.88Chapter 4Curvelet-based seismic dataprocessing4.1 IntroductionInthisletter, wedemonstratethatthediscretecurvelettransform(Cand`eset al., 2006a; Hennenfent and Herrmann, 2006b) can be used to reconstructseismic data from incomplete measurements, to separate primaries and mul-tiples and to restore migration amplitudes. The crux of the method liesin the combination of the curvelet transform, which attains a fast decayfor the magnitude-sorted curvelet coefficients, with a sparsity promotingprogram. By themselves sparsity-promoting programs are not new to thegeosciences (Sacchi et al., 1998). However, sparsity promotion with thecurvelet transform is new. The curvelet transform’s unparalleled ability todetect wavefront-like events that are locally linear and coherent means it isparticularly well suited to seismic data problems. In this paper, we showexamples including data regularization (Hennenfent and Herrmann, 2006a,A version of this chapter has been published. Herrmann, F.J., Wang, D., Hennenfent,G. and Moghaddam, P.P. (2008) Curvelet-based seismic data processing: a multiscale andnonlinear approach.Geophysics, 73(1):A1-A5, 2008.894.2. Curvelets2007a), primary-multiple separation (Herrmann et al., 2007) and migration-amplitude recovery (Herrmann et al., 2008). Application of this formalismto wavefield extrapolation is presented elsewhere (Lin and Herrmann, 2007).4.2 CurveletsCurvelets are localized ’little plane-waves’ (see Hennenfent and Her-rmann, 2006b, and the on-line ancillary material for an introduction onthis topic) that are oscillatory in one direction and smooth in the otherdirection(s). They are multiscale and multi-directional. Curvelets have ananisotropic shape – they obey the so-called parabolic scaling relationship,yielding a ”width” proportional to square of the ”length” for the supportof curvelets in the physical domain. This anisotropic scaling is necessaryto detect wavefronts and explains their high compression rates on seismicdata and images (Cand`es et al., 2006a; Herrmann et al., 2008), as long asthese datasets can be represented as functions with events on piece-wisetwice differentiable curves. Then, the events become linear at the fine scalesjustifying an approximation by the linearly shaped curvelets. Even seismicdata with caustics, pinch-outs, faults or strong amplitude variations fit thismodel, which amounts to a preservation of the sparsity attained by curvelets.Curvelets represent a specific tiling of the 2-D/3-D frequency domaininto strictly localized wedges. Because the directional sampling increasesevery-other scale doubling, curvelets become more anisotropic at finer scales.Curvelets compose multi-D data according to f = CTCf with C and CTthe forward and inverse discrete curvelet transform matrices (defined by the904.3. Common problem formulation by sparsity-promoting inversionfast discrete curvelet transform, FDCT, with wrapping, a type of periodicextenstion, see Cand`es et al., 2006a; Ying et al., 2005). The symbol Trepresents the transpose, which is equivalent to the inverse for this choiceof curvelet transform. This transform has a moderate redundancy (a factorof roughly 8 in 2-D and 24 in 3-D) and a computational complexity ofO(nlogn) with n the length of f. Even though CTC = I, with I theidentity matrix, the converse is not true, i.e., CCT negationslash= I.4.3 Common problem formulation bysparsity-promoting inversionOursolutionstrategyisbuiltonthepremisethatseismicdataandimageshave a sparse representation, x0, in the curvelet domain. To exploit thisproperty, our forward model readsy = Ax0 +n (4.1)with y a vector of noisy and possibly incomplete measurements; A themodeling matrix that includes CT; and n, a zero-centered white Gaussiannoise. Because of the redundancy of C and/or the incompleteness of thedata, the matrix A can not readily be inverted. However, as long as thedata, y, permits a sparse vector, x0, the matrix, A, can be inverted by asparsity-promoting program (Cand`es et al., 2006b; Donoho, 2006):914.3. Common problem formulation by sparsity-promoting inversionPepsilon1 :tildewidex = argminxbardblxbardbl1 s.t. bardblAx−ybardbl2≤epsilon1tildewidef = STtildewidex(4.2)in which epsilon1 is a noise-dependent tolerance level, ST the inverse transform andtildewidef the solution calculated from the vector tildewidex (the symbol tildewide denotes a vectorobtained by nonlinear optimization) minimizing Pepsilon1. The difference betweentildewidex and x0 is proportional to the noise level.Nonlinear programs Pepsilon1 are not new to seismic data processing as in spikydeconvolution (Taylor et al., 1979; Santosa and Symes, 1986) and Fouriertransform-based interpolation (Sacchi et al., 1998). The curvelets’ high com-pression rate makes the nonlinear program Pepsilon1 perform well when CT is in-cluded in the modeling operator. Despite its large-scale and nonlinearity,the solution of the convex problem Pepsilon1 can be approximated with a limited(< 250) number of iterations of a threshold-based cooling method derivedfrom work by Figueiredo and Nowak (2003); Daubechies et al. (2004); Eladet al. (2005). At each iteration the descent update (x←x+ATparenleftbigy−Axparenrightbig),minimizing the quadratic part of Equation 4.2, is followed by a soft thresh-olding (x← Tλ(x) with Tλ(x) := sgn(x)·max(0,|x|−|λ|)) for decreasingthreshold levels λ. This soft thresholding on the entries of the unknowncurvelet vector captures the sparsity and the cooling, which speeds up thealgorithm, allows additional coefficients to fit the data.924.4. Seismic data recovery4.4 Seismic data recoveryThe reconstruction of seismic wavefields from regularly-sampled datawith missing traces is a setting where a curvelet-based method will performwell. As with other transform-based methods, sparsity is used to reconstructthe wavefield by solving Pepsilon1.It is also shown that the recovery performancecan be increased when information on the major primary arrivals is includedin the modeling operator.4.4.1 Curvelet-based recoveryThe reconstruction of seismic wavefields from incomplete data corre-sponds to the inversion of the picking operator R. This operator modelsmissing data by inserting zero traces at source-receiver locations where datais missing passing recorded traces unchanged. The task of the recovery is toundo this operation by filling in the zero traces. Since seismic data is sparsein the curvelet domain, the missing data can be recovered by compoundingthe picking operator with the curvelet modeling operator, i.e., A := RCT.With this definition for the modeling operator, solving Pepsilon1 corresponds toseeking the sparsest curvelet vector whose inverse curvelet transform, fol-lowed by the picking, matches the data at the nonzero traces. Applyingthe inverse transform (with S := C in Pepsilon1) gives the interpolated data. Fordetails on the conditions that determine successful recovery, refer to Hen-nenfent and Herrmann (2007a,b) and Herrmann and Hennenfent (2007).An example of curvelet-based recovery is presented in Figure 4.1 whichshows the results of decimating, and then reconstructing, a seismic dataset.934.4. Seismic data recoveryThe original shot and receiver spacings were 25m, and 80% of the traces werethrown out at random (see Figure 4.1(b)). Comparing the ’ground truth’in Figure 4.1(a) with the recovered data in Figure 4.1(c) shows a successfulrecovery in case the high-frequencies are removed. Aside from sparsity in thecurvelet domain, no prior information was used during the recovery, whichis quite remarkable. Part of the explanation lies in the curvelet’s ability tolocally exploit the 3-D geometry of the data and this suggests why curveletsare successful for complex datasets where other methods may fail.4.4.2 Focused recoveryIn practice, additional information on the to-be-recovered wavefield is of-ten available. For instance, one may have access to the predominant primaryarrivals or to the velocity model. In that case, the recently introduced fo-cal transform (Berkhout and Verschuur, 2006), which ’deconvolves’ the datawith an estimate of the primaries, incorporates this additional informationinto the recovery process. Application of this primary operator, ∆P, addsa wavefield interaction with the surface, mapping primaries to first-ordersurface-related multiples (Verschuur and Berkhout, 1997; Herrmann, 2007).Inversion of this operator, strips the data off one interaction with the surface,focusing primary energy to (directional) sources. This focusing correpondsto a collapse of the 3-D primary events to an approximate line source whichhas a sparser representation in the curvelet domain.By compounding the non-adaptive, data-independent, curvelet trans-form with the data-adaptive focal transform, i.e., A := R∆PCT, the re-covery can be improved by solving Pepsilon1. The solution of Pepsilon1 now entails944.4. Seismic data recovery(a) (b)(c) (d)Figure 4.1: Comparison between 3-D curvelet-based recovery by sparsity-promoting inversion with and without focusing, (a)fully sampled real SAGAdata shot gather, (b) randomly subsampled shot gather from a 3-D datavolume with 80% of the traces missing in the receiver and shot directions,(c) curvelet-based recovery, (d) curvelet-based recovery with focusing. No-tice the improvement (denoted by the arrows) from the focusing with theprimary operator.954.5. Seismic signal separationthe inversion of ∆P, yielding the sparsest set of curvelet coefficients thatmatches the incomplete data when ’convolved’ with the primaries. Apply-ing the inverse curvelet transform, followed by ’convolution’ with ∆P yieldsthe interpolation, i.e. ST := ∆PCT. Comparing the curvelet recovery withthe focused curvelet recovery (Figure 4.1(c) and 4.1(d)) shows an overallimprovement in the recovered details.4.5 Seismic signal separationPredictive multiple suppression involves two steps, namely multiple pre-diction and primary-multiple separation. In practice, the second step ap-pears difficult and adaptive least-squares lscript2-matched-filtering techniques areknown to lead to residual multiple energy, high frequency jitter and de-terioration of the primaries (Herrmann et al., 2007). By employing thecurvelet’s ability to detect wavefronts with conflicting dips (e.g. caustics),a non-adaptive, independent of the total data, separation scheme can bedefined that is robust with respect to moderate errors in the multiple pre-diction. The nonlinear program, Pepsilon1, with y defined by the total data, canbe adapted to separate multiples from primaries by replacing the lscript1 normby a weighted lscript1 norm, i.e., bardblxbardbl1 mapsto→bardblxbardbl1,w = summationtextµ|wµxµ| with µ runningover all curvelets and w a vector with positive weights. By defining theseweights proportional to the magnitude of the curvelet coefficients of the 2-DSRME-predicted multiples, the solution of Pepsilon1 with A := CT removes mul-tiples. Primaries and multiples naturally separate in the curvelet domainand the weighting further promotes this separation while solving Pepsilon1. The964.5. Seismic signal separationweights that are fixed during the optimization penalize the entries in thecurvelet vector for which the predicted multiples are significant. The em-phasis on the weights versus the data misfit (the proportionality constant) isuser defined. The estimate for the primaries is obtained by inverse curvelettransforming the curvelet vector that minimizes Pepsilon1 for the weighted lscript1 norm(A = ST := CT).Figure 4.2 shows an example of 3-D curvelet-based primary-multipleseparation of a North Sea dataset with the weights set according to thecurvelet-domain magnitudes of the SRME-predicted multiples multiplied by1.25. Comparison between the estimates for the primaries from adaptivesubtraction by lscript2-matched filtering (Verschuur and Berkhout, 1997) andfrom our nonlinear and non-adaptive curvelet-based separation shows animprovement in (i) the elimination of the focused multiple energy belowshot location 1000m, induced by out-of-plane scattering due to small 3-Dvariations in the multiple-generating reflectors and (ii) an overall improvedcontinuity and noise reduction. This example demonstrates that the multi-scale and multi-angular curvelet domain can be used to separate primariesand multiples given an inaccurate prediction for the multiples. However,the separation goes at the expense of a moderate loss of primary energywhich compares favorably compared to the loss associated with lscript2-matchedfiltering (see also Herrmann et al., 2007).974.5. Seismic signal separation(a) (b)(c) (d)Figure 4.2: 3-D Primary-multiple separation with Pepsilon1 for the SAGA dataset,(a) near-offset section including multiples, (b) the SRME-predicted multi-ples, (c) the estimated primaries according to lscript2-matched filtering, (d) theestimated primaries obtained with Pepsilon1. Notice the improvement, in areaswith small 3-D effects (ellipsoid) and residual multiples.984.6. Migration-amplitude recovery4.6 Migration-amplitude recoveryRestoring migration amplitudes is another area where curvelets can beshown to play an important role. In this application, the purpose is toreplace computationally expensive amplitude recovery methods, such asleast-squares migration (Nemeth et al., 1999; Kuhl and Sacchi, 2003), byan amplitude scaling (Guitton, 2004). This scaling can be calculated froma demigrated-migrated reference vector close to the actual reflectivity.In order to exploit curvelet sparsity, we propose to scale in the curveletdomain. This choice seems natural because migrated images suffer fromspatially varying and dip-dependent amplitude deterioration that can beaccommodated by curvelets. The advantages of this approach are manifoldand include (i) a correct handling of reflectors with conflicting dips and (ii)a stable curvelet sparsity-promoting inversion of the diagonal that restoresthe amplitudes and removes the clutter by exploiting curvelet sparsity onthe model.The method is based on the approximate identity: KTKr≈CTDrCrwith K and KT the demigration, migration operators and Dr a reference-model specific scaling (Herrmann et al., 2008). By defining the modelingmatrix as A := CT√Dr, Pepsilon1 can be used to recover the migration am-plitudes from the migrated image. Possible spurious side-band effects anderroneously detected curvelets (Cand`es and Guo, 2002) are removed by sup-plementing the lscript1 norm in Pepsilon1 with an anisotropic diffusion norm (Fehmersand H¨ocker, 2003). This norm enhances the continuity along the imagedreflectors and removes spurious artifacts.994.7. Discussion and conclusionsResults for the SEG AA’ dataset (O’Brien and Gray, 1996; Aminzadehet al., 1997) are summarized in Figure 4.3. These results are obtained witha reverse-time ’wave-equation’ finite-difference migration code. To illustratethe recovery performance, idealized seismic data is generated by demigra-tion, followed by adding white Gaussian noise, yielding a signal-to-noise ratio(SNR) of only 3dB. This data is subsequently migrated and used as input.Despite the poor SNR, the image in Figure 4.3(a) contains most reflectors,which can be explained by the redundancy of the data, the migration op-erator’s sophistication (diffractions at the bottom of the salt are handledcorrectly) and the perfect ’match’ between the demigration and migrationoperators. However, the noise gives rise to clutter and there is dimming ofthe amplitudes, in particular for steep dips under the salt. Nonlinear recov-ery removes most of this clutter and more importantly the amplitudes forthe sub-salt steep-dipping events are mostly restored. This idealized exam-ple shows how curvelets can be used to recover the image amplitudes. Aslong as the background velocity model is sufficiently smooth and the reflec-tivity sufficiently sparse, this recovery method can be expected to performwell even for more complex images.4.7 Discussion and conclusionsThe presented examples show that problems in data acquisition andimaging can be solved with a common problem formulation during whichsparsity in the curvelet domain is promoted. For curved wavefront-like fea-tures that oscillate in one direction and that are smooth in the other di-1004.7. Discussion and conclusions(a)(b)Figure 4.3: Image amplitude recovery for a migrated image calculated fromnoisy data (SNR 3dB), (a) image with clutter, (b) image after nonlinearrecovery. The clearly visible non-stationary noise in (a) is mostly removedduring the recovery while the amplitudes are also restored. Steeply dippingreflectors (denoted by the arrows) under the salt are also well recovered.1014.8. Acknowledgmentsrection(s), curvelets attain high compression rates while other transformsdo not necessarily achieve sparsity for these geometries. Seismic imagesof sedimentary basins and seismic wave arrivals in the data both behavein this fashion, so that curvelets are particularly valuable for compression.It is this compression that underlies the success of our sparsity promotingformulation. First, we showed on real data that missing data can be re-covered by solving a nonlinear optimization problem where the data misfitand the lscript1-norm on the curvelet coefficients are simultaneously minimized.This recovery is improved further with a combined curvelet-focal transform.Sparsity also proved essential during the primary-multiple separation. Inthis case, it leads to a form of decorrelation of primaries and multiples, re-ducing the probability of having large overlapping curvelet entries betweenthese different events. Finally, the sparsity of curvelets on the image itselfwas exploited to recover the migration amplitudes of the synthetic subsaltimaging example. Through these three examples, the successful applicationof curvelets, enhanced with sparsity-promoting inversion, opens new per-spectives on seismic data processing and imaging. The ability of curveletsto detect wavefront-like features is key to our success and opens an excitingnew outlook towards future developments in exploration seismology.4.8 AcknowledgmentsThe authors would like to thank D.J. Verschuur and C. Stolk for theirinput in the primary-multiple separation and migration-amplitude recovery.We also would like to thank the authors of CurveLab (www.curvelet.org)1024.8. Acknowledgmentsand W. Symes for his reverse-time migration code. The examples wereprepared with Madagascar (rsf.sf.net), supplemented by SLIMpy operatoroverloading, developed by S. Ross Ross. Norsk Hydro is thanked for the fielddataset. M. O’Brien, S. Gray and J. Dellinger are thanked for the SEG AA’data. This work was in part financially supported by the NSERC Discovery(22R81254) and CRD Grants DNOISE (334810-05) of F.J.H. and was carriedout as part of the SINBAD project with support, secured through ITF, fromBG Group, BP, Chevron, ExxonMobil and Shell. Finally, the authors wouldlike to thank the anonymous reviewers whose constructive comments helpedimprove this letter.103BibliographyAminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D Salt and Overthrust Model.SEG/EAGE 3-D Modeling Series, No. 1: Society of Exploration Geophysi-cists, Tulsa.Berkhout, A. J. and D. J. Verschuur, 2006, Focal transformation, an imagingconcept for signal restoration and noise removal: Geophysics, 71, A55–A59.Cand`es, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2006a, Fast discretecurvelet transforms: SIAM Multiscale Modeling and Simulation, 5, 861–899.Cand`es, E. J. and F. Guo, 2002, New multiscale transforms, minimum totalvariation synthesis: Applications to edge-preserving image reconstruction:Signal Processing, 82, 1519–1543.Cand`es, E. J., J. K. Romberg, and T. Tao, 2006b, Stable signal recoveryfrom incomplete and inaccurate measurements: Communications on Pureand Applied Mathematics, 59, 1207–1223.Daubechies, I., M. Defrise, and C. De Mol, 2004, An iterative thresholdingalgorithm for linear inverse problems with a sparsity constraints: Com-munications on Pure and Applied Mathematics, 57, 1413–1457.Donoho, D. L., 2006, Compressed sensing: IEEE Transactions on Informa-104Chapter 4. Bibliographytion Theory, 52, 1289–1306.Elad, M., J. L. Starck, P. Querre, and D. L. Donoho, 2005, SimultaneousCartoon and Texture Image Inpainting using Morphological ComponentAnalysis (MCA): Journal of Applied and Computational Harmonic Anal-ysis, 19, 340–358.Fehmers, G. C. and C. F. W. H¨ocker, 2003, Fast structural interpretationwith structure-oriented filtering: Geophysics, 68, 1286–1293.Figueiredo, M. and R. Nowak, 2003, An EM algorithm for wavelet-basedimage restoration: IEEE Transactions on Image Processing, 12, 906–916.Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor nonunitary imaging operators: Geophysics, 69, 1017–1024.Hennenfent, G. and F. J. Herrmann, 2006a, Application of stable signalrecovery to seismic interpolation: Presented at the SEG InternationalExposition and 76th Annual Meeting.——–, 2006b, Seismic denoising with non-uniformly sampled curvelets:IEEE Computing in Science and Engineering, 8, 16–25.——–, 2007a, Irregular sampling: from aliasing to noise: Presented at theEAGE 69th Conference & Exhibition.——–, 2007b, Random sampling: new insights into the reconstruction ofcoarsely-sampled wavefields: Presented at the SEG International Exposi-tion and 77th Annual Meeting.Herrmann, F. J., 2007, Surface related multiple prediction from incompletedata: Presented at the EAGE 69th Conference & Exhibition.Herrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Nonlinear primary-multiple separation with directional curvelet frames: Geophysical Journal105Chapter 4. BibliographyInternational, 170, 781–799.Herrmann, F. J. and G. Hennenfent, 2007, Non-parametric seismic datarecovery with curvelet frames: Technical report, UBC Earth and OceanSciences Department. (TR-2007-1).Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2008, Sparsity- andcontinuity-promoting seismic imaging with curvelet frames: Applied andComputational Harmonic Analysis.Kuhl, H. and M. D. Sacchi, 2003, Least-squares wave-equation migration forAVP/AVA inversion: Geophysics, 68, 262–273.Lin, T. and F. J. Herrmann, 2007, Compressed wavefield extrapolation:Geophysics.Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration ofincomplete reflection data: Geophysics, 64, 208–221.O’Brien, M. and S. Gray, 1996, Can we image beneath salt?: The LeadingEdge, 15, 17–22.Sacchi, M. D., T. J. Ulrych, and C. Walker, 1998, Interpolation and extrap-olation using a high-resolution discrete Fourier transform: IEEE Trans-actions on Signal Processing, 46, 31–38.Santosa, F. and W. Symes, 1986, Linear inversion of band-limited reflectionseismogram: SIAM Journal on Scientific and Statistical Computing, 7,1307–1330.Taylor, H. L., S. C. Banks, and J. F. McCoy, 1979, Deconvolution with thelscript1 norm: Geophysics, 44, 39–52.Verschuur, D. J. and A. J. Berkhout, 1997, Estimation of multiple scatteringby iterative inversion, part II: practical aspects and examples: Geophysics,106Chapter 4. Bibliography62, 1596–1611.Ying, L., L. Demanet, and E. J. Cand´es, 2005, 3-D discrete curvelet trans-form: Wavelets XI, Expanded Abstracts, 591413, SPIE.107Chapter 5True amplitude depthmigration using curvelets5.1 IntroductionDespite the fact that modern depth migrations can accurately position re-flecting surfaces in the Earth subsurface, they often do not correctly recoverhigh-fidelity amplitude information. This is true even when accurate knowl-edge of the velocity model is available.Many approaches have been proposed to solve the “true-amplitude” mi-gration problem and they generally fall into three main categories. First,there are approaches that correct for the amplitude “during migration” byapplying imaging conditions that take reflector amplitudes into account (seee.g. J. C. Costa and Novais, 2009; Chattopadhyay and McMechan, 2008;J. Schleicher and Novais, 2008; Y. Zhang and Zhang, 2007).The second category involves linearized inversion where the linearizedscattering operator, i.e. the adjoint of migration (see e.g. Guitton, 2004;Claerbout, 1985; Gray, 1997; Symes, 2007), is inverted using iterative Lanc-A version of this chapter will be submitted for publication. Moghaddam, P.P., Her-rmann F. J. and Shahidi, R. (2010) True amplitude depth migration using curvelets.1085.1. Introductionzos methods (Herrmann et al. (2009)). In this category, a variety of algo-rithms for true amplitude migration have been introduced which iterativelyupdate the migrated image to minimize the misfit between observed and thesimulated data (see e.g. Tarantola, 1987; Nemeth et al., 1999; Kuehl andSacchi, 2003; Mulder and Plessix, 2004; Herrmann et al., 2009).Third, there are the so called image-to-image scaling methods, wherethe normal operator (i.e. migration followed by the scattering operator)is approximated through a transform-domain scaling method. The scalingis done by comparing the migrated image (with some possible pre/post-processings) with a second image obtained by remigration of the image (i.e.after applying the normal operator on the migrated image). Typically, thisdiagonal scaling is performed in some transform domain ( Fourier, curvelets,physical, etc, seee.g. Rickett,2003;PlessixandMulder,2004;Guitton,2004;Symes, 2008a; Herrmann et al., 2008).Our method falls in the third category, and we address the amplitude re-covery problem by exploiting the recently developed curvelet frames. Theseframe expansions compress seismic images (see e.g. Hennenfent and Her-rmann, 2006; E. J. Cand‘es and Ying, 2006) and consist of a collection offrameelements‘curvelets’thatareinvariantunderthe normaloperator(Her-rmann et al. (2008, 2009)). These two properties allow for the developmentof an approach where the normal operator is nonlinearly inverted, using theeigen-like behavior of curvelets, in combination with the ability of curveletsto handle conflicting dips and sparsely represent the seismic images.The main contributions of this chapter are three-fold. First, we replacethe linear least-squares formulation for the estimation of the curvelet-domain1095.2. Outline of the chaptercoefficients (Herrmann et al. (2008)) by a nonlinear least-squares formula-tion that approximates the symmetric positive definite normal operator byimposing positivity on the scaling coefficients. Second, we apply this methodto correct the amplitudes of images that suffer from migration noise causedby small number of shots, as well as noise in the data. Third, we comparetwo sparsity-promoting methods for recovery, a soft thresholding-correctiontechnique and a large scale sparsity solver.Our method is tested with reverse-time ‘wave equation’ migration codesimulating the acoustic wave equation on different synthetic models.5.2 Outline of the chapterIn the first section, we derive our formulation and propose two methods foramplitude recovery. Next, we derive the diagonal estimation for the normaloperator in the curvelet domain. The chapter ends with an example of eachmethod on the SEG salt model linearized Born dataset and comparison ofthe proposed methods followed by conclusion.5.3 Problem FormulationBy assuming our data consist of primaries only (i.e., the surface related andall the internal multiples are assumed to be removed before migration), thedata can be expressed as,d≈Km+n, (5.1)1105.3. Problem Formulationwith d data, K the linearized Born scattering operator (also known as thedemigration operator), m the perturbation in the velocity model (i.e., δvv3with v as velocity, see William W. Symes and Gao (2004); Symes (2008b) )and n zero-centered Gaussian noise with variance σ2n. Theoretically, migra-tion is the adjoint of the scattering operator. After applying the migrationto the data, we obtain the image,y = KTd (5.2)= Ψm+KTn (5.3)where y is the migrated data, KT is the migration operator and the Ψ is thenormal operator defined byΨ = KTK. In order to findmin above equation,we need to invert the normal operator (Ψ). To solve the above equation, wefollow recent results by Herrmann et al. (2008, 2009) and approximate thenormal operator by a curvelet-domain scaling, i.e.,we have,ΨmsimilarequalCTWCm, (5.4)with W a diagonal matrix. This diagonal approximation involves applyingthe curvelet transform C, followed by a scaling with the positive diagonalmatrix W and the inverse curvelet transform back to the physical domainby CT. [Note that the curvelets are a tight frame (E. J. Cand‘es and Ying,2006) so the adjoint of curvelet transform (CT) is equal to its pseudo-inverse1115.3. Problem Formulation(i.e., CTC = I)]. Inserting this approximation into Equation 5.2 leads to,ysimilarequalCTWCm+e, (5.5)with e = KTn colored noise with covariance,Cee = E(eeT) = E(KTnnTK) = KTE(nnT)K = σ2nKTK, (5.6)with E(.) the statistical expectation and Cnn = E(nnT) = σ2nI the covari-ance of the white Gaussian noise. To exploit curvelet-domain sparsity onthe recovered model, we rewrite Equation 5.5 as,y = CTWx+e (5.7)where x = Cm. In the context of stable signal recovery (see e.g. Cand`eset al., 2006; Daubechies et al., 2005), x can be found by solving the followingoptimization problem,P1 :˜x = argminx||x||1 = summationtextµ∈M|xµ|subject to||C−12ee (y−CTWx)||2≤υ˜m = CT˜x,(5.8)where υ is the noise dependent tolerance, ˜. stands for the estimated valueand µ is the index set of curvelets. The term C−12ee in front of misfit acts asthe whitening operator for the data misfit (Tarantola (1987)). During theoptimization, the vectorxis found by minimizing the lscript1 norm on the curveletrepresentation for the model and the lscript2 norm for data misfit. Equation 5.81125.3. Problem Formulationallows us to exploit curvelet-domain sparsity subject to fitting the migratedimage within a tolerance by solving the sparsity-promoting program.Inserting Equation 5.4, the approximation for the normal operator intoEquation 5.6 for the noise covariance yields,Cee = σ2nKTKsimilarequalσ2nCTWC. (5.9)Assuming the curvelet transform is nearly orthonormal, we can approximatethe action of C−12ee in Equation 5.8 by,C−12ee similarequalσ−1n CTW−12C. (5.10)With this approximation, Equation 5.8 becomes,P :˜x = argminx||x||1 subject to||b−CTW12x||2≤epsilon1˜m = CT˜x,(5.11)with b = CTW−12Cy and epsilon1 = σnυ. We propose two alternative methods forsolving above problem, one is a Basis Pursuit De-Noising problem (BPDN)and other one is a soft-thresholding solution. BPDN proposes a solution toproblem P equivalent to a weighted lscript1 problem (by changing the variableW12xmapsto→x) yielding ˜m = CTW−12˜x,BPDN : ˜x = argminx ||x||1,W−12 =summationdisplayi|xiw−1/2i |subject to||b−CTx||2≤epsilon1(5.12)We solve optimization BPDN with a projected gradient method (SPGL1),1135.4. Diagonal Approximation of the Normal Operatorrecently introduced by (van den Berg and Friedlander, 2008). An alternativeway can be evaluating P in unconstrained form. This is equivalent to aquadratic program yielding ˜m = CT˜x,QP : ˜x = argminx 12||b−CTx||22 +λ(epsilon1)||x||1,W−12(5.13)where λ(epsilon1) is a constant controlled by the noise level in the data. An ap-proximate solution for the optimization QP can be found with a weightedsoft-thresholding method, given by,˜m = CTsign(Cb)circledotmax(0,|(Cb)circledotw−12|−λ(epsilon1)), (5.14)where circledot is the element-wise product, w = diagW is a vectorized form ofthe diagonal elements for the matrix W, sign(·) is the element-wise signfunction and |·| the absolute value function. In the examples section, weexamine both methods in term of quality of their output.5.4 Diagonal Approximation of the NormalOperatorIn the previous section, we stated that we can diagonally approximate thenormal operator in the curvelet domain. This approximation (c.f., Equa-tion 5.4), replaces the construction of the normal operator, which is com-putationally expensive, with a curvelet domain diagonal scaling, which iscomputationally cheap. In this section, we provide the theoretical underpin-ning of such approximation by introducing the curvelets and their invariance1145.5. Curvelets and their invariance under normal operatorbound under the action of normal operator. We then propose a method toapproximate the normal operator with a diagonal scaling in curvelet domain.5.5 Curvelets and their invariance under normaloperatorA 2D curvelet φµ is defined by its index µ = (j,k,θ) with j = 0,1,2,... con-sisting of its scale subindex, θ = 0,1,...,2floorleftj/2floorright−1, its orientation subindex(floorleftxfloorright is the lower integer part of x ), and k = (kx,kz) the location in thewavenumber domain subindex which are scale and angle dependent. Fig-ure 5.1(a) shows an example of different curvelets with different scales andorientations and Figure 5.1(b) shows the corresponding Fourier spectrum.An important feature of curvelets is that the action of the normal oper-ator (i.e., Ψ = KTK ) on a curvelet φµ corresponds to its multiplication bya positive scalar. In Herrmann et al. (2008), we proved a theorem bound-ing the error with the normal operator diagonal approximation in curveletdomain, i.e., we showed,||Ψφµ−CTWCφµ||2≤β2−j/2, (5.15)with W a positive diagonal matrix as defined before and β a constant.According to this theorem, the approximation error decays with scale (i.e.,j), which is consistent with the high-frequency asymptotic behavior thatunderlies the diagonal approximation. Shahidi and Herrmann (2009a) showsthat the constant β is of moderate size.1155.5. Curvelets and their invariance under normal operatork10-0.25-0.50 0.25 0.50-0.50k2-0.2500.250.5001002003004005000 100 200 300 400 500SamplesSamplesFigure 5.1: Spatial and frequency representation of curvelets, (a) four dif-ferent curvelets in the spatial domain at three different scales, (b) dyadicpartitioning in the frequency domain, where each wedge corresponds to thefrequency support of a curvelet in the spatial domain. This figure illustratesthe microlocal correspondence between curvelets in the physical and Fourierdomain. Curvelets are characterized by rapid decay in the physical spaceand of compact support in the Fourier space. Notice the correspondencebetween the orientations of curvelets in the two domains. The 90o rotationis a property of the Fourier transform. Courtesy of Herrmann et al. (2008).1165.6. Diagonal approximationThe normal operator is pseudodifferential when no multipathing occurs(Beylkin, 1985; Rakesh, 1988; Ten-Kroode et al., 1998), and its order inspace dimension n is n−1 (Stolk, 2000). The inequality in Equation 5.15holds when the order of normal operator is zero (Herrmann et al., 2008).Ten-Kroode et al. (1998); Stolk (2000) show that in 2D applying an operatorof order−1/2 to both migration and modeling operators makes the operatorof zero order. As it is shown in Herrmann et al. (2009), this is equivalent tofractional integration of the source wavelet in both the modeling and migra-tion operators. Mathematically, this can be written as ∂−1/2t ·= F−∗|ω|−1/2·with F−∗ the inverse Fourier transform.Equation 5.15 states that curvelets are similar to eigenvectors of thenormal operator. Therefore, we propose a quasi-eigenvalue decompositionfor the normal operator with curvelets as quasi-eigenvectors. This typeof approach is similar to the Wavelet-Vaguellette Decomposition (WVD)proposed by (Donoho, 1995; Herrmann et al., 2008).5.6 Diagonal approximationIn Herrmann et al. (2008), we derive an approximation of the normal oper-ator by remigrating the migrated image again (i.e., mremig = KTKmmig).Subsequently, find a diagonal approximation for the normal operator usingthese image pairs. Mathematically,mremig = (KTK)mmig similarequalCTWCmmig (5.16)similarequal CTdiag(v)w, (5.17)1175.6. Diagonal approximationwith v = Cmmig and w representing the diagonal elements of W (i.e.,W = diag(w)).Equation5.16 aims to finda curvelet-domaindiagonalscaling vector(i.e.,w in Equation 5.16). Because the curvelet transform is redundant (i.e., ithas more rows than columns), inverting for w in Equation 5.16 needs solv-ing an underdetermined system of equations yielding non-unique solutions.To find a unique solution to this estimation problem, we exploit additionalproperties of the normal operator. First, the normal operator is positive def-inite, yielding positive entries for the scaling vector w. Second, the normaloperator is elliptic Pseudo-differential (Stolk and De-Hoop, 2006; Shahidiand Herrmann, 2009a) and is invertible when the symbols of this operatorare smooth in both physical and Fourier space, when the background veloc-ity model is smooth (Stolk and De-Hoop, 2006; Herrmann et al., 2008). Weuse these two properties to formulate the diagonal estimation in terms of anonlinear least-squares problem, which yields positivity of the entries of thediagonal and smoothness amongst neighboring curvelet coefficients in spaceand angle domain. Compared to our earlier work (Herrmann et al., 2008),the positivity assumption is new in the context of imaging (see related workon application of this method for multiple removal, Herrmann et al. (2007)) and the diagonal estimated by following,Pd :˜z = argminzJ(z) = 12||mmig−CTdiag(v)ez||2 +κ||Lz||2,˜w = (diag)(˜w) with ˜w = e˜z,(5.18)where w is replaced by ez to ensure positivity of the estimated diagonal1185.7. Practical Workflowelements (W). In above equation, phase-space smoothness is promoted byminimizing lscript2 norm of the coefficients after applying the sharpening oper-ator, L = [Dx Dz Dθ]. This operator penalizes fluctuations amongstneighboring coefficients in ez, i.e., coefficients that are close in phase spaceboth in position and angle (for more details on this topic see Shahidi andHerrmann (2009b)). The matrices Dx,z contain first-order differences ateach scale in the x,z directions, and Dθ the first-order difference in the θdirection. Because the curvelet grid differs for each scale and angle, thissharpening operator is implemented for each curvelet wedge separately. κ isa smoothness parameter that control the smoothness of the scaling versusdata misfit. In this paper, we select this parameter by exhaustive search forautomatic parameter selection (Shahidi and Herrmann (2009a)).Equation 5.18 is convex and can be minimized by an efficient numericaloptimization method. We use Limited-memory BFGS (see e.g. Nocedal andWright, 1999), using following gradient,∇zJ(z) = ezcircledotdiag(v)C(CTdiag(v)ez−mmig)+2κLTLz. (5.19)5.7 Practical WorkflowOur algorithm consists of two main steps, namely the calculation of thecurvelet scaling coefficients via an image-to-remigrated image matched fil-tering and the subsequent estimation of amplitude-corrected, denoised imageby sparsity promotion.1195.7. Practical WorkflowIn order to apply the curvelet-domain match filter, we require a curveletimplementation (E. Candes (2005)). Our algorithm also requires access tothe Born modeling operator and its adjoint, the migration operator. For thispurpose, we use Symes (2007) reverse-time migration with optimal check-pointing and its linearized modeling.We make the normal operator zero order by applying a zero-phase half-time integration to the source wavelet input to both operators (Ten-Kroodeet al., 1998; Symes, 2008a). Our algorithm consists of the following steps:1. Perform migration and depth correction–i.e., d mapsto→ DZKTd = mmig,with mmig depth corrected migrated image and DZ = diag(z) wherezi = itrianglez, for i = 1 : nz, trianglez is the depth grid spacing and nz is thenumber of samples in depth (Herrmann et al. (2009)). The reasonfor depth correction is to amplify the amplitudes of deep events in themigrated image. This ensures that after applying the normal operator,the amplitude of the deep events are partially calibrated so it can beused to extract information from the normal operator via estimationof its diagonal in curvelet domain.2. Remigrateandresimulatethemigratedimage–i.e., KTK˜mmig = mremig.3. FindthecurveletdomaindiagonalscalingWforwhich ˜mmig similarequalCTWCmremigby solving the optimization problem Pd (cf., Equation 5.18).4. Construct the original image by solving the minimization problem QPwith weighted soft-thresholding (Equation 5.14) or the optimizationproblem BPDN using a solver for large-scale sparse reconstruction.1205.8. Examples5.8 ExamplesIn this section, we apply our amplitude recovery methods on the (noisy)linearized Born dataset with noise using SEG-AA salt model (Aminzadehet al. (1997); O’Brien and Gray (1996)). The SEGAA model shown inFigure 5.2(a) is 5.28km in depth and 16.8km in offset with grid resolutionof 12m. The 5000 time-step linearized Born modeling takes about 24 minwith 25 CPUs on a MPI cluster, while the migration takes about 90 min.In this section, we compare the soft-thresholding and the BPDN methodsdescribed in this paper in terms of output quality and calculation speedfor both no noise in data (migration noise only) and when noise is addedto data. We also analyze the accuracy of the estimation by comparing therecovered image traces with the reflectivity.5.8.1 Noise-free caseSince the removal of the migration noise is our primary interest, we conductan imaging experiment with a small number of shots. This dimensionalityreduction results in cheaper acquisition and faster computation but goes atthe expense of introducing the migration noise. This experiment serves asan illustration on how well our method performs in both the removal ofmigration noise and restoration of the amplitudes.Weusealand-acquisitionscenariowith50shotsbetween50mand16.5kmat a 336m spacing and 680 receivers per shot at a 24m interval with offsetsbetween 200m and 16.5km. To generate the linearized Born dataset, wesmooth this velocity model as shown in Figure 5.2(b) and use the difference1215.8. Examples(a)(b)(c)Figure 5.2: Conflicting dips velocity model and reflectivity (a) velocitymodel, (b) smooth velocity model used to generate our example, (c) re-flectivity generated by subtracting smooth velocity model from original one1225.8. Examplesbetween the smooth and the original velocity as the reflectivity shown inFigure 5.2(c).Figure 5.3 summarizes our approximation of the normal operator withthe curvelet-domain scaling. We chose the reflectivity in Figure 5.2(c) andapplied the scattering operator to generate linearized Born dataset (i.e.,d = Km). We used the smoothed velocity model shown in Figure 5.2(b)as background velocity for the scattering operator. The depth correctedmigrated result is included in Figure 5.3(a). After another modeling andmigration of this image, remigrated image is shown in Figure 5.3(b). Thesetwo images serve as input to our curvelet scaling coefficient estimation. Ap-plying this curvelet-domain scaling to the reference image yields a good ap-proximation to the action of normal operator as can be seen in Figure 5.3(c).The relative approximation error in this case is 3.71%.To illustrate the importance of imposing curvelet domain smoothnessduring the estimation of diagonal scaling, we compare the results of twoexperiments, one without κ = 0 in Equation 5.19, and one with a smoothingregularization (κ = 0.01 found empirically and used to generate the examplein Figure 5.3. Figure 5.4 shows the plot of coefficients for each wedge of thereference vector (Figure 5.4(a)), the diagonal estimation without smoothness(Figure 5.4(b)) and with the smoothness (Figure 5.4(c)) constraints. Thesefigures show curvelet decomposition of a shot gather at different frequencyband (scale) and angles (dip). Five scales are used for the decomposition.The centre (coarsest scale) shows the DC and low frequency components ofthe shot gather. The second coarsest scale has 8 angles. The number ofangles is doubled to 16 at the third and fourth scale and doubled again to1235.8. Examples(a)(b)(c)Figure 5.3: Diagonal approximation of normal operator, (a) reference image(depth corrected migrated image), (b) normal operator on the referenceimage (i.e., Ψr) (c) approximate normal operator on the reference image(i.e., CTWCr), approximation error is %3.711245.8. Examples32 for the fifth scale. Note the portions of shot gather captured at variousangles and scales. As one can see, the curvelet coefficients for the non zerosmoothing regularization κ = .01 (Figure 5.4(c)) is smooth over differentangles and locations.Soft-thresholding solutionWe investigate the effect of the threshold λ on the soft-thresholding solutionfor the optimization problem QP. Figure 5.5 shows the result of solving theoptimization problem QP with soft thresholding for different values of thethreshold λ. The solution of soft thresholding shown in Figure 5.5(a)is forwhen λ = 0 . This solution is similar to approximately inverting the normaloperator in the curvelet domain and applying it to the migrated image (i.e.,in Equation 5.11 ˜m = CTW−12Cb). Figure 5.5(b) shows the result forλ = σb and Figure 5.5(c) for λ = 2σb). We clearly see from the figures thatby increasing λ the incoherent noise removed from the image at the expenseof dimming the amplitude of the bottom reflectors.BPDN solutionHere we investigate the effect of the tolerance epsilon1 on the solution for theoptimizationproblemBPDNfornoise-freedata. Figure5.6showthe resultsfor solving the optimization problem BPDN using SPGL1 solver (Berg andFriedlander, 2007; van den Berg and Friedlander, 2008) for different valueof the tolerance epsilon1. The solution of SPGL1 shown in Figure 5.6(a) is forwhen epsilon1 is small (epsilon1 = σb with σb is the standard deviation of b in BPDN).This solution is similar to Figure 5.5(a) and corresponds to applying the1255.8. ExamplesscalesAngles(a)scalesAngles(b)scalesAngles(c)Figure 5.4: Curvelet-domain representation of the curvelet vector obtainedfrom (a) reference vector, (b) scaling coefficients without smoothing con-straints, (c), scaling coefficients with smoothing constraint. The differentsub-images represent the curvelet coefficients at different scales (coarsest inthe center) and different angles.1265.8. Examples(a)(b)(c)Figure 5.5: Examples of the solution for optimization problem QP withsoft-thresholding, (a) λ = 0, (b) λ = σb (c) λ = 2σb1275.8. Examplesinverse of the diagonal approximation to b. Figure 5.6(b) shows the resultfor larger epsilon1 (epsilon1 = 50σb) and Figure 5.6(c) is the solution when the epsilon1 is veryhigh (epsilon1 = 100σb). As can be seen from the figures, when epsilon1 is small, the imagecontains many artifacts that increases with depth and when the epsilon1 is large,these artifacts removed at the expense of dimming the bottom reflectors.The best choice of epsilon1 is shown in Figure 5.6(b) which is the result of testinga range of different value in the solver. In this image, the deep artifacts arevisibly removed and the bottom reflectors are enhanced.Amplitude AnalysisTo investigate the performance and accuracy of our amplitude recoverymethod, we examine a number of traces in the recovered images and com-pare them to the original reflectivity. Figure 5.7 shows the location of thetraces. To examine the performance of our method in the presence of a saltbody, we choose a number of traces inside and outside the salt body. Theoffset location of traces are (3480, 5160, 5880, 8640, 10080, and 12720) m,respectively. Figure 5.8 shows the amplitude of the traces along the off-set lines (3480, 5160, and 5880) m, respectively. Except small amount ofresidual noise in the recovered images, all the strong phases well recovered.Early arrivals are generally well recovered; however, with increasing depththe amount of noise increases.Figure 5.9 shows the amplitude of the traces along the offset lines of(8640, 10080, and 12720) m, respectively. Note that these traces passthrough the salt body. Similar to Figure 5.8, all important phases in thereflectivity are well matched in the recovered image, except the phases that1285.8. Examples(a)(b)(c)Figure 5.6: Examples of the solution for optimization problem BPDN withSPGL1, (a) epsilon1 = σb, (b) epsilon1 = 50σb (c) epsilon1 = 100σb1295.8. Examples(a)Figure 5.7: Reflectivity model, the vertical lines are locations of investigationcome from the bottom salt. Salt bottom imaging is one of the most chal-lenging fields in the seismic imaging. One reason might be the difference inthe nature of the salt bottom reflector. This means the amplitude of thereflector should be interpreted as change in the acoustic impedance in thatlocation rather than the simply perturbation in the velocity model.5.8.2 Noisy dataWe examine the effectiveness of our recovery method for noisy data. Forthis purpose, we add white Gaussian noise to the linearized Born data withSNR = 10 dB and re-run our method on the data. Figure 5.10 shows databefore and after adding the noise. Figure 5.11 shows the depth correctedmigrated image for this dataset. Note that the noise is obscuring the deepreflectors in the migrated image.1305.8. Examples0 500 1000150020002500300035004000−10−8−6−4−202468Depth (m)Amplitudeoffset at 2880 (m)SPGL1Soft ThresholdingReflectivity(a)0 500 1000150020002500300035004000−10−8−6−4−202468Depth (m)Amplitudeoffset at 4560 (m)SPGL1Soft ThresholdingReflectivity(b)0 500 1000150020002500300035004000−15−10−50510Depth (m)Amplitudeoffset at 5280 (m)SPGL1Soft ThresholdingReflectivity(c)Figure 5.8: Amplitude analysis of the recovered images using proposedmeth-ods, (a) offset=3480 (m) , (b) offset=5160 (m), (c) offset=5880 (m)1315.8. Examples0 500 1000150020002500300035004000−12−10−8−6−4−202468Depth (m)Amplitudeoffset at 8040 (m)SPGL1Soft ThresholdingReflectivitySalt−bottom (a)0 500 1000150020002500300035004000−15−10−5051015Depth (m)Amplitudeoffset at 9480 (m)SPGL1Soft ThresholdingReflectivitySalt−bottom (b)0 500 1000150020002500300035004000−14−12−10−8−6−4−20246Depth (m)Amplitudeoffset at 12120 (m)SPGL1Soft ThresholdingReflectivitySalt−bottom(c)Figure 5.9: Amplitude analysis of the recovered images using proposedmeth-ods, (a)offset= 8640 (m), (b)offset=10080 (m) , (c)offset=12720 (m), notethe salt bottom amplitude mismatch1325.8. Examples(a)(b)Figure 5.10: Synthetic data, (a) before , (b) after adding the noise, SNR =10dB1335.8. Examples(a)Figure 5.11: Depth corrected migrated noisy imageSoft-thresholding solutionWe investigate the effect of the threshold λ on the soft-thresholding solutionfor the optimization problem QP. Figure 5.12 shows the result of solvingthe optimization problem QP with soft thresholding for different values ofthe threshold λ. The solution of soft thresholding shown in Figure 5.12(a)isfor when λ = σb . Figure 5.12(b) shows the result for λ = 2σb and Fig-ure 5.12(c) is for when λ = 3σb). As can be seen from the figures, byincreasing λ the incoherent noise is visibly removed from the image at theexpense of dimming the amplitude of the bottom reflectors. The best choiseof λ is in Figure 5.12(b).BPDN solutionFigure 5.13 shows the result of solving the optimization problem BPDN us-ingSPGL1solverfordifferentvalueofthe toleranceepsilon1. The solutionofBPDN1345.8. Examples(a)(b)(c)Figure 5.12: Examples of the solution for optimization problem QP withsoft-thresholding, (a) λ = σb, (b) λ = 2σb (c) λ = 3σb1355.9. Tolerance analysisis shown in Figure 5.13(a) for small epsilon1 (epsilon1 = 10σb with σb is the standard de-viation of b in BPDN). This solution is similar to approximately invertingthe approximation of normal operator in curvelet domain and applying it tomigrated image (i.e., in Equation 5.11 ˜m = CTW−12Cb). Figure 5.13(b)shows the result for when epsilon1 is larger (epsilon1 = 140σb) and Figure 5.13(c) is thesolution when the epsilon1 is high (epsilon1 = 200σb). We can clearly see from the figuresthat when epsilon1 is small, the image contains high amount of noise and when theepsilon1 is large, the noise removed with expense of dimming the bottom reflectors.The best choice of epsilon1 in Figure 5.13(b) is found through exhaustive processof running the solver on a wide range of epsilon1 and compare the results.5.9 Tolerance analysisIn the last section, we did not specifically analyze the value of the tolerance epsilon1in the optimization problem BPDN rather we empirically ran the problemwith different values of epsilon1 and compared the results. A candidate value for thetolerance can be the statistical expectation of the misfit function in BPDNin Equation 5.8,˜υ = E[||C−12ee (y−CTWx)||2] (5.20)with some manipulation ˜υ = N with N is the number of elements in yin lexicographical order. Indeed Figure 5.13(b) is obtained for υ = N orepsilon1 = Nσn.1365.9. Tolerance analysis(a)(b)(c)Figure 5.13: Examples of the solution for optimization problem BPDNwith SPGL1, (a) epsilon1 = 10σb, (b) epsilon1 = 140σb (c) epsilon1 = 200σb1375.10. Verification of approximation5.10 Verification of approximationThroughout our formulation, we made a subtle approximation from Equa-tion 5.8 to Equation 5.11,CTW−12CCTWxsimilarequalCTW12x. (5.21)In this approximation, we assume the projection CCT is approximatelyequal to identity. Here we investigate the accuracy of this approximation.We assume x = Cm with m is the true reflectivity. We call the lefthand of above expression as El = CTW−12CCTWx and the right hand asEr = CTW12x.Figure 5.14 show the comparison between the terms in Equation 5.21.Figure 5.14(a) is the left hand El expression and Figure 5.14(b) is the righthand Er expression. The relative error between two images is 4.09 percent.5.11 ComparisonComparing the output quality and the performance, we suggest using thesoft-thresholding technique for amplitude recovery. Soft thresholding inmost cases results in better or similar quality of the output with less com-putational burden comparing to BPDN.1385.11. Comparison(a)(b)Figure 5.14: Comparison between the approximation terms in Equation5.21, (a) left hand term El = CTW−12CCTWx , (b) right hand termEr = CTW12x, relative error is 4.09 percent.1395.12. Conclusion5.12 ConclusionThe method presented in this chapter, is a revisit of the earlier work (Her-rmann et al. (2008)), by taking into account the migration noise and the datanoise. The method relies on the curvelet transform in both approximationand inversion of the normal operator and has following features,• speeding up the calculation of the normal operator by diagonalizing itin curvelet domain.• efficiently using the compression of seismic images by curvelets to solvea curvelet-sparsity optimization problem.• bringing the amplitude correction problem within the context of stablesignal recovery.The results of applying our method on synthetic data suggest that ourmigration amplitude recovery can be both efficient and effective in eliminat-ing migration artifacts and recovering the amplitudes in the seismic image.The recovered images showed partial elimination of noise, improved spatialresolution, and enhanced reflectivity amplitude.5.13 AcknowledgmentsThe authors would like to thank the authors of CurveLab for making theircodes available. The authors would also like to thank W. Symes for pro-viding the migration code and TOTAL E&P for providing the synthetic BPdata. This work was in part financially supported by the Natural Sciences1405.13. Acknowledgmentsand Engineering Research Council of Canada Discovery Grant (22R81254)and Collaborative Research and Development Grant DNOISE (334810-05)of Felix J. Herrmann and was carried out as part of the SINBAD projectwith support, secured through ITF (the Industry Technology Facilitator),from the following organizations: BG Group, BP, Chevron, ExxonMobil andShell.141BibliographyAminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D Salt and Overthrust Model.SEG/EAGE 3-D Modeling Series, No. 1: Society of Exploration Geophysi-cists, Tulsa.Berg, E. V. D. and M. Friedlander, 2007, In pursuit of a root: OptimizationOnline:2007.06.1708, 8.Beylkin, G., 1985, Imaging of the discontinuities in the inverse scatteringproblem by inversion of casual generalized curvelet transform: Journal ofMathematical Physics, 26, 99–108.Cand`es, E. J., J. Romberg, and T. Tao, 2006, Stable signal recovery fromincomplete and inaccurate measurements: 59, 1207–1223.Chattopadhyay, S. and G. A. McMechan, 2008, Imaging conditions forprestack reverse-time migration: Geophysics, 73, S81–S89.Claerbout, J., 1985, Earth soundings analysis, processing versus inversion:Blackwell Scientific Publications.Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterative thresholdingalgorithm for linear inverse problems with a sparsity constraints: 57,1413–1457.Donoho, D., 1995, Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition: Applied and Computational Harmonic Analy-142Chapter 5. Bibliographysis 2, 101 – 126.E. Candes, L. Demanet, D. D. L. Y., 2005, Fast discrete curvelet transforms.E. J. Cand‘es, L. Demanet, D. L. D. and L. Ying, 2006, Fast discrete curvelettransforms: SIAM Multiscale Model. Simul., 5, 861899.Gray, S., 1997, True amplitude seismic migration: A comparison of threeapproaches: Geophysics, 62, 929–936.Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor nonunitary imaging operators: Geophysics, 69, 1017–1024.Hennenfent, G. and F. J. Herrmann, 2006, Seismic denoising with non-uniformly sampled curvelets: Comp. in Sci. and Eng., 8, 16–25.Herrmann, F., C. Brown, Y. Erlangga, and P. Moghaddam, 2009, Curvelet-based migration preconditioning and scaling: Geophysics, 74, A41–A45.Herrmann, F., P. Moghaddam, and C. Stolk, 2008, Sparsity- and continuity-promotingseismicimagerecoverywithcurveletframes: AppliedandCom-putational Harmonic Analysis 24, 2, 150–173.Herrmann, F. J., D. Wang, and D. J. Verschuur, 2007, Adaptive curvelet-domain primary-multiple separation. ((Submitted to Geophysics)).J. C. Costa, F. A. Silva Neto, M. R. M. A. J. S. and A. Novais, 2009,Obliquity-correction imaging condition for reverse time migration: Geo-physics, 74, S57–S66.J. Schleicher, J. C. C. and A. Novais, 2008, A comparison of imaging condi-tions for wave-equation shot-prole migration: Geophysics, 73, S219–S227.Kuehl, H. and M. Sacchi, 2003, Least-squares wave-equation migration forAVP/AVA inversion: Geophysics, 68, 262–273.Mulder, W. and R.-E. Plessix, 2004, A comparison between one-way and143Chapter 5. Bibliographytwo-way wave equation migration: Geophysics, 69, 14911504.Nemeth, T., C. Wu, and G. Schuster, 1999, Least-squares migration of in-complete reflection data: Geophysic, 64, 208–221.Nocedal, J. and W. Wright, 1999, Numerical optimization: Springer Verlag.O’Brien, M. and S. Gray, 1996, Can we image beneath salt?: The LeadingEdge, 15, 17–22.Plessix, R. and W. Mulder, 2004, Frequency-domain finite-differenceamplitude-preserving migration: Geophys. J. Int., 157, 975–987.Rakesh, A., 1988, A linearized inverse problem for the wave equation:Comm. on P.D.E., 13, 573–601.Rickett, J., 2003, Illumination based normalization for wave-equation depthmigration: Geophysics, 68, 208–221.Shahidi, R. and F. Herrmann, 2009a, Curvelet-domain matched filtering:under preparation.——–, 2009b, Curvelet-domain matched filtering with frequency-domainregularization: 28.Stolk, C., 2000, Microlocal analysis of a seismic linearized inverse problem:Wave Motion, 32, 267–290.Stolk, C. and M. De-Hoop, 2006, Seismic inverse scattering in the downwardcontinuation approach: Wave Motion, 43, 579–598.Symes, W. W., 2007, Reverse time migration with optimal checkpointing:Geophysics, 72, SM213–SM221.——–, 2008a, Approximate linearized inversion by optimal scaling ofprestack depth migration: Geophysics, 73, R23–R35.——–, 2008b, Migration velocity analysis and waveform inversion: Geophys-144Chapter 5. Bibliographyical Prospecting, 56, 765 – 790.Tarantola, A., 1987, Inverse problem theory: Elseweir.Ten-Kroode, A., D. Smith, and A. Verdel, 1998, A microlocal analysis ofmigration: Wave Motion, 28.van den Berg, E. and M. P. Friedlander, 2008, Probing the pareto frontierfor basis pursuit solutions: SIAM Journal on Scientific Computing, 31,890–912.William W. Symes, Christiaan C. Stolk, B. B. and F. Gao, 2004, Reversetime shot-geophone migration: Technical Report.Y. Zhang, S. Xu, N. B. and G. Zhang, 2007, True-amplitude, angle-domain,common-image gathers from one-way wave-equation migrations: Geo-physics, 72, S49–S58.145Chapter 6ConclusionsIn this chapter I summarize the main contributions of this thesis anddiscuss some limitations of the work presented. I also suggest follow-upwork as well as possible extensions.6.1 Main contributionsThe topic of this thesis is seismic image amplitude correction which isdirectly applied to seismic imaging. The approach I advocate is to view seis-mic imaging in a mathematical perspective. I identify a transform, calledthe curvelet transform (Cand`es and Donoho, 2004), to that effect and useit in a new formulation of the seismic images amplitude correction problem.The most challenging operator in the seismic imaging problem, called nor-mal operator, is a computationally exhausting. This operator needs to beinverted during the imaging process to give an accurate understanding aboutthe reflectors amplitude in a seismic image. Each element of the curvelettransform, namely curvelet, has an attractive property under the normaloperator. In this thesis, I proved that each curvelet remains approximatelyinvariant under the normal operator, which gives us an exceptional abilityto approximately invert the normal operator with only one iteration. I call1466.1. Main contributionsthis approach the curvelet match filtering (CMF).I leverage this approximation towards the development of a stable am-plitude recovery by coining it with sparsity-promoting inversion, which isbasically a large-scale one-norm solver.The remainder of this section provides more details about the contribu-tions.6.1.1 Sparsity and continuity promoting seismic imagerecovery with curvelet framesI combine the compression of images by curvelets with the invarianceof this transform under the normal operator. This combination allows fora formulation of a stable recovery method for seismic amplitudes. Duringthe recovery, the normal operator is approximately inverted. Compared toother approaches for migration scaling, the presented method (i) includes atheoretical bound on the L2−error for the diagonal approximation in thecurvelet domain; (ii) prescribes a procedure for the estimation of the diag-onal from numerical implementations of the imaging operators; (iii) formu-lates the amplitude recovery problem as a nonlinear optimization problem,where the inversion of the diagonalized normal operator is regularized byimposing sparsity in the curvelet domain and continuity along the imagedreflectors.The diagonal approximation is used to formulate the seismic amplituderecovery in terms of a constrained optimization problem. The amplitude-correctedimageisobtainedbysolvingthissparsity-andcontinuity-promotingoptimization problem. The invariance of curvelets under the normal opera-1476.1. Main contributionstor preserves the sparsity. The cost of computing the diagonal approxima-tion is one demigration-migration per reference vector, which is much lesscompared to the cost of Krylov-based least-squares inversion. The recov-ery results show an overall improvement of the image quality. The joinedsparsity- and continuity-enhanced image has diminished artifacts, improvedresolution and recovered amplitudes.6.1.2 Curvelet-based migration preconditioning and scalingI present a method to approximate the normal operator and use this ap-proximation as a preconditioner in the least squares framework. Becauseof the size of the seismic imaging problem, preconditioning of least-squaresmigration is an interesting topic. Traditionally, the first few iterations of theLSQR algorithm of least-squares migration are known to make significantprogress towards the solution. Unfortunately, even for this limited numberof iterations, the computational costs are often still prohibitively large forpractical problems. Our method, partly resolves this issue through a com-bination of left and right preconditioning together with a curvelet-domainscaling. Inclusion of the latter proved particularly important because it re-stores the amplitudes and leads to faster convergence, at a relatively smallcomputational overhead. Preconditioning also plays a pivotal role in mak-ing this approach numerically feasible and may extend to a solution of thefull-waveform inversion problem with sparsity promotion. Both approachesare justified by ample evidence that curvelets are sparse on the model.1486.1. Main contributions6.1.3 Curvelet-based seismic data processingBeside the seismic wavefield reconstruction problem, I recast a few otherprocessingsteps—signalseparation, migration-amplituderecovery, anddeconvolution—in a sparsity-promoting program that exploits the high degree of sparsityattained by curvelets on seismic data and images. The promising resultsobtained show that the insights gained from the developments of CMF canbe leveraged towards a much broader range of applications. This prospectopens an exciting new outlook towards future developments in explorationseismology.6.1.4 True amplitude depth migration using curveletsI presented a fast and robust approach for approximation of the normaloperator by revisiting of the earlier work and taking into account the migra-tion noise and the data noise. Compared to other approaches for migrationamplitude recovery, some improvements in the design of curvelet amplituderecovery method including, speeding up the calculation of the normal oper-ator by diagonalizing it in curvelet domain, designing an efficient approxi-mation that takes into account a laterally-variant velocity models and steepreflectors. The results of applying method on synthetic data suggest thatmigration amplitude recovery method can be both efficient and effective ineliminating migration artifacts and recovering the amplitudes in the seismicimage. The recovered images showed partial elimination of noise, improvedspatial resolution, and enhanced reflectivity amplitude.1496.2. Follow-up work6.2 Follow-up workI suggest a few ideas that go beyond the reported experiments.6.2.1 Full waveform inversionOur estimation for normal operator can be directly applied to the contextof full waveform inversion in which the normal operator is inverted and usediteratively. In this context, the velocity model (not just reflectivity in themigration context) is estimated directly from data (see e.g. Symes, 2008).6.2.2 3D true amplitude migrationOur 2D true amplitude migration can be easily extend to 3D. This canbe done in two ways, first, I apply the curvelet smoothing in 2D and for 3rddimensionI canexploit wavelet smoothing. Second, I canuse the 3Dcurvelettransform and extend the method to 3D curvelet diagonal estimation. Thelatter, needs a sophisticated algorithm to locate the closest curvelets in 3Drather than 2D for smoothing operator.6.3 Adding more priori knowledge to solutionIn this thesis, I propose some ideas for amplitude recovery which pre-sume that the seismic image are sparse in the curvelet domain. The commontheme of most ideas is that there is no priori knowledge of the geologicalstructure of the survey area. However, if a priori knowledge about the struc-ture already obtained through the well-log or other geological/geophysical1506.4. Current limitationsexperiments, it can be used within stable recovery framework in addition ofsparsity to reconstruct seismic wavefields.6.4 Current limitationsI examine both the practical and the fundamental weaknesses of thecurrent curvelet transform.6.4.1 Curvelet codeThe true amplitude migration results presented in this thesis were ob-tained using the FDCT based on the wrapping of specially selected Fouriersamples (Cand`es et al., 2005). This implementation breaks down the inputimage or volume into a number of scales depending on the length of theshortest axis. In other words, if one axis is much shorter than the others,the decomposition along the long axes is unnecessarily limited. Despite anincreased implementation complexity, an alternative would be to extend theshort axis to be sufficiently large enough. This extension can be done bypadding the image along the short axis or wrapping it.A more fundamental limitation of the FDCT is related to the redun-dancy of the transform. Indeed, the FDCT is around 8-redundant in 2Dand around 24-redundant in 3D, which precludes, at least for now, tractablehigher-dimensional FDCTs. Lu and Do (2007) propose a less redundant N-dimensional (N ≥2) implementation, termed surfacelet transform, by com-bining a directional filter bank with a multiscale pyramid. Another option isto combine the curvelet transform with another transform (Herrmann, 2003;1516.4. Current limitationsNeelamani et al., 2008) to reduce redundancy and reach higher dimensions.The different treatment of the axes is unsatisfactory in several applications(see, e.g., Neelamani et al., 2008), though. For interest, Kutyniok and La-bate (2005) propose yet another N-dimensional (N ≥2) transform, calledshearlet transform, but no discrete implementation is available at this pointto determine the redundancy and the effectiveness of shearlets for wavefieldreconstruction. In addition, in some cases, defining the smoothing opera-tor in the transform domain is not as straightforward as curvelets since insome cases finding the closest elements in transform domain might requirea significant amount of operations.152BibliographyCand`es, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005, Curvelab:software. (Available at http://www.curvelet.org).Cand`es, E. J. and D. L. Donoho, 2004, New tight frames of curvelets andoptimal representations of objects with C2 singularities: Communicationson Pure and Applied Mathematics, 57, no. 2, 219 – 266.Herrmann, F. J., 2003, Multifractional splines: Application to seismic imag-ing.Kutyniok, G. and D. Labate, 2005, Shearlets: website. (www.shearlet.org).Lu, Y. M. and M. N. Do, 2007, Multidimensional directional filter banksand surfacelets: Transactions on Image Processing, 16, no. 4.Neelamani, R., A. I. Baumstein, D. G. Gillard, M. T. Hadid, and W. L.Soroka, 2008, Coherent and random noise attenuation using curvelettransforms: The Leading Edge. (In press).Symes, W. W., 2008, Migration velocity analysis and waveform inversion:Geophysical Prospecting, 56, no. 6, 765–790.153Appendix AProofs of lemma 1 andtheorem 1We first prove Theorem 1 from Lemma 1, we then prove the Lemma whichappears in chapter 2.Proof: [Proof of Theorem 1] CT acting on a vector c ={cµ}µ∈M in lscript2 issimply given by CTc = summationtextν cνϕν. Hence,CTDΨCϕµ =summationdisplayν(Cϕµ)νa(xν,ξν)ϕν.The difference (Ψ(x,D)−CTDΨC)ϕµ can be written as, using that CTC =Id,summationdisplayν(Cϕµ)ν(Ψ(x,D)−a(xν,ξν))ϕνThe vector (Cϕµ)ν is bounded in lscript2(M) of vectors with curvelet indices, andits entries can only be nonzero if|µ|−2≤|ν|≤|µ|+ 2, since the supportof two curvelets in the Fourier domain is disjoint when curvelets are two orA version of this appendix has been published. Herrman, F.J., Moghaddam, P.P. andStolk, C. (2008) Sparsity- and continuity-promoting seismic image recovery with curveletframes. Applied and Computational Harmonic Analysis, Vol. 24, No. 2, pp. 150-173,2008.154Appendix A. Proofs of lemma 1 and theorem 1more scales apart. It follows thatbardbl(Ψ(x,D)−CTDΨC)ϕµbardbl2L2(Rn)≤C6summationdisplayνbardbl(Cϕµ)ν(A(x,D)−a(xν,ξν))ϕνbardbl2L2(Rd)≤C7summationdisplayν2−|ν||(Cϕµ)ν|2≤C82−|µ|/2.Proof: [Proof of Lemma 1] Possibly after a coordinate transformation,which does not affect the pseudodifferential nature of Ψ, we may assumethat ξµ = (0,...,0,λ), we will write ξprime = (ξ1,...,ξd−1).We may take|ν|greater than some minimum value, so that the supportin the wavenumber domain is away from ξ = 0. We can then assume that theprincipal symbol is homogeneous of order 0 in the support of the curvelet.On the support of the curvelet in the Fourier domain the symbol can bewritten as a(x, ξprimeξn). We apply a preparation theorem to this term. Thereare C∞ functions bn(x,ξprime/ξd), n = 1,..., 2d−1, such thata(x,ξprime/ξd)−a(xν,ξprimeν/ξν,d)=dsummationdisplayn=1(x−xν)nbn(x,ξprime/ξd)+d−1summationdisplayn=1parenleftbiggξnξd −ξν,nξν,dparenrightbiggbd+n(x,ξprime/ξd).Therefore, we can writeparenleftbigΨ(x,Dprime/Dd)−a(xν,ξprimeν/ξν,d)parenrightbigϕν=bracketleftBigg dsummationdisplayn=1(x−xν)nbn(x,Dprime/Dd)+d−1summationdisplayn=1bd+n(x,Dprime/Dd)DnDdbracketrightBiggϕν. (A.1)Consider first the contribution for one of the bn with 1 ≤ n ≤ d. The155Appendix A. Proofs of lemma 1 and theorem 1operator acting on the curvelet reads(x−xν)nbn(x,Dprime/Dd)ϕν= bn(x,Dprime/Dd)(x−xν)nϕν +[(x−xν)n,bn(x,Dprime/Dd)]ϕν.Because of the support and decay properties of curvelets we have bardbl(x−xν)nϕνbardblL2 ≤C−|ν|/29 , thereforebardblbn(x,Dprime/Dd)(x−xν)nϕνbardblL2(Rd)≤C102−|ν|/2Bythecalculusofpseudodifferentialoperators, theoperator[(x−xν)n,bn(x,Dprime/Dd)]is a pseudodifferential operator of order −1, which is continuous H−1(Rd)to L2(Rd) thereforebardbl[(x−xν)n,bn(x,Dprime/Dd)]ϕνbardblL2(Rd)≤C10bardblϕνbardblH−1(Rd)≤C112−|ν|.Adding the previous two estimates, we find thatbardbl(x−xν)nbn(x,Dprime/Dd)ϕνbardblL2(Rd)≤C122−|ν|/2. (A.2)Wenowconsiderthebn, d+1≤n≤2d−1. Wemustestimatebd+n(x,Dprime/Dd)DnDd ϕν.By (2.11) we havebardblDnDd ϕνbardblL2(Rn)≤C132−|ν|/2. Therefore,bardblbd+n(x,Dprime/Dd)DprimeDdϕνbardblL2(Rd)≤C142−|ν|/2. (A.3)From (A.1) and the estimates (A.2) and (A.3) for the summands, we find156Appendix A. Proofs of lemma 1 and theorem 1the result.157Appendix BReverse time wave equationmigrationThis appendix is part of chapter 5. Most migration operators are definedto be the adjoint of what is called the scattering operator. This assumptionis also true for reverse time migration.B.1 Single scatteringThe causal acoustic Green’s function G(x,t;xs),x∈R3 for a point sourceat x = xs is the solution of1v2(x)∂2G∂t2 (x,t;xs)−∇2xG(x,t;xs) = δ(x−xs)δ(t) (B.1)with G = 0,t < 0 and v is the acoustic wave velocity field.Denote by m(x) = δv(x)/v(x) a relative perturbation of the velocityfield. Then linearization of the wave equation yields for the correspondingperturbation of the Green’s functionA version of this appendix will be submitted for publication. Moghaddam, P.P.,Herrmann F. J. and Shahidi, R. (2010) True amplitude depth migration using curvelets.158B.2. Shot-geophone modeling and migration1v2(x)∂2δG∂t2 (x,t;xs)−∇2xδG(x,t;xs) =2m(x)v2(x)∂2∂t2G(x,t;xs) (B.2)whose solution has the integral representation at the source and receiverpoints xr,xsδG(xr,t;xs) = ∂2∂t2integraldisplaydxintegraldisplaydτ2m(x)v2(x) G(x,t−τ;xr)G(x,τ;xs) (B.3)B.2 Shot-geophone modeling and migrationThe single-scattered wave field is the time convolution of δG with a sourcewavelet. The main concern of this paper is the kinematic relationshipsbetween data and image, thus we ignore the filtering effect of the sourcefunctional and replace it by delta function. This replacement of the sourceby an impulse does not violate any of our assumptions regarding the adjointstate method, thus the Born modeling operator K[v] isK[v]m(x) = δG(xr,t;xs) (B.4)The crux of our amplitude recovery method relies on the shot-geophonemigration operator to be the adjoint of the shot-geophone modeling opera-tor. The derivation of the adjoint reverse time implementation is a minorvariation on the usual implementation of reverse time migration (the ’ad-joint state method’, ((see e.g. Tarantola, 1987; Whitmore, 1983; Lailly,1983; Yoon and Hong, 2003; Symes and Stolk, 2004; Symes, 2007)). The159B.2. Shot-geophone modeling and migrationresult isK∗[v]d(xr,t;xs) = ˆm(x) =−integraldisplaydxsintegraldisplay T0dt2v(x)q(x,t;xs)∂2G∂t2 (x,t;xs)(B.5)where the adjoint state or backpropagated field q(x,t;xs) satisfies q = 0,t >T and1v2(x)∂2q∂t2 (x,t;xs)−∇2xq(x,t;xs) =integraldisplaydxrd(xr,t;xs)δ(x−xr) (B.6)The migration operator defined by above equations is the reverse timemigration operator. Symes and Stolk (2004) showed that the migrationoperator which is defined in equations B.5 and B.6 is the adjoint of themodeling operator defined in equation B.3. The reverse time migration thatis used for this work is based on the Symes and Stolk (2004) implementation,forwhichtheadjointnessofthemodelingandmigrationoperatorsisproperlytested in the discrete sense. By having the migration and modeling operatorsproperly set, we can proceed with our amplitude recovery method.160BibliographyLailly, P., 1983, The seismic inverse problem as a sequence of before stackmigrations: SIAM, 206–220.Symes, W., 2007, Reverse time migration with optimal checkpointing: Geo-physics, 72, 213–221.Symes, W. and C. Stolk, 2004, Reverse time shot-geophone migration: TheRice Inversion Project, Department of Computational and Applied Math-ematics, Rice University, Houston, Texas, USA.Tarantola, A., 1987, Inverse problem theory: Elseweir.Whitmore, N. D., 1983, Iterative depth migration by backward time propa-gation: Presented at 53rd Annual International SEG Meeting, Las Vegas.Yoon, K. and S. Hong, 2003, 3D reverse-time migration using hte acous-tic wave equation: An experience with the SEG/EAGE data set.: TheLeading Edge, 22–38.161
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Curvelet-based migration amplitude recovery
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Curvelet-based migration amplitude recovery P. Moghaddam, Peyman 2010
pdf
Page Metadata
Item Metadata
Title | Curvelet-based migration amplitude recovery |
Creator |
P. Moghaddam, Peyman |
Publisher | University of British Columbia |
Date | 2010 |
Date Issued | 2010-05-04 |
Description | Migration can accurately locate reflectors in the earth but in most cases fails to correctly resolve their amplitude. This might lead to mis-interpretation of the nature of reflector. In this thesis, I introduced a method to accurately recover the amplitude of the seismic reflector. This method relies on a new transform-based recovery that exploits the expression of seismic images by the recently developed curvelet transform. The elements of this transform, called curvelets, are multi-dimensional, multi-scale, and multi-directional. They also remain approximately invariant under the imaging operator. I exploit these properties of the curvelets to introduce a method called Curvelet Match Filtering (CMF) for recovering the seismic amplitude in presence of noise in both migrated image and data. I detail the method and illustrate its performance on synthetic dataset. I also extend CMF formulation to other geophysical applications and present results on multiple removal. In addition of that, I investigate preconditioning of the migration which results to rapid convergence rate of the iterative method using migration. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | Eng |
Collection |
Electronic Theses and Dissertations (ETDs) 2008+ |
Date Available | 2010-05-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial 3.0 Unported |
DOI | 10.14288/1.0052987 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc/3.0/ |
URI | http://hdl.handle.net/2429/24421 |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_2010_fall_poormoghaddam_peyman.pdf [ 7.94MB ]
- Metadata
- JSON: 1.0052987.json
- JSON-LD: 1.0052987+ld.json
- RDF/XML (Pretty): 1.0052987.xml
- RDF/JSON: 1.0052987+rdf.json
- Turtle: 1.0052987+rdf-turtle.txt
- N-Triples: 1.0052987+rdf-ntriples.txt
- Original Record: 1.0052987 +original-record.json
- Full Text
- 1.0052987.txt
- Citation
- 1.0052987.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Japan | 7 | 0 |
China | 7 | 0 |
United States | 5 | 0 |
India | 3 | 0 |
Russia | 3 | 0 |
Egypt | 1 | 0 |
City | Views | Downloads |
---|---|---|
Tokyo | 7 | 0 |
Ashburn | 5 | 0 |
Unknown | 4 | 0 |
Qingdao | 3 | 0 |
Pune | 3 | 0 |
Shanghai | 2 | 0 |
Guangzhou | 2 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052987/manifest