UBC Faculty Research and Publications

Seismic imaging and processing with curvelets Herrmann, Felix J. 2007

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

Item Metadata


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

Full Text

Seismic imaging and processing with curvelets Felix J. Herrmann joint work with Deli Wang  Combinations of parsimonious signal representations with nonlinear sparsity promoting programs hold the key to the next-generation of seismic data processing algorithms ... Since they   allow for formulations that are stable w.r.t.     noise incomplete data moderate phase rotations and amplitude errors  Finding a sparse representation for seismic data & images is complicated because of     wavefronts & reflectors are multiscale & multidirectional the presence of caustics, faults and pinchouts  The curvelet transform  Representations for seismic data Transform  Underlying assumption  FK  plane waves  linear/parabolic Radon transform  linear/parabolic events  wavelet transform  point-like events (1D singularities)  curvelet transform  curve-like events (2D singularities)  Properties curvelet transform:     multiscale: tiling of the FK domain into dyadic coronae multi-directional: coronae subpartitioned into angular wedges, # of angle doubles every other scale    anisotropic: parabolic scaling principle    Rapid decay space    Strictly localized in Fourier    Frame with moderate redundancy (8 X in 2-D and 24 X in 3-D)  k2 2j/2  angular wedge  2j  fine scale data  k1  coarse scale data  2-D curvelets  curvelets are of rapid decay in space  x-t  curvelets are strictly localized in frequency  f-k  Oscillatory in one direction and smooth in the others! Obey parabolic scaling relation length ≈ width2  Curvelet tiling & seismic data Angular wedge  Curvelet tiling  # of angles doubles every other scale doubling!  Real data frequency bands example  Data is multiscale!  1  2  3  4  5  6  7  8  Single frequency band angular wedges  6  6th scale image  Data is multidirectional! Seismic Laboratory for Imaging and Modeling  Decomposition in angular wedges  Wavefront detection 0  -2000  Offset (m) 0  2000  Time (s)  0.5  1.0  1.5  2.0  Significant curvelet coefficient  Curvelet coefficient~0  Extenstion to 3-D Cartesian Fourier space [courtesy Demanet ‘05, Ying ‘05]  Curvelets live in a wedge in the 3 D Fourier plane...  3-D curvelets  Curvelets are oscillatory in one direction and smooth in the others.  Coefficients Amplitude Decay In Transform Domains  Fourier Wavelets Curvelets  Partial Reconstruction Curvelets (1% largest coefficients)  SNR = 6.0 dB  Curvelet sparsity promotion  Forward model Linear model for the measurements of a function m0:  y  = Km0 + n  with y  = data  K  = the modeling matrix  m0  = the model vector  n = noise    inversion of K either ill-posed or underdetermined. seek a prior on m.  Key idea ˜ = arg min x x x  1  sparsity enhancement  s.t.  Ax − y  2  ≤  data misfit  When a traveler reaches a fork in the road, the l1 -norm tells him to take either one way or the other, but the l2 -norm instructs him to head off into the bushes. John F. Claerbout and Francis Muir, 1973 New field “compressive sampling”: D. Donoho, E. Candes et. al., M. Elad etc. Preceded by others in geophysics: M. Sacchi & T. Ulrych and co-workers etc.  Inline 0  Linear quadratic (lsqr):  •  2  s.t.  Ax − y  2  ≤  model Gaussian  Non-linear  s.t.  Ax − y  2  ≤  0.5  model Cauchy (sparse)  Problem:  •  500  300  400  500  Inline  Time (s)  •  200  400  1.5  : 1  100  300  1.0  2.0 0  ˜ = arg min x x x  200  0.5  Time (s)  ˜ = arg min x x x  100  1.0  data does not contain point scatterers  1.5  not sparse  2.0  Our contribution Inline 0  Model as superposition of little plane waves. Time (s)  Compound modeling operator with curvelet synthesis:  0.5  1.0  1.5  K m0 ˜ m  T  →  KC  =  ˜ C x  →  x0 T  Exploit parsimoniousness of curvelets on seismic data & images ...  2.0  100  200  300  400  500  Sparsity-promoting program Problems boils down to solving for x0 signal  y  =  A  +  n  noise  x0 curvelet representation of ideal data  with  P :     ˜ = arg minx x x ˜ = CT x ˜ m  1  s.t.  Ax − y  2  ≤  exploit sparsity in the curvelet domain as a prior find the sparsest set of curvelet coefficients that T ˜ match the data, i.e., y ≈ KC x invert an underdetermined system  Solver Initialize: i = 0; x0 = 0; Choose: L, AT y while y − Axi  2  ∞  >  > λ1 > λ2 > · · · do  for l = 1 to L do xi+1 = Tλsi xi + AT y − Axi end for i = i + 1; end while f = CT xi .  Applications Problems in seismic processing can be cast in to P    stable under noise stable under missing data  Obtain a formulation that     explicitly exploits compression by curvelets is stable w.r.t. noise exploits the “invariance” of curvelets under imaging  Applications include     seismic data regularization primary-multiple separation seismic amplitude recovery  Seismic data regularization joint work with Gilles Hennenfent  Motivation  Irregular sub-sampling  incoherent noise  Noisy because of irregular sampling ...  Sparsity-promoting inversion* Reformulation of the problem signal  y  =  H  RC  +  n  noise  x0 curvelet representation of ideal data  Curvelet Reconstruction with Sparsity-promoting Inversion (CRSI)   P :  look for the sparsest/most compressible, physical solution KEY POINT OF THE   RECOVERY data misfit sparsity constraint sparsity constraint         = arg ˜ H H − PC minmin x 0x 0 s.t. s.t. y −yPC x 2x≤2!≤ ! x =x˜arg ˜0 )= arg minxx Wx x s.t. Ax − y 2 ≤ (P0 )(Px 1     T  ˜   ˜H f = C x  H ˜f = ˜fC= xC ˜ x˜  * inspired by Stable Signal Recovery (SSR) theory by E. Candès, J. Romberg, T. Tao, Compressed sensing by D. Donoho & Fourier Reconstruction with Sparse Inversion (FRSI) by P. Zwartjes  Original data  80 % missing  CRSI recovery with 3-D curvelets  Primary multiple separation Joint work with Eric Verschuur, Deli Wang, Rayan Saab and Ozgur Yilmaz  Motivation Primary-multiple separation step is crucial    moderate prediction errors 3-D complexity & noise  Inadequate separation leads to    remnant multiple energy deterioration primary energy  Introduce a transform-based technique    stable insensitive to moderate shifts & phase rotations  Exploit sparsity and parameterization transformed domain  Move-out error  Move-out error  now two sparsityThe matrices y =problem Ax0 one + n,for each signal co Sparse signal model:  y = Ax0 + n, with  A = [A1   A2 ]  and  x0 = [x01  augmented synthesis and sparsity vectors index11 <->2primary 0 01 index 2 <-> multiple  T  x02 ] T  A = [A A and ] and xvector, = [x respectivel x02 ] ynthesis matrix sparsity    erved for primaries and multiples. The above s  se weights drive the two signal components apart during the optimization ‡ accuracy .  The solution The w-weighted optimization problem becomes  to reasonable The weighted norm-one       minx x w,1      Pw : sˆ1 = A1 xˆ 1          given: s˘2  [w1 ,  optimization problem:  subject to  y − Ax  2  ≤ε  and sˆ2 = A2 xˆ 2 and w(y, s˘2 )  with T  w2 ] the weighting vectors withT strictly positive weights defined in  w  :=  w1 , w2  T primaries T multiples. TheA estimates for the and multiples are computed from := C , C  s.2During := thepredicted multiples minimizes Pw˘ optimization, the sparsity vector is recovered b ˘s1  ˘2 := S − S  Solution cont’d The weights  with ˘ 1 ≈ C˘s1 u ˘ 2 ≈ C˘s2 u     √ w1 := max σ · 2 log N , C1 |˘ u1 | √ w2 := max σ · 2 log N , C2 |˘ u2 |  during minimization signal components are driven apart curvelet compression helps separates on the basis of position, scale and direction  Synthetic example  total data  SRME predicted multiples  Synthetic example  SRME predicted primaries  curvelet-thresholded  Synthetic example  SRME predicted primaries  estimated  Real example total data  SRME predicted multiples  Real example curvelet thresholded  curvelet estimated  SRME predicted primaries  curvelet estimated primaries  Seismic amplitude recovery Joint work with Chris Stolk and Peyman Moghaddam  Motivation Migration generally does not correctly recover the amplitudes. Least-squares migration is computationally unfeasible. Amplitude recovery (e.g. AGC) lacks robustness w.r.t. noise. Existing diagonal amplitude-recovery methods    do not always correct for the order (1 - 2D) of the Hessian [see Symes ‘07] do not invert the scaling robustly  Moreover, these (scaling) methods assume that there  are no conflicting dips (conormal) in the model     is infinite aperture are infinitely-high frequencies etc.  Existing scaling methods Methods are based on a diagonal approximation of Ψ.      Illumination-based normalization (Rickett ‘02) Amplitude preserved migration (Plessix & Mulder ‘04) Amplitude corrections (Guitton ‘04) Amplitude scaling (Symes ‘07)  We are interested in an ‘Operator and image adaptive’ scaling method which       estimates the action of Ψ from a reference vector close to the actual image assumes a smooth symbol of Ψ in space and angle does not require the reflectors to be conormal <=> allows for conflicting dips stably inverts the diagonal  Our approach “Forward” model:  y  ≈ Ax0 + ε  with y  =  A  := C Γ  T  AA r K    = K Km + ε T  migrated data T  T  ≈  K Kr  =  the demigration operator  =  migrated noise.  diagonal approximation of the demigration-migration operator costs one demigration-migration to estimate the diagonal weighting  Solution  Solve  P: with    minx J(x) subject to  y − Ax  2    ˜ = (AH )† x ˜ m sparsity  J(x) = α x  1  +β Λ  1/2  A  H  †  continuity  x  p  .  ≤  Example SEGAA’ data:     “broad-band” half-integrated wavelet [5-60 Hz] 324 shots, 176 receivers, shot at 48 m 5 s of data  Modeling operator     Reverse-time migration with optimal check pointing (Symes ‘07) 8000 time steps modeling 64, and migration 294 minutes on 68 CPU’s  Scaling requires 1 extra migration-demigration  Migrated data  Amplitude-corrected & denoised migrated data  Noise-free data  Noisy data (3 dB)  Data from migrated image  Data from amplitude-corrected & denoised migrated image  Nonlinear data  Conclusions The combination of the parsimonious curvelet transform with nonlinear sparsity & continuity promoting program allowed us to...     recover seismic data from large percentages missing traces separate primaries & multiples recover migration amplitudes  This success is due to the curvelet’s ability to     detect wavefronts <=> multi-D geometry differentiate w.r.t. positions, angle(s) and scale diagonalize the demigration-migration operator  Because of their parsimoniousness on seismic data and images, curvelets open new perspectives on seismic processing ...  Acknowledgments The authors of CurveLab (Demanet,Ying, Candes, Donoho) William Symes for the reverse-time migration code. These results were created with Madagascar developed by Sergey Fomel. 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), from the following organizations: BG Group, BP, Chevron,ExxonMobil and Shell.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items