UBC Faculty Research and Publications

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

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

Item Metadata


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

Full Text

Robust seismic amplitude recovery using curveletsPeyman P. Moghaddamasteriskmath, Felix Herrmann, University of British Columbia, Canada and Chris Stolk, Uni-versity of Twente, NetherlandsSUMMARYIn this paper, we recover the amplitude of a seis-mic image by approximating the normal (demigration-migration) operator. In this approximation, we makeuse of the property that curvelets remain invariant underthe action of the normal operator. We propose a seis-mic amplitude recovery method that employs an eigen-value like decomposition for the normal operator usingcurvelets as eigen-vectors. Subsequently, we proposean approximate non-linear singularity-preserving solu-tion to the least-squares seismic imaging problem withsparseness in the curvelet domain and spatial continu-ity constraints. Our method is tested with a reverse-time?wave-equation? migration code simulating the acousticwave equation on the SEG-AA salt model.INTRODUCTIONIn the mid-90s, Hart Smith (Smith, 1997, 1998) intro-duced waveforms that are invariant under the action ofcertain Fourier and pseudo-differential operators. Basedon 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 andimage processing. Candes and Demanet (2002) showedhow curvelets behave under the action of Fourier inte-gral operators. They stated that under certain conditionsa curvelet maps to another curvelet under the action ofFourier 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 curveletsare discrete and hence move around on grids, this makesit a challenge to define fast migration operators. Curvelets,however, prove to be very useful for solving the seismicamplitude recovery problem, during which curvelets arebeing imaged. This seismic amplitude recovery probleminvolves 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 approximateand subsequently approximately invert the normal oper-ator, an approach followed in this abstract. During theinversion curvelet sparsity and continuity along imagedreflectors are jointly promoted leading to a stable recov-ery of the seismic amplitudes in particular along steeplydipping events and under the salt.THE PROBLEM FORMULATIONAfter linearization and by ignoring the source and re-ceiver signatures, the discretized forward model that gen-erates seismic data can be written asd =Km+n, (1)with d the first-order Born approximation data, K thescattering operator (or commonly called de-migration ormodeling operator), n is the noise and m represents the(singular) fluctuations in the earth?s acoustic propertieswith respect to an appropriately chosen smoothly vary-ing background velocity model (the density of mass isassumed constant). These fluctuations are referred to asthe model and seismic imaging aims to recover both thelocations and the relative amplitudes of the velocity fluc-tuations from seismic data. Applying the adjoint of thelinearized scattering operator K to the data vector ( d inequation 1) leads to the migrated image.KT d = KT Km+KT n (2)y = parenleftbigPsimparenrightbig +e,with Psi =KT K defined as normal operator and e =KT nis the colored noise.An extensive literature has emerged on restoring the mi-gration amplitude by inverting the normal matrix (Psi =KT K) (Nemeth. et al., 1999; Kuhl and Sacchi, 2003;Herrmann et al., 2006), and involves the computation ofthe pseudo inverse of the scattering or demigration ma-trix. Unfortunately, the normal operator it too big to beconstructed explicitly and is too expensive to be eval-uated as part of an iterative Krylov-subspace solver toinvert the Hessian (see e.g. Nemeth. et al., 1999).THEORY: DIAGONAL APPROXIMATION OF NOR-MAL OPERATORA curvelet phi? is defined by its index ? which is a triple( j, k, l) with j =0, 1, 2, ... a scale index; l =0, 1, ..., 2floorleft j/2floorright -1 an orientation index (floorleftxfloorright is the integer part of x ); andRobust Seismic Amplitude Recovery Using Curveletsk = (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., Psi ) on a curvelet phi? is ap-proximately mapped into the same curvelet . To be morespecific, a pseudo-differential operator induces a map-ping ? mapstoarrowright ?prime with property that significant curvelet co-efficients of Psi(phi?) are located very close to ? , itself.We propose following decomposition of the normal op-eratorPsi approxequalCT ?C, (3)where ? is a diagonal matrix, C and CT are the curveletand the transpose of curvelet transform. This decom-position states that we can decompose the normal oper-ator in a form of eigenvalue-like decomposition wherecurvelets play the role of eigenfunctions and the diago-nal ? the role of the eigenvalues.Since the eigenvalues are not known, we propose a methodthat estimates this diagonal from the action of the normaloperator on a reference vector (r) that is close enough tothe original model.Psir =CT ?Cr mapstoarrowright CT diag(u)lambda =Psir, (4)with u =Cr and lambda is the diagonal elements of ?.During this diagonal estimation method, we aim to finda curvelet-domain diagonal weighting vector (i.e., lambda)that approximates the action of the normal operator onthe reference vector. We choose the spherical-spreadingcorrected 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 imposeadditional smoothing constrains on the to-be-estimateddiagonal (lambda). These constraints reflect the behavior ofthe normal operator for smooth velocity models. In thecurvelet domain, this smoothness constraint correspondsto limiting the size of the differences amongst neigbor-ing curvelet coefficients. Mathematically, the diagonalestimation can be written in terms of the following min-inmization problem:minlambda12||PsiPsiPsir-CCCT diag(u)lambda||lscript2 +kappa||LLLlambda||lscript2 , (5)with L = [Dx Dz Dtheta ] a so-called sharpening opera-tor penalizing fluctuations amongst neighboring coeffi-cients in lambda which is in the curvelet domain. Dx,z containthe first-order differences at each scale in the x, z direc-tions, and Dtheta the first-order difference in the theta direction.These differences are scale dependent because the spa-tial grid and angles of the curvelet partitioning of phasespace both depend on the scale.SPARSITY-CONTINUITY ENHANCING AMPLI-TUDE RECOVERYNormal 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 approximationmentioned in the last section.We can now rewrite the normal equation 2 in the fol-lowing approximate formsimilarequal parenleftbigAAT mparenrightbig +e (6)= Ax0 +e,with A = CT G. Above equation is the direct result ofequation 3. G = radical? is the square-root of the ?eigen-values? and can be shown to be smooth in the curveletdomain (equation 5). x0 = AT m is called the sparsityvector.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, arejointly promoted. Both penalties are part of the follow-ing nonlinear optimization problem (see Herrmann et al.(2006))P:braceleftBigg? =minx J(x) subject to bardbly-Axbardbl2 <=<epsilon? =parenleftbigAT parenrightbig? ?,(7)where sparsity vector x is optimized with respect to thepenalty functional J(x) and the data misfit. We use theterm sparsity vector for x to point out that this vectorcorresponds 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 themodel by finding a coefficient vector x that minimizesthe penalty term subject to fitting the data to within auser-specified tolerance level epsilon . We reserved the ?tilde?Robust Seismic Amplitude Recovery Using Curveletssymbol to denote estimates that are found through opti-mization. The recovered model m is calculated by com-puting the inverse of ATTT , which represents the diagonally-weighted curvelet synthesis matrix.RECOVERY ALGORITHMOur algorithm for approximately inverting the normaloperator involves the following sequence of steps:1. Form the normal operator using one?s favoritenumerical implementation for the migration op-erator and its adjoint, (i.e., Psi = KTTT K) . Thisdiscrete normal operator needs to be made zeroorder;2. Select a reference vector (depth corrected migratedimage), r, that is close enough to the unknownimage;3. Estimate the diagonal approximation. This diag-onal approximation defines the synthesis matrixA;4. Estimate x by solving the nonlinear optimizationproblem P. This program inverts the synthesismatrix. The discretized model vector m is cal-culated from the recovered coefficient vector xthrough the pseudo inverse of AT ;EXAMPLESFigures 1(a) and 1(b) show a typical subsalt imagingmodel. Figure 1(a) shows a velocity model for SEG-AA?salt model and the reflectivity which can be obtainedby differentiating it from its smoothed version used isshown in Figure 1(b). The synthetic data generated us-ing finite difference code (see e.g. Symes, 2006) usingvelocity model shown in figure 1(a). This synthetic datais muted to avoid direct waves, tapered and bandpassfiltered to avoid dipping waves and further migrated asshown in figure 2(a). This migrated image is used in or-der to recover amplitudes according to our method. Thedepth 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 obtainthe enhanced image shown in Figure 2(b). For this ex-ample, we used a two-way wave-equation reverse-timemigration and modeling (see e.g. Symes, 2006). Thedataset consist of 324 shots and 176 receivers for eachshot, each trace contains 626 samples with 8 ms sam-pling intervals.CONCLUSIONThis 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-likedecomposition with curvelets as eigenvector. Our pro-posed method is faster than other Krylov-based meth-ods since our method evaluates the normal operator onlyonce. The method presented in this paper combines thecompression of images by curvelets with the invarianceof this transform under the normal operator. This com-bination allows us to formulate a stable recovery methodfor seismic amplitudes. Compared to other approachesfor migration preconditioning, our method (i) brings theamplitude correction problem within the context of sta-ble signal recovery; (ii) provides for a diagonal approx-imation of normal operator.ACKNOWLEDGMENTSThe authors would like to thank the authors of Curve-Lab for making their codes available. The authors wouldalso like to thank W. Symes for providing the migrationcode. This work was in part financially supported by theNatural Sciences and Engineering Research Council ofCanada Discovery Grant (22R81254) and CollaborativeResearch and Development Grant DNOISE (334810-05)of Felix J. Herrmann and was carried out as part of theSINBAD project with support, secured through ITF (theIndustry Technology Facilitator), from the following or-ganizations: BG Group, BP, Chevron, ExxonMobil andShell.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 CurveletsREFERENCESCandes, E. J. and L. Demanet, 2002, Curvelets and fourier integral operators: Comptes-Rendus de l?Academie desSciences, I.Chauris, H., 2006, Seismic imaging in the curvelet domain and its implications for the curvelet design: SEG/NewOrleans 2006 Annual Meeting.E. Candes, D. D., 2000, Curvelets - a surprisingly effective nonadaptive representation for objects with edges, in curveand surface fitting: Vanderbilt University Press.Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2006, Sparsity- and continuity-promoting seismic image recoverywith 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 ofComputational and Applied Mathematic, Rice University, Houston, Texas, USA.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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