UBC Faculty Research and Publications

Compressive sampling meets seismic imaging 2008

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

Item Metadata

Download

Media
herrmann07PIMS.pdf
herrmann07PIMS.pdf [ 10.1MB ]
Metadata
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
Citation
1.0102592.ris

Full Text

Compressive sampling meets seismic imaging Felix J. Herrmann joint work with Tim Lin*, Peyman Moghaddam*, Gilles Hennenfent*, Deli Wang* & Chris Stolk (Universiteit Twente) *Seismic Laboratory for Imaging and Modeling slim.eos.ubc.ca PIMS-CINVESTAV 2007, Monterrey, October 19 Today’s challenges Aside 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 operators Today’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 Motivation Exploit two aspects of curvelets, namely their  parsimoniousness  invariance under certain operators Formulate  data-adaptive scaling algorithms  non-adaptive wavefield reconstruction algorithms Applications  nonlinear migration-amplitude recovery  nonlinear sampling for wavefields  nonlinear sampling for operators Today’s topics Sparsity-promoting seismic-image amplitude recovery  curvelet-domain diagonal approximation of PsDO’s  stable sparsity-promoting inversion Directional frame-based wavefield reconstruction by sparsity promotion  curvelet parsimoniousness  jitter sampling Compression of FIO’s through compressive sampling  measurement basis diagonalizes operator The problem Minimization: After linearization (Born app.) forward model with noise: Conventional imaging:     is prohibitively expensive to invert evaluation of         involves expensive wavefield extrapolators ( KT d ) (x) = ( KTKm ) (x) + ( KTn ) (x) y(x) = ( Ψm ) (x) + e(x) Ψ d(xs, xr, t) = ( K[c̄]m ) (xs, xr, t) + n(xs, xr, t) K[c̄] c̃ = arg min c ‖d− F [c]‖22 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 Curvelets are oscillatory in one direction and smooth in the others. 3-D curvelets Approximate linearized inversion by curvelet scaling & sparsity promotion Joint work with Chris Stolk* and Peyman Moghaddam Mathematics Department, Twente University, the Netherlands “Sparsity- and continuity-promoting seismic imaging with  curvelet frames” to appear in ACHA Related work Wavelet-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 models Corresponds to a spatially-varying dip filter after appropriate preconditioning (=> zero-order PsDO). Ψ( Ψf ) (x) := ( KT Kf ) (x) = ∫ Rd e−ix·ξa(x, ξ)f̂(ξ)dξ c̄ • Substitutions: So let Ψ = Ψ(x,D) be a pseudodifferential operator of order 0, with homo- geneous principal symbol a(x, ξ). in Rd. Lemma 1. With C ′ some constant, the following holds ‖(Ψ(x,D)− a(xν , ξν))ϕν‖L2(Rn) ≤ C ′2−|ν|/2. (14) To approximate Ψ, we define the sequence u := (uµ)µ∈M = a(xµ, ξµ). Let DΨ be the diagonal matrix with entries given by u. Next we state our result on the approximation of Ψ by CTDΨC. Theorem 1. The following estimate for the error holds ‖(Ψ(x,D)− CTDΨC)ϕµ‖L2(Rn) ≤ C ′′2−|µ|/2, (15) where C ′′ is a constant depending on Ψ. This main result proved in Appendix A shows that the approximation error for the diagonal approximation goes to zero for increasingly finer scales. The approximation derives from the property that the symbol is slowly varying over the support of a curvelet, an approximation that becomes more accurate as the scale increases. Decomposition of the normal operator By virtue of Theorem 1, the normal operator can be factorized ( Ψϕµ ) (x) $ (CTDΨCϕµ)(x) (16) = ( AATϕµ ) (x) with A := √ DΨC and AT := CT √ DΨ. Because the seismic reflectivity can be written as a superposition of curvelets, we can replace ϕµ in the above equation with the model m. We 15 leading behavior for their composition, the normal operator Ψ, corresponds to that of an order-one invertible elliptic PsDO . To make this PsDOamenable to an approximation by curvelets, the following sub- stitutions are made for the scattering operator and the model: K !→ K (−∆)−1/2 and m !→ (−∆)1/2 m with ((−∆)αf)∧(ξ) = |ξ|2α · f̂(ξ). Alternatively, these operators can be made zero-order by composing the data side with a 1/2-order fractional integration along the time coordinate, i.e., K !→ ∂−1/2t K (see e.g. 3). After these substitutions, the normal operator Ψ becomes zero-order. Remark that these subsitutions are similar to the substi- tution made in the WVD methods, where vaguelettes are introduced according the same mappings. Before detailing the approximate diagonalization of the normal operator, we first discuss the properties of continuous curvelets under this operator. APPROXIMATION OF THE NORMAL OPERATOR In this section, a diagonal approximation of the normal operator in the curvelet domain is presented. Invariance properties of curvelets under the normal operator (see also Fig. 2) are used. The approximation leads to a SVD-like decomposition of the normal operator and makes large-scale seismic image recovery amenable to optimization. To understand our approximation, we first list the important properties of continuous curvelets. An upper bound for the L2-error of the diagonal approximation is discussed next, followed by the diagonal decomposition of the normal operator and a method to numerically estimate the diagonal from discrete implementations of the normal operator. We conclude this section by discussing the empirical performance of the approximation on a synthetic data set. 11 l adi g behavior for their composition, the normal operator Ψ, corresponds to that of an order-one invertible elliptic PsDO . To make this PsDOamenable to an approximation by curvelets, the following sub- stitutions are made for the scattering operator and the model: K !→ K (−∆)−1/2 and m !→ (−∆)1/2 m with ((−∆)αf)∧(ξ) = |ξ|2α · f̂(ξ). Alternatively, these operators can be ade zero-order by c mposing the data side with a 1/2-order fractional integration along the time coordinate, i.e., K !→ ∂−1/2t K (see e.g. 3). After these substitutions, the normal operator Ψ becomes zero-order. Remark that these subsitutions are similar to the substi- tution made in the WVD methods, where vaguelettes are introduced according the same mappings. Before detailing the approximate diagonalization of the normal operator, we first discuss the properties of continuous curvelets under this operator. APPROXIMATION OF THE NORMAL OPERATOR In this section, a diagonal approximation of the normal operator in the curvelet domain is presented. Invariance properties of curvelets under the normal operator (see also Fig. 2) are used. The approximation leads to a SVD-like decomposition of the normal operator and makes large-scale seismic image recovery amenable to optimization. To understand our approximation, we first list the important properties of continuous curvelets. An upper bound for the L2-error of the diagonal approximation is discussed next, followed by the diagonal decomposition of the normal operator and a method to numerically estimate the diagonal from discrete implementations of the normal operator. We conclude this section by discussing the empirical performance of the approximation on a synthetic data set. 11 leading behavior for their composition, the normal operator Ψ, corresponds to that of an order-one invertible elliptic PsDO . To make this PsDOamenable to an approximation by curvelets, the following sub- stitutions are made for the scattering operator and the model: K !→ K (−∆)−1/2 and m !→ (−∆)1/2 m with ((−∆)αf)∧(ξ) = |ξ|2α · f̂(ξ). Alternatively, these operators can be made z o-order by c mposing the data side with a 1/2-order fractional integration along the time coordinate, i.e., K !→ ∂−1/2t K (see e.g. 3). After these substitutions, the normal operator Ψ becomes zero-order. Remark that these subsitutions are similar to the substi- tution made in the WVD methods, where vaguelettes are introduced according the same map ings. Before detailing the approximate diagonalization of the normal operator, we first discuss the properties of continuous curvelets under this operator. APPROXIMATION OF THE NORMAL OPERATOR In this section, a diagonal approximation of the normal operator in the curvelet domain is presented. Invariance properties of curvelets under the normal operator (see also Fig. 2) are used. The approximation leads to a SVD-like decomposition of the normal operator and makes large-scale seismic image recovery amenable to optimization. To understand our approximation, we first list the important properties of continuous curvelets. An upper bound for the L2-error of the diagonal approximation is discussed next, followed by the diagonal decomposition of the normal operator and a method to numerically estimate the diagonal from discrete implementations of the normal operator. We conclude this section by discussing the empirical performance of the approximation on a synthetic data set. 11 leading behavior for their composition, the normal operator Ψ, corresponds to that of an order-one invertible elliptic PsDO . To make this PsDOamenable to an approximation by curvelets, the following sub- stitutions are made for the scattering operator and the model: K !→ K (−∆)−1/2 and m !→ (−∆)1/2 m with ((−∆)αf)∧(ξ) = |ξ|2α · f̂(ξ). Alternatively, these operators can be made zero-order by composing the d ta side wi h a 1/2-order fractional integration along the time coordinate, i.e., K !→ ∂−1/2t K (see e.g. 3). After these substitutions, the normal operator Ψ becomes zero-order. Remark that these subsitutions are similar to the substi- tution made in the WVD methods, where vaguelettes are introduced according the same mappings. Befor d tailing the approximate diagonalization of the normal operator, we first discuss the properties of continuous curvelets under this operator. APPROXIMATION OF THE NORMAL OPERATOR In this section, a diagonal approximation of the normal operator in the curvelet domain is presented. Invariance properties of curvelets under the normal operator (see also Fig. 2) a sed. The approximation leads to a SVD-like decomposition of the normal operator a d makes large-scale seismic image recovery amenable to optimization. To understand our approximation, we first list e important roperties of continuous curvelets. An upper bound for the L2-error of the diagonal approximation is discussed next, followed by the diagonal decomposition of the normal operator and a method to numerically estimate the diagonal from discrete implementations of the normal operator. We conclude this section by discussing the empirical performance of the approximation on a synthetic data set. 11 or with App oximation Tiling the ξ space ~ 2 ~ 2 j j/2 µ! In red, the essential frequency support of a curvelet φµ. 10 Allows for decomposition of the normal operator in Rd. Lemma 1. With C ′ some constant, the following holds ‖(Ψ(x,D)− a(xν , ξν))ϕν‖L2(Rn) ≤ C ′2−|ν|/2. (14) To approximate Ψ, we define the sequence u := (uµ)µ∈M = a(xµ, ξµ). Let DΨ be the diagonal matrix with entries given by u. Next we state our result on the approximation of Ψ by CTDΨC. Theorem 1. The following estimate for the error holds ‖(Ψ(x,D)− CTDΨC)ϕµ‖L2(Rn) ≤ C ′′2−|µ|/2, (15) where C ′′ is a constant depending on Ψ. This main result proved in Appendix A shows that the approximation error for the diagonal approximation goes to zero for increasingly finer scales. The approximation derives from the property that the symbol is slowly varying over the support of a curvelet, an approximation that becomes more accurate as the scale increases. Decomposition of the normal operator By virtue of Theorem 1, the normal operator can be factorized ( Ψϕµ ) (x) $ (CTDΨCϕµ)(x) (16) = ( AATϕµ ) (x) with A := √ DΨC and AT := CT √ DΨ. Because the seismic reflectivity can be written as a superposition of curvelets, we can replace ϕµ in the above equation with the model m. We 15 ( Ψϕµ ) (x) ! (CTDΨCϕµ)(x) = ( AATϕµ ) (x) with A := √ DΨC and AT := CT √ DΨ. Scaling 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. Key idea Estimation curvelet-domain scaling  inversion of an underdetermined system  over fitting  positivity and reasonable scaling Solution:  use smoothness of the symbol  formulate nonlinear estimation problem that minimizes with  solve with l-BFGS [Noccedal, Symes ‘07] 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 Our approach “Forward” model:  diagonal approximation of the demigration-migration operator  costs one demigration-migration to estimate the diagonal weighting with y = migrated data A := CTΓ AAT r ≈ KTKr K = the demigration operator ! = migrated noise. y = KTKm + ε ≈ Ax0 + ε Solution Solve P :  minx J(x) subject to ‖y −Ax‖2 ≤ ! m̃ = (AH)†x̃ with J(x) = sparsity︷ ︸︸ ︷ α‖x‖1 +β ‖Λ1/2 ( AH )† x‖p︸ ︷︷ ︸ continuity .    Migrated data Amplitude-corrected & denoised migrated data Imaging 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’sKT Observations Curvelet-domain scaling  handles conflicting dips (conormality assumption)  exploits invariance under the PsDO Diagonal approximation  exploits smoothness of the symbol  uses “neighbor” structure of curvelets Results on the SEG AA’ show  recovery of amplitudes beneath the Salt  successful recovery from clutter  improvement of the continuity  robust w.r.t. noise Curvelet-domain matched filter ... A primer on compressive sampling Compressive 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 cooling Seismic Laboratory for Imaging and Modeling Simple example Fourier transform ✓ ✗ 3-fold under-sampling significant coefficients detected ambiguity few significant coefficients Fourier transform Fourier transform Seismic Laboratory for Imaging and Modeling Forward problem x0 A A := RFH y = Fourier coefficients (sparse) with Fourier transform restriction operator signal Seismic Laboratory for Imaging and Modeling Naive sparsity-promoting recovery inverse Fourier transform detection + data-consistent amplitude recovery Fourier transform y A H = A y = detection Ar data-consistent amplitude recovery y A † r = x0 Seismic Laboratory for Imaging and Modeling  “noise” – due to AHA ≠ I – defined by AHAx0-αx0 = AHy-αx0 Undersampling “noise” less acquired data 3 detectable Fourier modes 2 detectable Fourier modes 1 out of 2 1 out of 4 1 out of 6 1 out of 8 Seismic Laboratory for Imaging and Modeling Sparsity-promoting wavefield reconstruction x0 Ay = with sparsifying transform for seismic data restriction operator A := RS H [Sacchi et al ‘98] [Xu et al ‘05] [Zwartjes and Sacchi ‘07] [Herrmann and Hennenfent ‘07] complete wavefield  (transform domain) acquired data Interpolated data given by                 withf̃ = SH x̃ x̃ = arg min x ||x||1 s.t. y = Ax Seismic Laboratory for Imaging and Modeling Observations  bla bla  generalized to A=RMS^H  depends on solver, sampling strategy and sparsity transform Compressive sampling of wavefields joint 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 of recovers the function f. General form compressive sampling P! : { x̃ = arg minx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ = ST x̃ with A = RMST R = restriction matrix M = measurement matrix ST = sparsity synthesis matrix y = RMf The problem Total data 85 % traces missing Requirements Sparsifying transform (S)  curvelet  focussed curvelets Sampling scheme (RM)  random sampling  random jittered sampling => control largest gaps Sparsity promoting solver (P)  Iterative thresholding (Landweber + soft threshold) Seismic Laboratory for Imaging and Modeling Discrete random jittered undersampling receiver positions receiver positions PDF receiver positions PDF receiver positions PDF [Hennenfent and Herrmann ‘07] Typical spatial convolution kernel (amplitudes) Averaged spatial convolution kernel (amplitudes) Sampling schemeType po or ly jit ter ed op tim all y jit ter ed ra nd om re gu lar Solution of recovers the wavefield f. Curvelet-based recovery P! : { x̃ = arg minx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ = ST x̃ with A = RICT R = jitter sampling I = Dirac basis CT = curvelet synthesis y = Rf Seismic Laboratory for Imaging and Modeling Model Seismic Laboratory for Imaging and Modeling Regular 3-fold undersampling Seismic Laboratory for Imaging and Modeling CRSI from regular 3-fold undersampling SNR = 20× log10 ( ‖model‖2 ‖reconstruction error‖2 ) SNR = 6.92 dB Seismic Laboratory for Imaging and Modeling Random 3-fold undersampling Seismic Laboratory for Imaging and Modeling CRSI from random 3-fold undersampling SNR = 20× log10 ( ‖model‖2 ‖reconstruction error‖2 ) SNR = 9.72 dB Seismic Laboratory for Imaging and Modeling Optimally-jittered 3-fold undersampling Seismic Laboratory for Imaging and Modeling CRSI from opt.-jittered 3-fold undersampling SNR = 10.42 dB Seismic Laboratory for Imaging and Modeling Model Seismic Laboratory for Imaging and Modeling Regular 3-fold undersampling SNR = 12.98 dB Seismic Laboratory for Imaging and Modeling SNR = 12.98 dB Regular 3-fold undersampling Seismic Laboratory for Imaging and Modeling Optimally-jittered 3-fold undersampling SNR = 15.22 dB Seismic Laboratory for Imaging and Modeling Optimally-jittered 3-fold undersampling SNR = 15.22 dB Solution of recovers the wavefield f. Focussed recovery with A = R∆PCT ∆P = main primaries y = Rf P! : { x̃ = arg minx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ =∆PCT x̃ 80 % missing Focused curvelet recovery Curvelet-based processing 3 SPARSITY-PROMOTING INVERSION Our solution strategy is built on the premise that seismic data and images have a sparse representation, x0, in the curvelet domain. To exploit this property, our forward model reads y = Ax0 + n (1) with y a vector with noisy and possibly incomplete mea- surements; A the modeling matrix that includes CT ; and n, a zero-centered white Gaussian noise. Because of the redundancy of C and/or the incompleteness of the data, the matrix A can not readily be inverted. However, as long as the data, y, permits a sparse vector, x0, the ma- trix, A, can be inverted by a sparsity-promoting program (Candès et al., 2006b; Donoho, 2006) of the following type: P! : { x̃ = argminx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ = ST x̃ (2) in which ! is a noise-dependent tolerance level, ST the inverse transform and f̃ the solution calculated from the vector x̃ (the symbol ˜ denotes a vector obtained by non- linear optimization) that minimizes P!. Nonlinear programs such as P! are not new to seismic data processing and imaging. Refer, for instance, to the extensive literature on spiky deconvolution (Taylor et al., 1979) and transform-based interpolation techniques such as Fourier-based reconstruction (Sacchi and Ulrych, 1996). By virtue of curvelets’ high compression rates, the non- linear program P! can be expected to perform well when CT is included in the modeling operator. Despite its large- scale and nonlinearity, the solution of the convex problem P! can effectively be approximated with a limited (< 250) number of iterations of a threshold-based cooling method derived from work by Figueiredo and Nowak (2003) and Elad et al. (2005). Each step involves a descent projection, followed by a soft thresholding. SEISMIC DATA RECOVERY The reconstruction of seismic wavefields from regularly- sampled data with missing traces is a setting where a curvelet-based method will perform well (see e.g. Herr- mann, 2005; Hennenfent and Herrmann, 2006a, 2007). As with other transform-based methods, sparsity is used to reconstruct the wavefield by solving P!. It is also shown that the recovery performance can be increased when in- formation on the major primary arrivals is included in the modeling operator. Curvelet-based recovery The reconstruction of seismic wavefields from incomplete data corresponds to the inversion of the picking operator R. This operator models missing data by inserting zero traces at source-receiver locations where the data is miss- ing. The task of the recovery is to undo this operation by filling in the zero traces. Since seismic data is sparse in the curvelet domain, the missing data can be recovered by compounding the picking operator with the curvelet modeling operator, i.e., A := RCT . With this defini- tion for the modeling operator, solving P! corresponds to seeking the sparsest curvelet vector whose inverse curvelet transform, followed by the picking, matches the data at the nonzero traces. Applying the inverse transform (with S := C in P!) gives the interpolated data. An example of curvelet based recovery is presented in Figure 1, where a real 3-D seismic data volume is recov- ered from data with 80% traces missing (see Figure 1(b)). The missing traces are selected at random according to a discrete distribution, which favors recovery (see e.g. Hen- nenfent and Herrmann, 2007), and corresponds to an av- erage sampling interval of 125m . Comparing the ’ground truth’ in Figure 1(a) with the recovered data in Figure 1(c) shows a successful recovery in case the high-frequencies are removed (compare the time slices in Figure 1(a) and 1(c)). Aside from sparsity in the curvelet domain, no prior information was used during the recovery, which is quite remarkable. Part of the explanation lies in the curvelet’s ability to locally exploit the 3-D structure of the data and this suggests why curvelets are successful for complex datasets where other methods may fail. Focused recovery In practice, additional information on the to-be-recovered wavefield is often available. For instance, one may have access to the predominant primary arrivals or to the ve- locity model. In that case, the recently introduced focal transform (Berkhout and Verschuur, 2006), which ’decon- volves’ the data with the primaries, incorporates this addi- tional information into the recovery process. Application of this primary operator,∆P, adds a wavefield interaction with the surface, mapping primaries to first-order surface- related multiples (see e.g. Verschuur and Berkhout, 1997; Herrmann, 2007). Inversion of this operator, strips the data off one interaction with the surface, focusing pri- maries to (directional) sources, which leads to a sparser curvelet representation. By compounding the non-adaptive curvelet transform with the data-adaptive focal transform, i.e., A := R∆PCT , the recovery can be improved by solving P!. The solution of P! now entails the inversion of ∆P, yielding the spars- est set of curvelet coefficients that matches the incomplete data when ’convolved’ with the primaries. Applying the inverse curvelet transform, followed by ’convolution’ with ∆P yields the interpolation, i.e. ST :=∆PCT. Compar- ing the curvelet recovery with the focused curvelet recov- ery (Fig ?? and ??) shows an overall improvement in the recovered details. SEISMIC SIGNAL SEPARATION Predictive multiple suppression involves two steps, namely multiple prediction and the primary-multiple separation. In practice, the second step appears difficult and adap- Curvelet recovery SEISMIC DATA RECOVERY The reconstruction of seismic wavefields from regularly-sampled data with missing traces is a setting where a curvelet-based method will perform well (see e.g. Herrmann, 2005; Hennenfent and Herrmann, 2006a, 2007). As with other transform-based methods, sparsity is used to reconstruct the wavefield by solving P!. It is also shown that the recovery performance can be increased when information on the major primary arrivals is included in the modeling operator. Curvelet-based recovery The reconstruction of seismic wavefields from incomplete data corresponds to the inversion of the picking operator R. This operator models missing data by inserting zero traces at source-receiver locations where the data is missing. The task of the recovery is to undo this operation by filling in the zero traces. Since seismic data is sparse in the curvelet domain, the missing data can be recovered by compounding the picking operator with the curvelet modeling operator, i.e., A := RCT . With this definition for the modeling operator, solving P! corresponds to seeking the sparsest curvelet vector whose inverse curvelet transform, followed by the picking, matches the data at the nonzero traces. Applying the inverse transform (with S := C in P!) gives the interpolated data. An example of curvelet based recovery is presented in Figure 1, where a real 3-D seismic data volume is recovered from data with 80% traces missing (see Figure 1(b)). The missing traces are selected at random according to a discrete distribution, which favors recovery (see e.g. Hennenfent and Herrmann, 2007), and corresponds to an average sampling interval of 125m . Comparing the ’ground truth’ in Figure 1(a) with the recovered data in Figure 1(c) 5 Original data Observations Regular subsampling is unfavorable  random sampling favorable but suffers from gaps  jitter sampling favorable and controls gaps Focal transform  is reminiscent of an imaging operator  improves recovery <=> additional compression Solver  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 extrapolation Compressed wavefield extrapolation joint work with Tim Lin “Compressed wavefield extrapolation” in Geophysics Motivation Synthesis of the discretized operators form bottle neck of imaging Operators have to be applied to multiple right-hand sides Explicit operators are feasible in 2-D and lead to an order-of-magnitude performance increase Extension towards 3-D problematic  storage of the explicit operators  convergence of implicit time-harmonic approaches First go at the problem using CS techniques to compress the operator ... Related work Curvelet-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 domain signal in space domain F R L1 L1 incomplete signal in Fourier domain incomplete and shifted signal in Fourier domain shifted signal in space domain signal in space domain Compressed Sensing Compressed Processing e−j ∆x 2pi Λe j x 2piRF Suppose we want to shift a sparse spike train, i.e.,  Eigen modes <=> Fourier transform.  Can this operation be compressed by compressive sampling? Inspiration u = Tτv = e−τDv = Le−jτΩLHv where D = LΩLH L = The Fourier Transform Calculate instead  Take compressed measurements in Fourier space.  Recover with sparsity promotion  Shift operator is compressed by the restriction yielding compressed rectangular operators.  Extend this idea to wavefield extrapolation? Operators on spikes [Candes et. al, Donoho]  y′ = RejΩτFv A = RF ũ = arg minu ‖u‖1 s.t. Au = y′ R ∈ Rm×N with m" N Representation for seismic data [Berkhout] s+ W- R+ W+ p- Different representations diagonalization operator parsimony wavefield SVD/Lanczos/ modal ✓ ✕ curvelets ✕ ✓ Different representations diagonalization operator parsimony wavefield SVD/Lanczos/ modal ✓ ✕ curvelets ✕ ✓ If incoherent this may actually work .... Sparsity promoting formulation Buys us stability w.r.t. missing data  provided measurement and sparsity representations are mutually incoherent  sufficient mixing <=> random restriction Different strategy:  Let the physics define the measurement basis  Use the modal domain (domain of eigenfunctions) to define the measurement basis  See what you can recover Study eigenfunctions:  mutual coherence with sparsity representation  modal spectrum on the to-be-extrapolated wavefield One-Way Wave Operator  Structure of     confounds the meaning of its exponentiation, due to it being an operator     contains information about medium velocity A = ( 0 ωρ1ωρ1/2 (H2ρ−1/2) 0 ) Two-way Wave Operator A H2 = k2(x,ω) + ∂µ∂µ (Simon & Reed; Dessing ‘97; Grimbergen ‘98) H2 One-Way Wave Operator  Solution of the one-way wave equation  After discretization solve eigenproblem on  Helmholtz operator is Hermitian  monochromatic  velocity   varies laterally H2 (Claerbout, 1971; Wapenaar and Berkhout, 1989) W(x3;x′3) = exp(−j(x3 − x′3)H1) H2 =  ( ω c̄1 )2 0 · · · 0 0 ( ω c̄2 )2 · · · 0 ... ... . . . ... 0 0 · · · ( ω c̄n1 )2 +D2 c̄ 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 diagonalized H1 = LΛ1/2LH L W = FHLe−jΛ1/2(x3−x′3)LHF Compressed wavefield extrapolation Recorded DataOriginal events Forward model Reconstruct point scatterers  from recorded data .... u = Le−jΛ 1/2 ∆x3LHv 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 recovery  y = RLHu A = RejΛ 1/2 ∆x3LH x̃ = arg minx ‖x‖1 s.t. Ax = y ṽ = x̃ Compressed wavefield extrapolation Recorded Data Reconstructed events Reconstruction Only 1 % of original modes were used ... (a) (b) (c) (d) (e) (f) (g) (h) Figure 2: Illustration of the dip limitation during inverse extrapolation. (a) a bandwidth- limited impulsive source. (b) the forward extrapolated wavefield (cf. Eq. (24)). (c) The refocused pulse through inverse extrapolation with matched filtering (cf. Eq. (31)). (d) the same but with least-squares inverse extrapolation (cf. Eq. (32)); (e-h) The time-spatial spectra of (a-d). Notice the lack of spatial frequencies corresponding to steep dips in (f-h). 56  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 ... Observations Compressed wavefield extrapolation (a) (b) Figure 10: Initial wavefields used for the stylized extrapolation examples. (a) A chain of horizonally-oriented fine-scale curvelets playing the role of a “plane wave”. (b) A fan of fine-scale curvelets with different angle. Herrmann et.al. – 64  Extend to general wavefields  Use curvelets as the sparsity representation  Use the full & compressed forward operator operator  Compressively extrapolate back 600m to the source (a) (b) Figure 8: Lateral velocity profiles with background velocity 2000ms−1. (a) Profile with velocity low of 1200ms−1. (b) Profile with velocity high of 3500ms−1. Spatial sampling interval of the profiles is set to 4m with 256 samples, while the sigma of both Gaussian functions are set to 80m. Herrmann et.al. – 62 Restriction & sparsity strategies  Forward extrapolation:  Inverse extrapolation: W1 :  y′ = RejΛ 1/2 ∆x3LH A := RLHFCT x̃ = arg minx ‖x‖1 s.t. Ax = y′ ũ = CT x̃, F1 :  y = RLHFu A′ = RejΛ 1/2 ∆x3LHCH x̃ = arg minx ‖x‖1 s.t. A′x = y ṽ = CT x̃. Forward Extrapolation (a) (b) (c) (d) Figure 11: Compressed forward extrapolation according to W1 (cf. Eq. (42)) for different restrictions. The velocity model corresponds to the velocity low and is plotted in Fig. 8(a). The initial source wavefield v is plotted in Fig. 10(a). (a) The full extrapolated wavefield u = Wv is included for reference; (b) The compressed forward propagated wavefield with pf = 0.2 and pµ = 0.0.2; (c) The same as (b) but with pf = 0.4 and pµ = 0.4; (d) The same as (b) but with pf = 0.6 and pµ = 0.4. Observe that the forward propagated wavefield is largely recovered for the restriction in (c). Herrmann et.al. – 65  (a) is Full extrapolation  (b)-(d) is compressed extrapolation, (b) p = 0.04, (c) p = 0.16, (d) p = 0.24 Inverse Extrapolation (a) (b) (c) (d) (e) (f) Figure 13: Compressed inverse extrapolation according F1 (cf. Eq. (45)) for different re- strictions. For (a-c) the velocity model is given by the Gaussian low (Fig. 8(a)) and in (d-f) by the Gaussian high (Fig. 8(b)). The initial source wavefield v is plotted in Fig. 10(a). (a) Inverse extrapolated wavefield for pf = 0.2 and pµ = 0.2; (b) The same as (a) but with pf = 0.4 and pµ = 0.4;(c) The same as (a) but with pf = 0.6 and pµ = 0.4; (d-f) the same as (a-c) but for the velocity high. Observe that the recovery for the velocity high is slightly better. Herrmann et.al. – 67  (a) p = 0.04  (b) p = 0.16, (c) pf=0.4, px=0.4 Evanescent Recovery (a) (b) (c) Figure 14: Inversion of the evanescent wavemodes according ṽm = WTu or ṽ = F1[u] (cf. Eq. 45). The velocity model is constant at 2000ms−1. The initial source wavefield, v, is defined in terms of a the curvelet fan plotted in Fig. 10(b). (a) The full forward propgated wavefield u = Wv; (b) The matched filter; (c) The !1 recovery. Observe that the steep evanescent angles are fully recovered. Herrmann et.al. – 68  (a) is downward extrapolated wavefield  (b) is matched filter  (c) is “compressed” inverse extrapolation Velocity model Figure 13: Lateral velocity profile for the overthrust examples.Herrmann et.al. – 68 Compressed inverse extrapolation Overthrust exploding reflector Full forward extrapolation Matched filter Recovered from p=0.25 Multiscale and angular compressed wavefield extrapolation  Propose a scheme motivated by extensions of CS  adapt discretization & restriction  parallel implementation Fj1 :  yj = RjMju A′j := RjM ′ jC T j x̃j = arg minxj ‖xj‖1 s.t. A′jxj = yj ṽ = ∑ j C T j x̃j, with j = {j, l} the scale and angle. (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 ... Acknowledgments The audience for listening and the organizers for putting this great workshop together .... The authors of CurveLab (Demanet, Ying, Candes, Donoho) Dr.  W. W. Symes for 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), from the following organizations: BG Group, BP, Chevron,ExxonMobil and Shell.

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 4 0
China 3 0
Germany 1 1
Japan 1 0
City Views Downloads
Beijing 3 0
Sunnyvale 2 0
Unknown 1 1
Tokyo 1 0
Redmond 1 0
Ashburn 1 0

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

Share

Share to:

Comment

Related Items