UBC Faculty Research and Publications

Seismic Amplitude Recovery with Curvelets Moghaddam, Peyman P.; Herrmann, Felix J.; Stolk, Christiaan C. 2007-03-07

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

Item Metadata


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

Full Text

Seismic Amplitude Recovery with CurveletsPeyman P. Moghaddam, Felix J. Herrmann and Chris StolkABSTRACTA non-linear singularity-preserving solution to the least-squares seismic imaging problem withsparseness and continuity constraints is proposed. The applied formalism explores curvelets asa directional frame that, by their sparsity on the image, and their invariance under the imagingoperators, allows for a stable recovery of the amplitudes. Our method is based on the estimationof the normal operator in the form of an ?eigenvalue? decompsoition with curvelets as the?eigenvectors?. Subsequently, we propose an inversion method that derives from estimationof the normal operator and is formulated as a convex optimization problem. Sparsity in thecurvelet domain as well as continuity along the reflectors in the image domain are promoted aspart of this optimization. Our method is tested with a reverse-time ?wave-equation? migrationcode simulating the acoustic wave equation.INTRODUCTIONMotivated by recent results on stable signal recovery for natural images from incomplete and noisydata (see e.g. Candes et al., 2005), the seismic image recovery problem is formulated as a nonlinearoptimization problem. After linearization and by ignoring the source and receiver signatures, thediscretized forward model that generates seismic data can be written asd = Km. (1)In this single-scattering expression, m(x) represents the (singular) fluctuations in the earth?s acous-tic 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 andseismic 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 to the datavector (d in Eq. (1)) leads to the migrated image,y = KT d. (2)An extensive literature has emerged on restoring the migration amplitude by inverting the nor-mal matrix Psi = KT K(Nemeth. et al., 1999; Kuhl and Sacchi, 2003) and involves the computationof the Moore-Penrose pseudo inverse (denoted by the symbol ?) of the scattering or demigrationmatrix K,m=PsiDObracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownrightparenleftbigKT Kparenrightbig-1 KTbracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipuprightFIOd = K?d. (3)Unfortunately, the normal operator it too big to be constructed explicitly and is too expensive tobe evaluated as part of an iterative Krylov-subspace solver to invert the Hessian (see e.g. Nemeth.2et al., 1999). Instead, we argue that a stable amplitude recovery scheme can be obtained whena curvelet decomposition is used that not only sparsely represents the unknown image, but alsowhose curvelets are invariant under the action of the normal operator. We show that in that case anamplitude weighting (Symes, 2006a) can be performed as part of a nonlinear optimization procedurethat exploits these properties and hence adds stability. This stability includes insensitivity to noiseand the ability to recover the true amplitudes.SPARSITY- AND CONTINUITY-PROMOTING IMAGINGWe address the above issues by exploiting recently developed curvelet frames. These frame ex-pansions compress seismic images and consist of a collection of frame elements ?curvelets? thatare invariant under pseudo-differential operators. These properties allow us to develop an approachsimilar to the so-called wavelet-vaguelette method (WVD), as proposed by Donoho (1995) and laterby Cand` and Donoho (2000), where scale-invariant homogeneous operators are inverted using theeigenfunction-like behavior of multiscale transforms such as the wavelet and curvelet transform.The solution is formulated in terms of a sparsity-promoting nonlinear optimization problemand can be seen as a formalization of earlier ideas on stable seismic image recovery. During theoptimization, sparsity in the transformed domain as well as continuity along imaged reflectors, arejointly promoted. Both penalties are part of the following nonlinear optimization problem (see e.g.Herrmann et al., 2006)P :braceleftBiggtildewide = minx J(x) subject to bardbly -Axbardbl2 <=< epsilon1tildewide = parenleftbigAT parenrightbig? tildewide, (4)where the sparsity vector x is optimized with respect to the penalty functional J(x) and the datamisfit. We use the term sparsity vector for x to point out that this vector corresponds to the coeffi-cients 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 opti-mization problem solves for the model by finding a coefficient vector x that minimizes the penaltyterm subject to fitting the data to within a user-specified tolerance level epsilon1. We reserved the ?tilde?symbol (?) to denote estimates that are found through optimization. The recovered model m iscalculated by computing the pseudo inverse (denoted by the symbol ?) of AT , which represents thediagonally-weighted curvelet synthesis matrix. This synthesis matrix is designed such thatAAT r similarequal Psir, (5)with r an appropriately chosen discrete reference vector and Psithe discrete normal operator formedby compounding the discrete scattering and its transpose the migration operator. With A = CT Gand C the curvelet transform, Eq. (5) expresses the normal operator as a form of ?eigenvalue?decomposition with curvelets as ?eigenvectors?. G is the square-root of the ?eigenvalues? and can beshown to be smooth in the curvelet domain. Our algorithm for approximately inverting the normaloperator involves the following sequence of steps:1. Form the normal operator using one?s favorite numerical implementation for the migrationoperator and its adjoint, i.e., Psi = KTK with the symbol. This discrete normal operatorneeds to be made zero order;2. Select a relevant reference vector, r, that is close enough to the unknown image;3. Estimate the diagonal approximation. This diagonal approximation defines the synthesis ma-trix A;34. Estimate x by solving the nonlinear optimization problem P. This program inverts the syn-thesis matrix. The discretized model vector m is calculated from the recovered coefficientvector x through the pseudo inverse of AT ;EXAMPLEFigures 1-2 show the evaluation of our method on a typical subsalt imaging problem. Fig. 1(a)shows a sufficiently smooth background velocity for SEG-AA? salt model, with reflectivity shownin Fig.1(b). This reflectivity is de-migrated and migrated (KT K) using the smooth SEG-AA? back-ground velocity, Fig. 2(a) shows the result of this process. Fig. 2(b) shows the recovered imageusing our proposed recovery algorithm. Fig.2 (c) shows the trace comparision between the migratedimage (properly scaled) (Fig.2(a)), original reflectivity (Fig.1(b)) and recovered image (Fig.2(b))along the horizontal line in the bottom of image (3500 m depth) for the offsets from 4.3 to 7.2 km.For this example, we used a two-way wave-equation reverse-time migration and modeling (Symes(2006b)). The dataset consist of 324 shots and 176 receivers for each shot.(a) (b)Figure 1: (a) SEG-AA Sufficiently smooth background model, (b) Reflectivity model.CONCLUSIONThe method presented in this paper combines the compression of images by curvelets with theinvariance of this transform under the normal operator. This combination allows us to formulatea stable recovery method for seismic amplitudes. Compared to other approaches for migrationpreconditioning, our method (i) brings the amplitude correction problem within the context of stablesignal recovery; (ii) provides for a diagonal approximation of normal operator. The recovery resultsshow an overall improvement of the image quality. The joined sparsity- and continuity-enhancedimage has diminished artifacts, improved resolution and recovered amplitude.ACKNOWLEDGMENTSThe authors would like to thank the authors of CurveLab for making their codes available. The authors would alsolike to thank W. Symes for providing the migration code. This work was in part financially supported by the NaturalSciences and Engineering Research Council of Canada Discovery Grant (22R81254) and Collaborative Research andDevelopment Grant DNOISE (334810-05) of Felix J. Herrmann and was carried out as part of the SINBAD project withsupport, secured through ITF (the Industry Technology Facilitator), from the following organizations: BG Group, BP,Chevron, ExxonMobil and Shell.Candes, E., J. Romberg, and T. Tao, 2005, Stable signal recovery from incomplete and inaccurate measurements. toappear in Comm. Pure Appl. Math.4(a)(b)(c)Figure 2: Example of seismic amplitude recovery. (a) Normal operator (modeling followed by migration) appliedon the SEG-AA? reflectivity model, (b) recovered image using our proposed method, (c) amplitude comparison for thebottom reflector between original, migrated-demigrated and the recovered image.5Cand` E. J. and D. L. Donoho, 2000, Recovering Edges in Ill-posed Problems: Optimality of Curvelet Frames: Ann.Statist., 30, 784?842.Donoho, D. L., 1995, De-noising by soft thresholding: 41, 613?627.Herrmann, F. J., P. P. Moghaddam, and C. Stolk, 2006, Sparsity- and continuity-promoting seismic image recovery withcurvelet frames: under revision.Kuhl, H. and M. Sacchi, 2003, Least-squares wave-equation migration for AVP/AVA inversion: Geophysics, 68, 262?273.Nemeth., T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysics, 64,208?221.Symes, W. W., 2006a, Optimal scaling for reverse time migration: Technical Report TR 06-19, Department of Computa-tional and Applied Mathematics, Rice University, Houston, Texas, USA.???, 2006b, Reverse time migration with optimal checkpointing: Technical Report 06-18, Department of Computa-tional and Applied Mathematics, 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