Stable seismic data recovery Felix J. Herrmann* joint work with Peyman Moghaddam*, Gilles Hennenfent* & Chris Stolk (Universiteit Twente) *Seismic Laboratory for Imaging and Modeling slim.eos.ubc.ca AIP 2007, Vancouver, June 26 Combinations of parsimonious signal representations with nonlinear sparsity promoting programs hold the key to the next-generation of seismic inversion 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 presence of operators (FIO’s & PsDO’s) The seismic method Seismic data acquisition 0 time [s] 1 2 3 4 -3000 offset [m] -2000 -1000 Exploration seismology 0 1 Depth (km) 2 • • • create images of the subsurface 3 4 5 need for higher resolution/deeper 6 clutter and data incompleteness are problems 7 0 1 2 3 km Exploration seismology 0 1 Depth (km) 2 • • • create images of the subsurface 3 4 5 need for higher resolution/deeper 6 clutter and data incompleteness are problems 7 0 1 2 3 km Forward problem F [c]u := • • d 1 ∂ ∂ · − c2 (x) ∂t2 i=1 ∂x21 2 2 second order hyperbolic PDE interested in the singularities of m = c − c¯ u(x, t) = f (x, t) Inverse problem Minimization: m = arg min d − F [m] m 2 2 After linearization (Born app.) forward model with noise: d(xs , xr , t) = Km (xs , xr , t) + n(xs , xr , t) Conventional imaging: K d (x) = T y(x) = K Km (x) + K n (x) T Ψm (x) + e(x) Ψ is prohibitively expensive to invert requires regular sampling ... T Sparsity promoting inversion Formulate as inverse problem signal y = A + n noise x0 curvelet representation of ideal data ˜ = 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. Sparsity promoting inversion x0 can be recovered by solving P : ˜ = arg minx x x ˜f = ST x ˜ 1 s.t. Ax − y with y = (incomplete) data A = modeling matrix, e.g. A = RS ˜ = recovered sparsity vector x T = a number dependent on the noise level ST = the synthesis matrix ˜f = the recovered function f Crux lies in finding the sparse representation! 2 ≤ Curvelets & seismology Wish list Transform that is parsimonious detects the wavefronts localized in space and frequency (phase space) some invariance under “wave propagation” Events correspond to curved singularities with conflicting dips caustics faults & pinch outs Need a transform that is multiscale multidirectional exactly reconstructs 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 k2 2j/2 angular wedge 2j fine scale data k1 coarse scale data 2-D curvelets [Candes, Donoho, Demanet, Ying] 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! Wavefront detection 0 -2000 Offset (m) 0 2000 Time (s) 0.5 1.0 1.5 2.0 Significant curvelet coefficient Curvelet coefficient~0 curvelet coefficient is determined by the dot product of the curvelet function with the data Compression [From Demanet ‘05] 3-D curvelets Curvelets live in wedges in the 3 D Fourier plane... Nonlinear approximation Nonlinear approximation Curvelet-based seismic data recovery joint work with Gilles Hennenfent 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 85 % missing operation by filling in the zero traces. Since seismic the missing data can be Curvelet recovered by compounding recovery modeling operator, i.e., A := RCT . With this definit P corresponds to seeking the sparsest curvelet vec followed by the picking, matches the data at the n transform (with S := C in P ) gives the interpolated An example of curvelet based recovery is presente data volume is recovered from data with 80 % traces traces are selected at random according to a discrete d Observations Inverted a rectangular matrix worked because the curvelet transform is sparse exploits the higher dimensional geometry of seismic wavefields curvelets are incoherent with the Dirac measurement basis Data is recovered for large percentages of traces missing Is an example of an inverse problem with incomplete data Can these ideas be extended to recover migration amplitudes? approximately invert a PsDO diagonalize zero-order PsDO’s Stable seismic amplitude recovery “Sparsity- and continuitypromoting seismic image recovery with curvelet frames” by F.H, P. Moghaddam & C. Stolk to appear in special issue on imaging in ACHA Migrated data Seismic Laboratory for Imaging and Modeling Amplitude-corrected & denoised migrated data 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 need sparsity on the model invariance under the normal operator x p . ≤ Nonlinear approximation Migrated mobil data set Nonlinear approximation Recovery from largest 3 % Nonlinear approximation Difference Diagonal approximation of the Hessian Normal/Gramm operator [Stolk 2002, ten Kroode 1997, de Hoop 2000, 2003] In high-frequency limit Ψ is a PsDO Ψf (x) = • • R −ix·ξ e d a(x, ξ)fˆ(ξ)dξ pseudolocal singularities are preserved Inversion corrects for the ‘Hessian’ Invariance under Gramm matrix (a) (b) (c) (d) (a) (b) (c) (d) curvelets remain invariant (e) (f) (g) (h) (e) (g) (h) • • (f) approximation improves for higher frequencies Approximation stitutions are made for the scattering operator and the leading behavior for their composition, the normal operator Ψ, c leading behavior their composition, the normal operator Ψ, correspond e normal operator Ψ, corresponds to that for of an α f )∧ (ξ) = |ξ|2α · fˆ(ξ). Alter m → order-one (−∆)1/2invertible m withelliptic ((−∆) PsDO . order-one invertible elliptic PsDO . So let Ψ = Ψ(x, D) be a pseudodifferential operator of order 0, with homogeneous principalmade symbolzero-order a(x, ξ). this by composing thetodata side with aby 1/2-o To make PsDO amenable an approximation curv this subPsDO amenable to an approximation by curvelets, the approximation by curvelets, To the make following −1/2 operator and the model: stitutions are for made for the scattering −1/2 stitutions are made the scattering operator thee.g. model: → K erator and the model: K K (−∆) and the→ time coordinate, i.e., K → ∂t Kand (see 3).K After d 2αin ˆ R . · f (ξ). or 1/2 α f )∧ (ξ) 2α · fˆ(ξ). Alternativel 1/2 α ∧ 2α ˆ m → (−∆) m with ((−∆) = |ξ| with m → operators (−∆) mcan with ((−∆) f ) (ξ) = |ξ| · f (ξ). Alternatively, these op Alternatively, these be operator Ψ becomes zero-order. Remark that these subsit made zero-order by composing with a 1/2-order madeintegration zero-order by composing the datathe sidedata withside a 1/2-order fractional fra in de with a 1/2-order fractional along Lemma 1. With C some constant, the following holds tution made in the WVD−1/2 methods, −1/2 where vaguelettes are the time i.e., K → ∂t K → K∂ (see e.g. 3). After these substitutio thecoordinate, time coordinate, i.e., K (see e.g. 3). After these t e e.g. 3). After these substitutions, the normal −|ν|/2 n ≤ C 2 (Ψ(x, D) − a(x , ξ ))ϕ . these subsitutions (14) 2 ν ν ν L ( R ) operator Ψ becomes zero-order. Remark that are simila mappings. Before detailing the approximate diagonalizat operator Ψ becomes zero-order. Remark that these subsitutions hat these subsitutions are similar to the substi- tution made in the WVD methods, where vaguelettes are introduced acco tution made in the WVD methods, where vaguelettes are introd vaguelettes are introduced according the same first discuss the properties of continuous curvelets under t To approximate Ψ, we define the sequence u := (uµ )µ∈M = a(xµ , ξµ ). Let DΨ be the mappings. Before detailingdetailing the approximate diagonalization of the norm mappings. the approximate diagonalization of ate diagonalization of the normal operator, Before we diagonal matrix with entries given by u. Next we state our result on the approximation of first discuss the properties of continuous curvelets under this operator. urvelets under this operator. first discuss the properties of continuous curvelets under this ope T Ψ by C DΨ C. APPROXIMATION OF THE NORM Ψ by C T DΨ C. Approximation Theorem 1. The following estimate for the error holds (Ψ(x, D) − C T DΨ C)ϕµ L2 (R ) n ≤ C 2−|µ|/2 , where C is a constant depending on Ψ. Allows for the decomposition This main result proved in Appendix A shows that the approxima Ψϕµ (x) C DΨ Cϕµ (x) T diagonal approximation goes to zero for increasingly finer scales. The app T = AA ϕµ (x) from the property that the symbol is slowly varying over the support with A := √ DΨ C and AT := C √ T DΨ . approximation that becomes more accurate as the scale increases. Approximation y(x) = Ψm (x) + e(x) AA m (x) + e(x) T = Ax0 + e, Wavelet-vagulette like Amenable to nonlinear recovery Estimation of the diagonal scaling Diagonal estimation • • Define a reference vector (say conventional image). • Define the matrix • Calculate ‘data’ b = Ψr P := C diag(v) T with v = Cr Invert 1 ˜ = arg min b − Pu u u 2 2 2 + η Lu 2 2 2 Diagonal estimation Impose smoothness in phase space extended to include q different reference vectors by making the following substitutions L = [D1 D2 Dθ ] b → [b1 · · · bq ]T and P → [P1 · · · Pq ]T with the “data” vector and “modeling” matrix defined by the different reference vectors r1 , . . . , rq . Calculate: b = Ψr and v = Cr. Set: η = ηmin ; while ∃ (˜ uµ )µ∈M < 0 do Solve u = arg minu 1 2 b − Pu 2 2 + η 2 Lu Increase the Lagrange multiplier λ = η + ∆η end while 2 2 Diagonal estimation (a) (b) Seismic amplitude recovery Recovery Final form y = Ax0 + ε with x0 = ΓCm and Solve P: with = Ae. minx J(x) subject to y − Ax 2 H †˜ ˜ m = (A ) x sparsity J(x) = α x 1 +β Λ 1/2 A H † continuity x p . ≤ Image recovery anisotropic diffusion sotropic-diffusion penalty term (see e.g. 24) is given by [Black et. al ’98, Fehmers et. al. ’03 and Shertzer ‘03] Jc (m) = Λ 1/2 Define ∇m 2 2 T T D2 . discretized gradient matrix defined as ∇ = The block 1/2 Jc (m) = Λ ∇m p s location dependent (see Fig. 10, which plots the gradients) and ro DT1 with p=2 wards the tangents of the reflecting surfaces. This rotation matrix is g +D r 2 1 +D r −D r + υId Λ[r] = 2 1 ∇r 22 + 2υ −D1 r e discretized derivative in the ith coordinate direction and υ a param Gradient of the reference vector 500 1000 depth (m) 1500 2000 2500 3000 3500 2000 4000 6000 8000 lateral (m) 10000 12000 14000 Recovery Step 1: Update of the Jacobian of 1 2 y − Ax 22 : x ← x + AT (y − Ax) ; Step 2: projection onto the 1 ball S = { x 1 ≤ x0 1 } by soft thresholding x ← Tλw (x); Step 3: projection onto the anisotropic diffusion ball C = {x : J(x) ≤ J(x0 )} by x ← x − κ∇x Jc (x) Initialize: m = 0; x0 = 0; y = KT d; Choose: M and L AT y ∞ > λ1 > λ2 > · · · while y − Ax 2 > do m = m + 1; xm = xm−1 ; for l = 1 to L do xm = Tλm xm + AT (y − xm ) {Iterative thresholding} end for Anisotropic descent update; xm = xm − β∇xm Jc (xm ); end while x = xm ; m = AT † x. Table 2: Sparsity-and continuity-enhancing recovery of seismic amplitudes. Application to the SEG AA’ model 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 linearized modeling 64, and migration 294 minutes on 68 CPU’s Scaling required 1 extra migration-demigration Seismic Laboratory for Imaging and Modeling Seismic Laboratory for Imaging and Modeling Seismic Laboratory for Imaging and Modeling Migrated data Seismic Laboratory for Imaging and Modeling Amplitude-corrected & denoised migrated data Noise-free data Noisy data (3 dB) Seismic Laboratory for Imaging and Modeling Data from migrated image Data from amplitude-corrected & denoised migrated image 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 full modeling Scaling required 1 extra migration-demigration Conclusions Curvelet-domain scaling handles conflicting dips (conormality assumption) exploits invariance under the PsDO robust w.r.t. noise Diagonal approximation exploits smoothness of the symbol uses “neighbor” structure of the curvelet transform Results on the SEG AA’ show recovery of amplitudes beneath the Salt successful recovery of clutter improvement of the continuity Acknowledgments 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.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- Stable seismic data recovery
Open Collections
UBC Faculty Research and Publications
Stable seismic data recovery Herrmann, Felix J. 2007
pdf
Page Metadata
Item Metadata
Title | Stable seismic data recovery |
Creator |
Herrmann, Felix J. |
Contributor | University of British Columbia. Seismic Laboratory for Imaging and Modeling |
Date Issued | 2007 |
Description | In this talk, directional frames, known as curvelets, are used to recover seismic data and images from noisy and incomplete data. Sparsity and invariance properties of curvelets are exploited to formulate the recovery by a `1-norm promoting program. It is shown that our data recovery approach is closely linked to the recent theory of “compressive sensing” and can be seen as a first step towards a nonlinear sampling theory for wavefields. The second problem that will be discussed concerns the recovery of the amplitudes of seismic images in clutter. There, the invariance of curvelets is used to approximately invert the Gramm operator of seismic imaging. In the high-frequency limit, this Gramm matrix corresponds to a pseudo-differential operator, which is near diagonal in the curvelet domain. |
Extent | 8790195 bytes |
Subject |
curvelets sparsity Gramm matrix |
Genre |
Presentation |
Type |
Text Still Image |
File Format | application/pdf |
Language | eng |
Date Available | 2008-03-14 |
Provider | Vancouver : University of British Columbia Library |
Rights | All rights reserved |
DOI | 10.14288/1.0107417 |
URI | http://hdl.handle.net/2429/584 |
Affiliation |
Science, Faculty of Earth and Ocean Sciences, Department of |
Citation | Herrmann, Felix J. 2007. Stable seismic data recovery. Conference on Applied Inverse Problems 2007. |
Peer Review Status | Unreviewed |
Scholarly Level | Faculty |
Copyright Holder | Herrmann, Felix J. |
Aggregated Source Repository | DSpace |
Download
- Media
- 52383-herrmann07AIP.pdf [ 8.38MB ]
- Metadata
- JSON: 52383-1.0107417.json
- JSON-LD: 52383-1.0107417-ld.json
- RDF/XML (Pretty): 52383-1.0107417-rdf.xml
- RDF/JSON: 52383-1.0107417-rdf.json
- Turtle: 52383-1.0107417-turtle.txt
- N-Triples: 52383-1.0107417-rdf-ntriples.txt
- Original Record: 52383-1.0107417-source.json
- Full Text
- 52383-1.0107417-fulltext.txt
- Citation
- 52383-1.0107417.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
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.52383.1-0107417/manifest