UBC Faculty Research and Publications

Compressive sampling meets seismic imaging Herrmann, Felix J. 2007-12-31

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

Item Metadata


herrmann07PIMS.pdf [ 10.1MB ]
JSON: 1.0102592.json
JSON-LD: 1.0102592+ld.json
RDF/XML (Pretty): 1.0102592.xml
RDF/JSON: 1.0102592+rdf.json
Turtle: 1.0102592+rdf-turtle.txt
N-Triples: 1.0102592+rdf-ntriples.txt
Original Record: 1.0102592 +original-record.json
Full Text

Full Text

Compressive sampling meets seismic imagingFelix J. Herrmann joint work withTim Lin*, Peyman Moghaddam*, Gilles Hennenfent*,Deli Wang* & Chris Stolk (Universiteit Twente)*Seismic Laboratory for Imaging and Modeling slim.eos.ubc.caPIMS-CINVESTAV 2007, Monterrey, October 19Today?s challengesAside from spurious local minima seismic waveform inversion is difficult because of? lack of control on the image amplitudes? missing data and noise? computational cost to form the operatorsToday?s agenda is to leverage recent insights from applied harmonic analysis and information theory to? restore amplitudes => affordable q-Newton updates? stably reconstruct wavefields ? compress wavefield-extrapolation operators MotivationExploit two aspects of curvelets, namely their? parsimoniousness? invariance under certain operatorsFormulate? data-adaptive scaling algorithms ? non-adaptive wavefield reconstruction algorithmsApplications? nonlinear migration-amplitude recovery? nonlinear sampling for wavefields? nonlinear sampling for operatorsToday?s topicsSparsity-promoting seismic-image amplitude recovery? curvelet-domain diagonal approximation of PsDO?s? stable sparsity-promoting inversionDirectional frame-based wavefield reconstruction by sparsity promotion? curvelet parsimoniousness? jitter samplingCompression of FIO?s through compressive sampling? measurement basis diagonalizes operatorThe problemMinimization:After linearization (Born app.) forward model with noise:Conventional imaging:    is prohibitively expensive to invertevaluation of         involves expensive wavefield extrapolators2-D curveletscurvelets are of rapid decay in spacecurvelets are strictly localized in frequencyx-t f-kOscillatory in one direction and smooth in the others!Obey parabolic scaling relationCoefficients Amplitude Decay In Transform DomainsFourierWaveletsCurveletsPartial ReconstructionFourier (1% largest coefficients)SNR = 2.1 dBPartial ReconstructionCurvelets (1% largest coefficients)SNR = 6.0 dBCurvelets are oscillatory in one direction and smooth in the others.3-D curveletsApproximate linearized inversion by curvelet scaling & sparsity promotionJoint work with Chris Stolk* and Peyman MoghaddamMathematics Department, Twente University, the Netherlands?Sparsity- and continuity-promoting seismic imaging with  curvelet frames? to appear in ACHARelated workWavelet-Vaguelette/Quasi-SVD methods based on? homogeneous operators? absorb ?square-root? of the Gramm matrix in WVD?s? Wavelets/curvelets near diagonalize the operator and are sparse on the model? Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition (Donoho ?95)? Recovering Edges in Ill-posed Problems: Optimality of curvelet Frames (Candes & Donoho ?00)Scaling methods based on a diagonal approximation of    , assuming? smoothness on the symbol and conormality reflectors? Illumination-based normalization (Rickett ?02)? Amplitude preserved migration (Plessix & Mulder ?04)? Amplitude corrections (Guitton ?04)? Amplitude scaling (Symes ?07)Hessian/Normal operator[Stolk 2002, ten Kroode 1997, de Hoop 2000, 2003]Alternative to expensive least-squares migration.In high-frequency limit     is a pseudo-differential operator? composition of two Fourier integral operators? pseudolocal (near unitary)? singularities are preserved? symbol is smooth for smooth velocity modelsCorresponds to a spatially-varying dip filter after appropriate preconditioning (=> zero-order PsDO).? Substitutions: orwithApproximationAllows for decomposition of the normal operatorScalingMatching procedureCompute reference vector <=> defines g? migrate data ? apply spherical-divergence correctionCreate ?data? <=> defines f? demigrate? migrateEstimate scaling by inversion procedureDefine scaled curvelet transformRecover migration amplitudes by sparsity promotion.Key ideaEstimation curvelet-domain scaling? inversion of an underdetermined system? over fitting? positivity and reasonable scalingSolution:? use smoothness of the symbol ? formulate nonlinear estimation problem that minimizeswith? solve with l-BFGS [Noccedal, Symes ?07]Key ideaEast quadrantsWest quadrantsNorth quadrantsSouth quadrants16 angles/quad8 angles/quadFine scalescoarserscalesKey ideaImpose smoothness via following system of equationswithfirst-order differences in space and angle directions for each scale. Equivalent towithSmoothness penaltyincreasing smoothness? reduces overfitting? scaling is positive and reasonableSmoothness penaltySmoothness penaltyOur approach?Forward? model:? diagonal approximation of the demigration-migration operator? costs one demigration-migration to estimate the diagonal weightingSolutionSolveMigrated data Amplitude-corrected & denoised migrated dataImaging example? two-way reverse time wave-equation migration with checkpointing [Symes ?07]? adjoint state method with 8000 time steps? evaluation       takes 6 h on 60 CPU?sObservationsCurvelet-domain scaling? handles conflicting dips (conormality assumption)? exploits invariance under the PsDODiagonal approximation? exploits smoothness of the symbol? uses ?neighbor? structure of curveletsResults on the SEG AA? show? recovery of amplitudes beneath the Salt? successful recovery from clutter? improvement of the continuity ? robust w.r.t. noiseCurvelet-domain matched filter ...A primer on compressive samplingCompressive sensing[Candes, Romberg & Tao, Donoho, many others]Three key ingredients? existence of a sparsifying transform? handle wavefronts & reflectors with conflicting dips? existence of a sub-Nyquist sampling strategy that reduces coherent aliases? incoherence ? random sampling scheme? existence of a large-scale (norm-one) solver? sparsity promotion by iterative thresholding and coolingSeismic Laboratory for Imaging and ModelingSimple exampleFouriertransform??3-fold under-samplingsignificant coefficients detectedambiguityfew significant coefficientsFouriertransformFouriertransformSeismic Laboratory for Imaging and ModelingForward problem=Fourier coefficients(sparse)withFouriertransformrestrictionoperatorsignalSeismic Laboratory for Imaging and ModelingNaive sparsity-promoting recoveryinverseFouriertransformdetection +data-consistentamplitude recoveryFouriertransform==detectiondata-consistent amplitude recovery=Seismic Laboratory for Imaging and Modeling? ?noise?? due to AHA ? I? defined by AHAx0-?x0 = AHy-?x0Undersampling ?noise?less acquired data3 detectable Fourier modes 2 detectable Fourier modes1 out of 2 1 out of 4 1 out of 6 1 out of 8Seismic Laboratory for Imaging and ModelingSparsity-promoting wavefield reconstruction= withsparsifying transformfor seismic datarestriction operator[Sacchi et al ?98][Xu et al ?05][Zwartjes and Sacchi ?07][Herrmann and Hennenfent ?07]complete wavefield (transform domain)acquireddataInterpolated data given by                 withSeismic Laboratory for Imaging and ModelingObservations? bla bla? generalized to A=RMS^H? depends on solver, sampling strategy and sparsity transformCompressive sampling of wavefieldsjoint work with Deli Wang (visitor from Jilin university) and Gilles Hennenfent?Curvelet-based seismic data processing: a multiscale and nonlinear approach? & to appear in Geophysics, ?Non-parametric seismic data recovery with curvelet frames? and ?Simply denoise: wavefield reconstruction via jittered undersampling?Solution ofrecovers the function f. General form compressive samplingThe problemTotal data 85 % traces missingRequirementsSparsifying transform (S)? curvelet? focussed curvelets Sampling scheme (RM)? random sampling? random jittered sampling => control largest gapsSparsity promoting solver (P)? Iterative thresholding (Landweber + soft threshold)Seismic Laboratory for Imaging and ModelingDiscrete random jittered undersamplingreceiverpositionsreceiverpositionsPDFreceiverpositionsPDFreceiverpositionsPDF[Hennenfent and Herrmann ?07]Typical spatial convolution kernel(amplitudes)Averaged spatial convolution kernel(amplitudes)Sampling schemepoorlyjitteredoptimallyjitteredrandomregularSolution ofrecovers the wavefield f. Curvelet-based recoverySeismic Laboratory for Imaging and ModelingModelSeismic Laboratory for Imaging and ModelingRegular 3-fold undersamplingSeismic Laboratory for Imaging and ModelingCRSI from regular 3-fold undersamplingSNR = 6.92 dBSeismic Laboratory for Imaging and ModelingRandom 3-fold undersamplingSeismic Laboratory for Imaging and ModelingCRSI from random 3-fold undersamplingSNR = 9.72 dBSeismic Laboratory for Imaging and ModelingOptimally-jittered 3-fold undersamplingSeismic Laboratory for Imaging and ModelingCRSI from opt.-jittered 3-fold undersamplingSNR = 10.42 dBSeismic Laboratory for Imaging and ModelingModelSeismic Laboratory for Imaging and ModelingRegular 3-fold undersamplingSNR = 12.98 dBSeismic Laboratory for Imaging and ModelingSNR = 12.98 dBRegular 3-fold undersamplingSeismic Laboratory for Imaging and ModelingOptimally-jittered 3-fold undersamplingSNR = 15.22 dBSeismic Laboratory for Imaging and ModelingOptimally-jittered 3-fold undersamplingSNR = 15.22 dBSolution ofrecovers the wavefield f. Focussed recovery80 % missingFocused curvelet recoveryCurvelet recoveryOriginal dataObservationsRegular subsampling is unfavorable? random sampling favorable but suffers from gaps? jitter sampling favorable and controls gapsFocal transform ? is reminiscent of an imaging operator? improves recovery <=> additional compressionSolver? solves norm one problem for 200-300 matrix-vector multiplications for 230 unknowns ...Outlook? Migration-based wavefield reconstruction? sparsity on the image? focussing of the image (extra constraint)? or a more ?blue sky? approach of compressive one-way wavefield extrapolationCompressed wavefield extrapolationjoint work with Tim Lin?Compressed wavefield extrapolation? in GeophysicsMotivationSynthesis of the discretized operators form bottle neck of imagingOperators have to be applied to multiple right-hand sidesExplicit operators are feasible in 2-D and lead to an order-of-magnitude performance increaseExtension towards 3-D problematic? storage of the explicit operators? convergence of implicit time-harmonic approachesFirst go at the problem using CS techniques to compress the operator ...Related workCurvelet-domain diagonalization of FIO?s? The Curvelet Representation of Wave Propagators is Optimally Sparse (Candes & Demanet ?05)? Seismic imaging in the curvelet domain and its implications for the curvelet design (Chauris ?06)? Leading-order seismic imaging using curvelets (Douma & de Hoop ?06)Explicit time harmonic methods? Modal expansion of one-way operators in laterally varying media (Grimbergen et. al. ?98)? A new iterative solver for the time-harmonic wave equation (Riyanti ?06) Fourier restriction? How to choose a subset of frequencies in frequency-domain finite-difference migration (Mulder & Plessix ?04)signal in space domainsignal in space domainF L1L1incomplete signal in Fourier domainincomplete and shifted signal in Fourier domain shifted signal in space domainsignal in space domainCompressed SensingCompressed ProcessingFSuppose we want to shift a sparse spike train, i.e.,? Eigen modes <=> Fourier transform.? Can this operation be compressed by compressive sampling?InspirationCalculate instead? Take compressed measurements in Fourier space.? Recover with sparsity promotion? Shift operator is compressed by the restrictionyielding compressed rectangular operators.? Extend this idea to wavefield extrapolation?Operators on spikes[Candes et. al, Donoho]Representation for seismic data[Berkhout]s+W-R+W+p-Different representationsdiagonalizationoperatorparsimonywavefieldSVD/Lanczos/modal? ?curvelets ? ?Different representationsdiagonalizationoperatorparsimonywavefieldSVD/Lanczos/modal? ?curvelets ? ?If incoherent this may actually work ....Sparsity promoting formulationBuys us stability w.r.t. missing data? provided measurement and sparsity representations are mutually incoherent? sufficient mixing <=> random restrictionDifferent strategy:? Let the physics define the measurement basis? Use the modal domain (domain of eigenfunctions) to define the measurement basis? See what you can recoverStudy eigenfunctions:? mutual coherence with sparsity representation? modal spectrum on the to-be-extrapolated wavefieldOne-Way Wave Operator? Structure of     confounds the meaning of its exponentiation, due to it being an operator?    contains information about medium velocityTwo-wayWave Operator(Simon & Reed; Dessing ?97; Grimbergen ?98)One-Way Wave Operator? Solution of the one-way wave equation? After discretization solve eigenproblem on ? Helmholtz operator is Hermitian? monochromatic? velocity   varies laterally(Claerbout, 1971; Wapenaar and Berkhout, 1989)Modal transform? Solve eigenproblem & take square root?    is orthonormal & defines the modal transform that diagonalizes one-way wavefield extrapolation? Eigenvalues play role of vertical wavenumbers? Extrapolation operator is diagonalizedCompressed wavefield extrapolationRecorded DataOriginal eventsForward modelReconstruct point scatterers  from recorded data ....Compressed wavefield extrapolation? Randomly subsample & phase rotate in Modal domain? Recover by norm-one minimization? Capitalize on ? the incoherence modal functions and point scatterers? reduced explicit matrix size? constant velocity <=> Fourier recoveryCompressed wavefield extrapolationRecorded Data Reconstructed eventsReconstructionOnly 1 % of original modes were used ...? Despite the existence of evanescent (exponentially decaying) waves modes recovery is successful? If you are looking for point-scatterers, we have a proof of concept that is fast? Earth is more complex ...ObservationsCompressed wavefield extrapolation? Extend to general wavefields? Use curvelets as the sparsity representation? Use the full & compressed forward operator operator? Compressively extrapolate back 600m to the sourceRestriction & sparsity strategies? Forward extrapolation:? Inverse extrapolation:Forward Extrapolation? (a) is Full extrapolation? (b)-(d) is compressed extrapolation, (b) p = 0.04, (c) p = 0.16, (d) p = 0.24Inverse Extrapolation? (a) p = 0.04? (b) p = 0.16, (c) pf=0.4, px=0.4Evanescent Recovery? (a) is downward extrapolated wavefield? (b) is matched filter? (c) is ?compressed? inverse extrapolationVelocity modelCompressed inverse extrapolationOverthrust exploding reflector Full forward extrapolation Matched filter Recovered from p=0.25Multiscale and angular compressed wavefield extrapolation? Propose a scheme motivated by extensions of CS ? adapt discretization & restriction? parallel implementation(Tsaig and Donoho ?06)Conclusions? Curvelets sparsity on the model and near diagonalization yields stable inversion Gramm matrix? Jittered sampling and focussing in combination with curvelets leads to wavefield recovery? Compressed wavefield extrapolation? reduction in synthesis cost? inverse extrapolation works well when focussed? mutual coherence curvelets and modes? performance of norm-one solver? keep the constants under control ...? Double-role CS matrix is cool ... upscaling to ?real-life? is a challenge ....Open problems? What deeper insights can CS give?? inversion near unitary operators? coherency generalized to frames to study? cols modeling operator <=> curvelets? radiation vs guided modes <=> curvelets? Norm-one solver for reduced system as fast a LSQR on the full system? Fast random eigenvalue solver does not exist yet ...? Extension of CS to waveform inversion & to compressed computations ...? Many more ...AcknowledgmentsThe audience for listening and the organizers for putting this great workshop together ....The authors of CurveLab (Demanet,Ying,  Donoho)Dr. W.W. or his reverse-time migration code This work was in part financially supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant (22R81254) and the Collaborative Research and Development Grant DNOISE (334810-05) of F.J.H.This research was carried out as part of the SINBAD project with support,secured through ITF (the Industry Technology Facilitator), om the following organizations: oup, , vron,ExxonMobil and Shell.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics

Country Views Downloads
China 13 44
Brazil 8 0
United States 8 0
Russia 6 0
Germany 5 3
Canada 1 0
Japan 1 0
City Views Downloads
Unknown 14 3
Shenzhen 10 44
Saint Petersburg 6 0
Ashburn 4 0
Beijing 3 0
Sunnyvale 2 0
Redmond 1 0
Ottawa 1 0
Tokyo 1 0

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



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items