UBC Faculty Research and Publications

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

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

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 Curvelets Peyman P. Moghaddam, Felix J. Herrmann and Chris Stolk  ABSTRACT A non-linear singularity-preserving solution to the least-squares seismic imaging problem with sparseness and continuity constraints is proposed. The applied formalism explores curvelets as a directional frame that, by their sparsity on the image, and their invariance under the imaging operators, allows for a stable recovery of the amplitudes. Our method is based on the estimation of the normal operator in the form of an ’eigenvalue’ decompsoition with curvelets as the ’eigenvectors’. Subsequently, we propose an inversion method that derives from estimation of the normal operator and is formulated as a convex optimization problem. Sparsity in the curvelet domain as well as continuity along the reflectors in the image domain are promoted as part of this optimization. Our method is tested with a reverse-time ’wave-equation’ migration code simulating the acoustic wave equation.  INTRODUCTION Motivated by recent results on stable signal recovery for natural images from incomplete and noisy data (see e.g. Candes et al., 2005), the seismic image recovery problem is formulated as a nonlinear optimization problem. After linearization and by ignoring the source and receiver signatures, the discretized forward model that generates seismic data can be written as d = Km. (1) In this single-scattering expression, m(x) 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 to the data vector (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 normal matrix Ψ = KT K (Nemeth. et al., 1999; Kuhl and Sacchi, 2003) and involves the computation of the Moore-Penrose pseudo inverse (denoted by the symbol † ) of the scattering or demigration matrix K, ΨDO  m = KT K  −1  KT d = K† d.  (3)  FIO  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.  2  et al., 1999). Instead, we argue that a stable amplitude recovery scheme can be obtained when a curvelet decomposition is used that not only sparsely represents the unknown image, but also whose curvelets are invariant under the action of the normal operator. We show that in that case an amplitude weighting (Symes, 2006a) can be performed as part of a nonlinear optimization procedure that exploits these properties and hence adds stability. This stability includes insensitivity to noise and the ability to recover the true amplitudes.  SPARSITY- AND CONTINUITY-PROMOTING IMAGING We address the above issues by exploiting recently developed curvelet frames. These frame expansions compress seismic images and consist of a collection of frame elements ’curvelets’ that are invariant under pseudo-differential operators. These properties allow us to develop an approach similar to the so-called wavelet-vaguelette method (WVD), as proposed by Donoho (1995) and later by Cand`es and Donoho (2000), where scale-invariant homogeneous operators are inverted using the eigenfunction-like behavior of multiscale transforms such as the wavelet and curvelet transform. The 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 e.g. Herrmann et al., 2006) P:  x = minx J(x) subject to † m = AT x,  y − Ax  2  ≤  (4)  where the 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’ symbol (˜) to denote estimates that are found through optimization. The recovered model m is calculated by computing the pseudo inverse (denoted by the symbol †) of AT , which represents the diagonally-weighted curvelet synthesis matrix. This synthesis matrix is designed such that AAT r  Ψr,  (5)  with r an appropriately chosen discrete reference vector and Ψ the discrete normal operator formed by compounding the discrete scattering and its transpose the migration operator. With A = CT Γ and C the curvelet transform, Eq. (5) expresses the normal operator as a form of ’eigenvalue’ decomposition with curvelets as ’eigenvectors’. Γ is the square-root of the ’eigenvalues’ and can be shown to be smooth in the curvelet domain. 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., Ψ = KT K with the symbol. This discrete normal operator needs 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 matrix A;  3  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 AT ;  EXAMPLE Figures 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 shown in Fig.1(b). This reflectivity is de-migrated and migrated (KT K) using the smooth SEG-AA’ background velocity, Fig. 2(a) shows the result of this process. Fig. 2(b) shows the recovered image using our proposed recovery algorithm. Fig.2 (c) shows the trace comparision between the migrated image (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.  CONCLUSION 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. The recovery results show an overall improvement of the image quality. The joined sparsity- and continuity-enhanced image has diminished artifacts, improved resolution and recovered amplitude.  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.  Candes, E., J. Romberg, and T. Tao, 2005, Stable signal recovery from incomplete and inaccurate measurements. to appear in Comm. Pure Appl. Math.  4  (a)  (b)  (c)  Figure 2: Example of seismic amplitude recovery. (a) Normal operator (modeling followed by migration) applied on the SEG-AA’ reflectivity model, (b) recovered image using our proposed method, (c) amplitude comparison for the bottom reflector between original, migrated-demigrated and the recovered image.  5  Cand`es, 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 with curvelet 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 Computational and Applied Mathematics, Rice University, Houston, Texas, USA. ——–, 2006b, Reverse time migration with optimal checkpointing: Technical Report 06-18, Department of Computational 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