Seismic data processing with curvelets: a multiscale and nonlinear approach Felix J. Herrmann joint work with Deli Wang, Gilles Hennenfent and Peyman Moghaddam Motivation Exploit two aspects of curvelets, namely their  parsimoniousness  invariance under certain operators Formulate  non-adaptive wavefield reconstruction algorithms  data-adaptive matching algorithms Applications  nonlinear sampling theory for wavefields  nonlinear migration-amplitude recovery  nonlinear primary-multiple separation Approach Employ parsimoniousness by sparsity promotion. Exploit behavior of certain operators in phase space  diagonalization <=> curvelet domain scaling  smoothness <=> structure of phase space Combine parsimoniousness with structure in phase space  diagonal approximation operators  stable amplitude recovery  improved adaptive separation Migration-amplitude recovery methods are based on  diagonal approximation of Pseudo’s  estimate scaling from a reference vector and demigrated-migrated reference vector  Illumination-based normalization (Rickett ‘02)  Amplitude corrections (Guitton ‘04)  Amplitude scaling (Symes ‘07) Primary-multiple separation methods are based on  diagonal approximation in the Fourier domain  estimate scaling from mismatch pred. multiples & data  adaptive subtraction (Verschuur and Berkhout ‘97) We are interested in a formulation that  estimates the scaling with smoothness control  prevents overfitting  allows for conflicting dips  incorporates curvelet-domain sparsity promotion The curvelet transform 2-D curvelets curvelets are of rapid decay in space curvelets are strictly localized in frequency x-t f-k Oscillatory in one direction and smooth in the others! Obey parabolic scaling relation length ≈ width2 Coefficients Amplitude Decay In Transform Domains Fourier Wavelets Curvelets Partial Reconstruction Fourier (1% largest coefficients) SNR = 2.1 dB Partial Reconstruction Curvelets (1% largest coefficients) SNR = 6.0 dB Non-adaptive curvelet- domain sparsity promotion 00.5 1.0 1.5 2.0 Tim e ( s) 100 200 300 400 500 Inline Linear quadratic (lsqr): • model Gaussian Non-linear: • model Cauchy (sparse) Problem: • data does not contain point scatterers • not sparse 0 0.5 1.0 1.5 2.0 Tim e ( s) 100 200 300 400 500 Inline x̃ = arg min x ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! x̃ = arg min x ‖x‖2 s.t. ‖Ax− y‖2 ≤ ! Our contribution Model as superposition of little plane waves. Compound modeling operator with curvelet synthesis: Exploit parsimoniousness of curvelets on seismic data & images ... 0 0.5 1.0 1.5 2.0 Ti m e (s ) 100 200 300 400 500 Inline K !→ KCT m0 !→ x0 m̃ = CT x̃ Sparsity-promoting program Problems boil down to solving for with  exploit sparsity in the curvelet domain as a prior  find the sparsest set of curvelet coefficients that match the data, i.e.,  invert an underdetermined system signal =y + n noise curvelet representation of ideal data x0 A x0 P! : { x̃ = arg minx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! m̃ = CT x̃ y ≈ KCT x̃ Seismic wavefield reconstruction with CRSI Sparsity-promoting inversion* Reformulation of the problem Curvelet Reconstruction with Sparsity-promoting Inversion (CRSI)  look for the sparsest/most compressible, physical solution KEY POINT OF THE RECOVERY * 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 signal =y + n noise curvelet representation of ideal data x0 RCH (P0)   x̃= arg sparsity constraint︷ ︸︸ ︷ min x ‖x‖0 s.t. ‖y−PC Hx‖2 ≤ ! f̃= CH x̃ (P0)   x̃= arg sparsity constraint︷ ︸︸ ︷ min x ‖x‖0 s.t. data misfit︷ ︸︸ ︷ ‖y−PCHx‖2 ≤ ! f̃= CH x̃ P! : { x̃ = arg minx ‖Wx‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ = CT x̃ Original data 80 % missing CRSI recovery with 3-D curvelets Adaptive curvelet- domain matched filtering Forward model Linear model for amplitude mismatch:  spatially-varying dip filter  zero-order Pseudo After discretization  linear operator  f and g known  matrix B is full and not known .... f = Bg ( Bf ) (x) = ∫ x∈Rd ejk·xb(x, k)f̂(k)dk B = Pseudodifferential operator b(x, k) = the symbol Forward model Diagonal approximation in the curvelet domain:  curvelet domain scaling  opens the way to an estimation of w Examples: f = Bg ≈ CTdiag{w}Cg B f g migration migrated “image” “reflectivity” multiple removal obliquity factor total data predicted multiples KTK Key idea Problems with estimating w  inversion of an underdetermined system  over fitting  positivity and reasonable scaling by w Solution:  use smoothness of the symbol  formulate nonlinear estimation problem that minimizes with  solve with l-BFGS Jγ(z) = 1 2 ‖d− Fγez‖22, gradJ(z) = diag{ez}[FT (Fez − d)] Key idea East quadrants West quadrants North quadrants South quadrants 16 angles/ quad 8 angles/ quad x1 x2 θ Fine scales coarser scales D1 D2 Dθ Key idea Impose smoothness via following system of equations with first-order differences in space and angle directions for each scale. Equivalent to with f = CTdiag{Cg}w 0 = γLw L = [ DT1 D T 2 D T θ ]T w̃ = arg minw 1 2 ‖b−P[w]‖22 + γ2‖Lw‖22 P = CTdiag{Cg} Smoothness penalty increasing smoothness  reduces overfitting  scaling is positive and reasonable Smoothness penalty 50 100 150 200 20 40 60 80 100 120 140 160 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 γ = 0 Smoothness penalty 50 100 150 200 20 40 60 80 100 120 140 160 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 γ = 1/2 Smoothness penalty 50 100 150 200 20 40 60 80 100 120 140 160 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 γ = 5 Seismic amplitude recovery Matching procedure Compute reference vector <=> defines g  migrate data  apply spherical-divergence correction Create “data” <=> defines f  demigrate  migrate Estimate scaling by inversion procedure Define scaled curvelet transform Recover migration amplitudes by sparsity promotion. Primary-multiple separation Matching procedure Predict multiples <=> defines g  apply conventional Fourier matched filtering Consider total data as “true” multiples <=> defines f  do not know true multiples  use total data instead  minimize energy mismatch Estimate scaling by an inversion procedure. Define scaled curvelet-domain threshold. Separate primaries & multiples by sparsity promotion. Problem formulation Signal model for total data Multiple prediction by e.g. SRME may contain amplitude errors, i.e., Solve with s the total data. Use z to correct the predicted multiples, i.e., or correct the thresholding s = s1 + s2 s2 = Bs̆2 s2 ≈ CTdiag{w}Cs̆2 s̆2 !→ CTdiag{w̃}Cs̆2 with w̃ = ez̃ t = diag{w̃}|Cs̆2| Jγ(z) = 1 2 ‖s− Fγez‖22, Synthetic example Total data SRME predicted multiples s̆2s Synthetic example s̆1 Curvelet estimated primariesSRME predicted primaries s̃1 = CTTt ( Cp ) t = Cs̆2 Synthetic example Corrected multiplesCorrected multiples s̆corr.2 = C T diag{w}Cs̆2 for γ = 0 s̆corr.2 = CT diag{w}Cs̆2 for γ = 0.5 Synthetic example Scaled thresholded primariesScaled thresholded primaries s̃1 = CTTt ( Cp ) t = diag{w}|Cs̃2| s̃1 = CTTt ( Cp ) t = diag{w}|Cs̃2| Synthetic example Curvelet estimated primaries s̃1 = CTTt ( Cp ) t = Cs̆2 Scaled thresholded primaries s̃1 = CTTt ( Cp ) t = diag{w}|Cs̃2| Real example SRME predicted primariesSRME predicted multiples s̆2 s̆1 Real example Thresholded primaries s̃1 = CTTt ( Cp ) t = Cs̆2 s̃1 = CTTt ( Cp ) t = diag{w}|Cs̃2| Scaled thresholded primaries Conclusions Combining the parsimonious curvelet transform with phase-space structure allows us to control diagonal estimation <=> over fitting handle data with conflicting dips stably recover & separate Application improved migration-amplitude recovery improved primary-multiple separations Future 3-D non-smooth symbols Acknowledgments The authors of CurveLab (Demanet, Ying, Candes, Donoho) Christiaan C. Stolk for his contribution to phase-space smoothness. The SLIM team Sean Ross Ross, Cody Brown and Henryk Modzeleweski for developing SLIMPy: operator overloading in python 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