UBC Faculty Research and Publications

Robust seismic amplitude recovery using curvelets Moghaddam, Peyman P.; Herrmann, Felix J.; Stolk, Christiaan C. 2007-12-31

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

Item Metadata


moghaddam07seg.pdf [ 166.73kB ]
JSON: 1.0107414.json
JSON-LD: 1.0107414+ld.json
RDF/XML (Pretty): 1.0107414.xml
RDF/JSON: 1.0107414+rdf.json
Turtle: 1.0107414+rdf-turtle.txt
N-Triples: 1.0107414+rdf-ntriples.txt
Original Record: 1.0107414 +original-record.json
Full Text

Full Text

Robust seismic amplitude recovery using curvelets Peyman P. Moghaddam∗ , Felix Herrmann, University of British Columbia, Canada and Chris Stolk, University of Twente, Netherlands SUMMARY In this paper, we recover the amplitude of a seismic image by approximating the normal (demigrationmigration) operator. In this approximation, we make use of the property that curvelets remain invariant under the action of the normal operator. We propose a seismic amplitude recovery method that employs an eigenvalue like decomposition for the normal operator using curvelets as eigen-vectors. Subsequently, we propose an approximate non-linear singularity-preserving solution to the least-squares seismic imaging problem with sparseness in the curvelet domain and spatial continuity 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) introduced waveforms that are invariant under the action of certain Fourier and pseudo-differential operators. Based on his theory, E. Candes (2000) designed and implemented 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 integral 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 theoretical invariance property towards a diagonalization of migration 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 nonstationary convolution (i.e., pseudodifferential operator). Since curvelets in this case no longer move, their relative invariance can be used to diagonally approximate and subsequently approximately invert the normal operator, 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 receiver signatures, the discretized forward model that generates seismic data can be written as d = K m + 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 varying 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 fluctuations 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. KT d = KT Km + KT n y = Ψ m + e,  (2)  with Ψ = K T K defined as normal operator and e = K T n is the colored noise. An extensive literature has emerged on restoring the miΨ= gration amplitude by inverting the normal matrix (Ψ K T K ) (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 matrix. Unfortunately, the normal operator it too big to be constructed explicitly and is too expensive to be evaluated as part of an iterative Krylov-subspace solver to invert the Hessian (see e.g. Nemeth. et al., 1999). THEORY: DIAGONAL APPROXIMATION OF NORMAL 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, ..., 2 j/2 − 1 an orientation index ( x is the integer part of x ); and  Robust Seismic Amplitude Recovery Using Curvelets k = (k1 , k2 ) is the location index. One of the important features of curvelet is that the action of a pseudodifferential operator (i.e., Ψ ) on a curvelet φµ is approximately mapped into the same curvelet . To be more specific, a pseudo-differential operator induces a mapping µ → µ with property that significant curvelet coefficients of Ψ (φµ ) are located very close to µ , itself. We propose following decomposition of the normal operator T  Ψ ≈ C ΛC ,  (3)  where Λ is a diagonal matrix, C and C T are the curvelet and the transpose of curvelet transform. This decomposition states that we can decompose the normal operator in a form of eigenvalue-like decomposition where curvelets play the role of eigenfunctions and the diagonal Λ the role of the eigenvalues.  the first-order differences at each scale in the x, z directions, and D θ the first-order difference in the θ direction. These differences are scale dependent because the spatial grid and angles of the curvelet partitioning of phase space both depend on the scale.  SPARSITY-CONTINUITY ENHANCING AMPLITUDE RECOVERY Normal operator is an important operator since its inverse can correct for the inaccuracy created by migration process. We introduce a solution for seismic amplitude recovery which utilizes the diagonal approximation mentioned in the last section. We can now rewrite the normal equation 2 in the following approximate form AAT m + e = A x0 + e,  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 = C T ΛC r → C T diag(u)λ = Ψ r,  (4)  with u = C r 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. Because the curvelet transform is redundant, the estimation of this diagonal in equation refeqn:refvec is an underdetermined inverse problem. To find an unique solution 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 neigboring curvelet coefficients. Mathematically, the diagonal estimation can be written in terms of the following mininmization problem: 1 Ψr −C C T diag(u)λ || 2 + κ||L Lλ || 2 , min ||Ψ λ 2  (5)  Dx D z D θ ] a so-called sharpening operawith L = [D tor penalizing fluctuations amongst neighboring coefficients in λ which is in the curvelet domain. D x,z contain  (6)  with A = C T Γ. Above equation is the direct result of √ equation 3. Γ = Λ is the square-root of the ’eigenvalues’ and can be shown to be smooth in the curvelet domain (equation 5). x0 = A T m is called the sparsity vector. We propose an amplitude recovery solution. This solution is formulated in terms of a sparsity-promoting nonlinear optimization problem and can be seen as a formalization of earlier ideas on stable seismic image recovery. During the optimization, sparsity in the transformed domain as well as continuity along imaged reflectors, are jointly promoted. Both penalties are part of the following nonlinear optimization problem (see Herrmann et al. (2006))  P:  xˆ = minx J(x) subject to † ˆ = AT xˆ , m  y − Ax  2  ≤ε  (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 designed to be sparse on the model. The penalty functional J(x) is designed to promote sparsity and continuity. 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 optimization. The recovered model m is calculated by computing the inverse of A T , which represents the diagonallyweighted 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 operator and its adjoint, (i.e., Ψ = K T K ) . 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 diagonal 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 calculated from the recovered coefficient vector x through the pseudo inverse of A T ;  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 using 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 order to recover amplitudes according to our method. The depth corrected and bandpass filtered (to remove dipping wave) migrated image is used as the reference vector. The reference vector is used for diagonal approximation of normal operator according our method. Subsequently, the diagonal approximation is used according to our nonlinear amplitude recover method to obtain the enhanced image shown in Figure 2(b). For this example, 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 sampling intervals.  CONCLUSION This paper presents a fast and robust approach for approximation of normal operator and utilizes it in amplitude recovery in seismic image. We formulated the approximation of the normal operator an eigenvalue-like decomposition with curvelets as eigenvector. Our proposed method is faster than other Krylov-based methods 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 combination 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 stable signal recovery; (ii) provides for a diagonal approximation of normal operator. ACKNOWLEDGMENTS The authors would like to thank the authors of CurveLab 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 organizations: 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:


Citations by CSL (citeproc-js)

Usage Statistics

Country Views Downloads
China 24 0
Brazil 13 3
Germany 12 0
United States 10 0
Japan 2 0
Indonesia 2 0
Canada 1 1
Russia 1 0
City Views Downloads
Unknown 27 3
Shenzhen 23 0
Ashburn 9 0
Tokyo 2 0
University Park 1 0
Ottawa 1 1
Beijing 1 0
Saint Petersburg 1 0

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



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