UBC Faculty Research and Publications

Robust seismic amplitude recovery using curvelets Moghaddam, Peyman P. 2008

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

Item Metadata


moghaddam07seg.pdf [ 166.73kB ]
JSON: 1.0107426.json
JSON-LD: 1.0107426+ld.json
RDF/XML (Pretty): 1.0107426.xml
RDF/JSON: 1.0107426+rdf.json
Turtle: 1.0107426+rdf-turtle.txt
N-Triples: 1.0107426+rdf-ntriples.txt

Full Text

Robust seismic amplitude recovery using curvelets Peyman P. Moghaddam∗, Felix Herrmann, University of British Columbia, Canada and Chris Stolk, Uni- versity of Twente, Netherlands SUMMARY In this paper, we recover the amplitude of a seis- mic image by approximating the normal (demigration- migration) operator. In this approximation, we make use of the property that curvelets remain invariant under the action of the normal operator. We propose a seis- mic amplitude recovery method that employs an eigen- value like decomposition for the normal operator using curvelets as eigen-vectors. Subsequently, we propose an approximate non-linear singularity-preserving solu- tion to the least-squares seismic imaging problem with sparseness in the curvelet domain and spatial continu- ity constraints. Our method is tested with a reverse-time ’wave-equation’ migration code simulating the acoustic wave equation on the SEG-AA salt model. INTRODUCTION In the mid-90s, Hart Smith (Smith, 1997, 1998) intro- duced waveforms that are invariant under the action of certain Fourier and pseudo-differential operators. Based on his theory, E. Candes (2000) designed and imple- mented curvelets which behave following similar ideas. Since then it has been used in different area of signal and image processing. Candes and Demanet (2002) showed how curvelets behave under the action of Fourier inte- gral operators. They stated that under certain conditions a curvelet maps to another curvelet under the action of Fourier Integral Operator (FIO). Recently efforts have been made to exploit this theoret- ical invariance property towards a diagonalization of mi- gration operators (see e.g. Chauris, 2006). Since curvelets are discrete and hence move around on grids, this makes it a challenge to define fast migration operators. Curvelets, however, prove to be very useful for solving the seismic amplitude recovery problem, during which curvelets are being imaged. This seismic amplitude recovery problem involves the inversion of the normal (migration-demigration) operator whose behavior corresponds to that of a non- stationary convolution (i.e., pseudodifferential operator). Since curvelets in this case no longer move, their rela- tive invariance can be used to diagonally approximate and subsequently approximately invert the normal oper- ator, an approach followed in this abstract. During the inversion curvelet sparsity and continuity along imaged reflectors are jointly promoted leading to a stable recov- ery of the seismic amplitudes in particular along steeply dipping events and under the salt. THE PROBLEM FORMULATION After linearization and by ignoring the source and re- ceiver signatures, the discretized forward model that gen- erates seismic data can be written as d= Km+n, (1) with d the first-order Born approximation data, K the scattering operator (or commonly called de-migration or modeling operator), n is the noise and m represents the (singular) fluctuations in the earth’s acoustic properties with respect to an appropriately chosen smoothly vary- ing background velocity model (the density of mass is assumed constant). These fluctuations are referred to as the model and seismic imaging aims to recover both the locations and the relative amplitudes of the velocity fluc- tuations from seismic data. Applying the adjoint of the linearized scattering operator K to the data vector ( d in equation 1) leads to the migrated image. KTd = KTKm+KTn (2) y = ( Ψm ) + e, withΨ =KTK defined as normal operator and e=KTn is the colored noise. An extensive literature has emerged on restoring the mi- gration amplitude by inverting the normal matrix (Ψ = KTK) (Nemeth. et al., 1999; Kuhl and Sacchi, 2003; Herrmann et al., 2006), and involves the computation of the pseudo inverse of the scattering or demigration ma- trix. Unfortunately, the normal operator it too big to be constructed explicitly and is too expensive to be eval- uated as part of an iterative Krylov-subspace solver to invert the Hessian (see e.g. Nemeth. et al., 1999). THEORY: DIAGONALAPPROXIMATIONOFNOR- MAL OPERATOR A curvelet φµ is defined by its index µ which is a triple ( j,k, l)with j= 0,1,2, ... a scale index; l= 0,1, ...,2b j/2c− 1 an orientation index (bxc is the integer part of x ); and Robust Seismic Amplitude Recovery Using Curvelets k = (k1,k2) is the location index. One of the impor- tant features of curvelet is that the action of a pseudo- differential operator (i.e., Ψ ) on a curvelet φµ is ap- proximately mapped into the same curvelet . To be more specific, a pseudo-differential operator induces a map- ping µ 7→ µ ′ with property that significant curvelet co- efficients of Ψ(φµ) are located very close to µ , itself. We propose following decomposition of the normal op- erator Ψ ≈CTΛC, (3) where Λ is a diagonal matrix, C and CT are the curvelet and the transpose of curvelet transform. This decom- position states that we can decompose the normal oper- ator in a form of eigenvalue-like decomposition where curvelets play the role of eigenfunctions and the diago- nal Λ the role of the eigenvalues. Since the eigenvalues are not known, we propose a method that estimates this diagonal from the action of the normal operator on a reference vector (r) that is close enough to the original model. Ψr=CTΛCr 7→CTdiag(u)λ =Ψr, (4) with u=Cr and λ is the diagonal elements of Λ. During this diagonal estimation method, we aim to find a curvelet-domain diagonal weighting vector (i.e., λ ) that approximates the action of the normal operator on the reference vector. We choose the spherical-spreading corrected migrated image as the reference vector. Be- cause the curvelet transform is redundant, the estima- tion of this diagonal in equation refeqn:refvec is an un- derdetermined inverse problem. To find an unique solu- tion to this diagonal approximation problem we impose additional smoothing constrains on the to-be-estimated diagonal (λ ). These constraints reflect the behavior of the normal operator for smooth velocity models. In the curvelet domain, this smoothness constraint corresponds to limiting the size of the differences amongst neigbor- ing curvelet coefficients. Mathematically, the diagonal estimation can be written in terms of the following min- inmization problem: min λ 1 2 ||Ψr−CTdiag(u)λ ||`2 +κ||Lλ ||`2 , (5) with L = [Dx Dz Dθ ] a so-called sharpening opera- tor penalizing fluctuations amongst neighboring coeffi- cients in λ which is in the curvelet domain. Dx,z contain the first-order differences at each scale in the x,z direc- tions, andDθ the first-order difference in the θ direction. These differences are scale dependent because the spa- tial grid and angles of the curvelet partitioning of phase space both depend on the scale. SPARSITY-CONTINUITY ENHANCING AMPLI- TUDE RECOVERY Normal operator is an important operator since its in- verse can correct for the inaccuracy created by migra- tion process. We introduce a solution for seismic ampli- tude recovery which utilizes the diagonal approximation mentioned in the last section. We can now rewrite the normal equation 2 in the fol- lowing approximate form ' (AATm)+ e (6) = Ax0+ e, with A = CTΓ. Above equation is the direct result of equation 3. Γ = √ Λ is the square-root of the ’eigen- values’ and can be shown to be smooth in the curvelet domain (equation 5). x0 = ATm is called the sparsity vector. We propose an amplitude recovery solution. This solu- tion is formulated in terms of a sparsity-promoting non- linear optimization problem and can be seen as a formal- ization of earlier ideas on stable seismic image recovery. During the optimization, sparsity in the transformed do- main as well as continuity along imaged reflectors, are jointly promoted. Both penalties are part of the follow- ing nonlinear optimization problem (see Herrmann et al. (2006)) P : { x̂=minx J(x) subject to ‖y−Ax‖2 ≤ ε m̂= ( AT )† x̂, (7) where sparsity vector x is optimized with respect to the penalty functional J(x) and the data misfit. We use the term sparsity vector for x to point out that this vector corresponds to the coefficients of a transform that is de- signed to be sparse on the model. The penalty func- tional J(x) is designed to promote sparsity and conti- nuity. The above optimization problem solves for the model by finding a coefficient vector x that minimizes the penalty term subject to fitting the data to within a user-specified tolerance level ε . We reserved the ’tilde’ Robust Seismic Amplitude Recovery Using Curvelets symbol to denote estimates that are found through opti- mization. The recovered modelm is calculated by com- puting the inverse of AT , which represents the diagonally- weighted curvelet synthesis matrix. RECOVERY ALGORITHM Our algorithm for approximately inverting the normal operator involves the following sequence of steps: 1. Form the normal operator using one’s favorite numerical implementation for the migration op- erator and its adjoint, (i.e., Ψ = KTK) . This discrete normal operator needs to be made zero order; 2. Select a reference vector (depth corrected migrated image), r, that is close enough to the unknown image; 3. Estimate the diagonal approximation. This diag- onal approximation defines the synthesis matrix A; 4. Estimate x by solving the nonlinear optimization problem P. This program inverts the synthesis matrix. The discretized model vector m is cal- culated from the recovered coefficient vector x through the pseudo inverse of AT ; EXAMPLES Figures 1(a) and 1(b) show a typical subsalt imaging model. Figure 1(a) shows a velocity model for SEG-AA’ salt model and the reflectivity which can be obtained by differentiating it from its smoothed version used is shown in Figure 1(b). The synthetic data generated us- ing finite difference code (see e.g. Symes, 2006) using velocity model shown in figure 1(a). This synthetic data is muted to avoid direct waves, tapered and bandpass filtered to avoid dipping waves and further migrated as shown in figure 2(a). This migrated image is used in or- der to recover amplitudes according to our method. The depth corrected and bandpass filtered (to remove dip- ping wave) migrated image is used as the reference vec- tor. The reference vector is used for diagonal approxi- mation of normal operator according our method. Sub- sequently, the diagonal approximation is used accord- ing to our nonlinear amplitude recover method to obtain the enhanced image shown in Figure 2(b). For this ex- ample, we used a two-way wave-equation reverse-time migration and modeling (see e.g. Symes, 2006). The dataset consist of 324 shots and 176 receivers for each shot, each trace contains 626 samples with 8 ms sam- pling intervals. CONCLUSION This paper presents a fast and robust approach for ap- proximation of normal operator and utilizes it in ampli- tude recovery in seismic image. We formulated the ap- proximation of the normal operator an eigenvalue-like decomposition with curvelets as eigenvector. Our pro- posed method is faster than other Krylov-based meth- ods since our method evaluates the normal operator only once. The method presented in this paper combines the compression of images by curvelets with the invariance of this transform under the normal operator. This com- bination allows us to formulate a stable recovery method for seismic amplitudes. Compared to other approaches for migration preconditioning, our method (i) brings the amplitude correction problem within the context of sta- ble signal recovery; (ii) provides for a diagonal approx- imation of normal operator. ACKNOWLEDGMENTS The authors would like to thank the authors of Curve- Lab for making their codes available. The authors would also like to thank W. Symes for providing the migration code. This work was in part financially supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant (22R81254) and Collaborative Research and Development Grant DNOISE (334810-05) of Felix J. Herrmann and was carried out as part of the SINBAD project with support, secured through ITF (the Industry Technology Facilitator), from the following or- ganizations: BG Group, BP, Chevron, ExxonMobil and Shell. Robust Seismic Amplitude Recovery Using Curvelets (a) (b) Figure 1: (a) SEG-AA velocity model used to generate synthetic data (b) Reflectivity of the SEG-AA model . Robust Seismic Amplitude Recovery Using Curvelets (a) (b) Figure 2: (a)Migrated synthetic data (b) Enhanced image using our recovery method. Robust Seismic Amplitude Recovery Using Curvelets REFERENCES Candes, E. J. and L. Demanet, 2002, Curvelets and fourier integral operators: Comptes-Rendus de l’Academie des Sciences, I. Chauris, H., 2006, Seismic imaging in the curvelet domain and its implications for the curvelet design: SEG/New Orleans 2006 Annual Meeting. E. Candes, D. D., 2000, Curvelets - a surprisingly effective nonadaptive representation for objects with edges, in curve and surface fitting: Vanderbilt University Press. Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2006, Sparsity- and continuity-promoting seismic image recovery with curvelet frames: in revision. Kuhl, H. and M. Sacchi, 2003, Least-squares wave-equation migration for AVP/AVA inversion: Geophysic, 68, 262– 273. Nemeth., T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysic, 64, 208–221. Smith, H., 1997, A hardy space for fourier integral operaotrs: Geom. Anal., 7, 784–842. ——–, 1998, A parametrix consrtuction for wave equations with c1,1 coefficients: , Ann. Inst. Fourier, Grenoble, 48, 797–835. Symes, W. W., 2006, Reverse time migration with optimal checkpointing: Technical Report TR06-18, Department of Computational and Applied Mathematic, Rice University, Houston, Texas, USA.


Citation Scheme:


Usage Statistics

Country Views Downloads
Japan 2 0
China 1 0
City Views Downloads
Tokyo 2 0
Beijing 1 0

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


Share to:


Related Items