UBC Faculty Research and Publications

Curvelet-domain least-squares migration with sparseness constraints. Herrmann, Felix J.; Moghaddam, Peyman P. 2004

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

Item Metadata


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

Full Text

2881 CURVELET-DOMAIN LEAST-SQUARES MIGRATION WITHSPARSENESS CONSTRAINTSFELIX J. HERRMANN AND PEYMAN MOGHADDAMDepartment of Earth and Ocean Sciences, University of British Columbia, Vancouver, BC, CanadaAbstractA non-linear edge-preserving solution to the least-squares migration problem with sparseness constraints is introduced. The appliedformalism explores Curvelets as basis functions that, by virtue of their sparseness and locality, not only allow for a reduction of thedimensionality of the imaging problem but which also naturally lead to a non-linear solution with significantly improved signal-to-noise ratio. Additional conditions on the image are imposed by solving a constrained optimization problem on the estimatedCurvelet coefficients initialized by thresholding. This optimization is designed to also restore the amplitudes by (approximately)inverting the normal operator, which is like-wise the (de)-migration operators, almost diagonalized by the Curvelet transform.IntroductionLeast-squares migration and migration deconvolution has been a topic that received a recent flare of interest [9, 10]. This interestis for a good reason because inverting for the normal operator (the demigration-migration operator) restores many of the amplitudeartifacts related to acquisition and illumination imprints. However, the downside to this approach is that least-squares tends tosmear the energy leading to a loss of resolution. By imposing certain sparseness constraints on the imaged reflectivity, progress hasbeen made to boost the frequency content of the image [13, 10].This paper comes up with an alternative formulation for the imaging problem designed to (i) deal with substantial amounts ofnoise; (ii) use the optimal (sparse & local) representation properties of Curvelets and their almost diagonalization of the imagingoperators; (iii) use non-linear thresholding techniques, supplemented by constrained optimization on the estimated coefficients,imposing additional sparseness on the model. This paper borrows from ideas by [4] and is an extension of earlier work by [7], inwhich Contourlets [5] were used to denoise and approximately least-square migrate without having access to demigration operators.The paper is organized as follows. First, we briefly discuss the imaging problem and motivate why Curvelets are the appropriatechoice for seismic imaging and processing. We proceed by introducing the non-linear estimation procedure with thresholdingin the image space. We show that Monte-Carlo sampling techniques can be used to compute a correction for the threshold thatincorporates the coloring of the noise due to migration. We conclude by introducing a constrained optimization approach aimedto (approximately, via the diagonal/symbol of the normal operator in the Curvelet domain) invert for the normal operator whileimposing sparseness. The optimization is designed to reduce possible imaging and estimation artifacts, such as side-band effects,erroneous thresholding and bad illumination. The method will be illustrated by a synthetic example using a Kirchoff migrationoperator.The seismic imaging problemIn all generality, the seismic imaging problem can after linearization and high-frequency approximation be cast into the followingform for the forward model describing our datad = Km + n, (1)where K is the (Kirchoff) (de)-migration/scattering operator, m the model with the reflectivity and n white Gaussian noise. Thepertaining inverse problem has the following general form [see e.g. 11, 14]? : minm 12bardbld -Kmbardbl22 + ?J(m), (2)where J(m) is an additional penalty function that contains prior information on the model, such as particular sparseness constraints.The control-parameter ? rules how much emphasis one would like to give to the prior information on the model. How can we findthe appropriate domain to solve this inverse problem?By hitting the data with the migration operator (the adjoint of the scattering operator denoted by asteriskmath), followed by sandwiching thenormal operator between basis-function (de)-compositions, collapses the energy onto a limited number of coefficients. Question ishow to recover these coefficients from the now colored noise,Bu = BFIObracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownrightKasteriskmathKBasteriskmathBm + BPsiDObracehtipdownleftbracehtipuprightbracehtipupleftbracehtipdownrightKasteriskmath n or ? = ? ? + ?. (3)At this point, the cruciality of the right choice for the basis functions becomes apparent. Not only do we wish to represent the modelFigure 1: Properties of the Curvelet Transform under imaging. Curvelet and its demigration (first 2 panels). Clearly, the Curvelet remainsCurvelet-like, which can be explained from its excellent frequency-localization. FK-spectrum is included in the third panel. The rapid decay for the?off-diagonal? coefficients of the Curvelet-images under the (de)-migration and normal operators are shown in the fourth panel.sparsely, but we would also like to correct for the coloring of the noise. Essentially what we want is high non-linear approximationrate for the model and an almost diagonalization of the (de)-migration and normal operators. Mathematically these operatorscorrespond to Fourier Integral (FIO))and Pseudo Differential Operators (PsiDO), respectively. Only under these conditions, we canhope for a significant improvement in the image quality by solving? : minm 12bardbl? - ? ?bardbl22 + ?J(m) or approximately ? : minm 12bardbl? -diag{?} ?bardbl22 + ?J(m). (4)CurveletsCurvelets as proposed by [2, 4], constitute a relatively new family of non-separable wavelet bases that are designed to effectivelyrepresent seismic data with reflectors that generally lie on piece-wise smooth curves (m element C2). Since it has also been shown thatthese basis functions are well behaved under the action of the (de)-migration and normal operators [1] Curvelets - as opposed toWavelets - are the right choice. The well-behavedness corresponds to the property that these basis function remain approximatelyinvariant under the operators as we can see from Fig. 1, where not only Curvelets remain Curvelet-like but where also a rapiddecay for the coefficients of images of Curvelets are observed. Evidently, these properties make Curvelets suitable to be used inimaging since they not only represent the model very well (they obtain near optimal non-linear approximation rates), but they alsoalmost diagonalize the normal operator and hence the Covariance of to the colored noise on the image, Cov?? = E [??asteriskmath]. Withthe appropriate correction for the threshold, we are in the position to relatively easily separate noise from model. How do Curveletsobtain such a high non-linear approximation rate? Without being all inclusive [see for details 2, 4, 5], the answer to this questionlies in the fact that Curvelets are? multi-scale, i.e. they live in different dyadic corona in the 2D Fourier-domain.? multi-directional, i.e. they live on wedges within these corona.? anisotropic, i.e. they obey the following scaling law width proportional length2.? directional selective with #orientations proportional 1radicalscale .? local both in (x,y)) and (kx, ky).? almost orthogonal, they are tight frames.Migration denoising by thresholdingAfter applying migration, the estimate for the u (denoted by ) involves the solution of a denoising problem in the presence ofcolored Gaussian noise. For white Gaussian noise, thresholding on the coefficients solves inverse problems of the type given inEq. 2 by setting K = I. When using Curvelets,? = B-1theta? (Bu) (5)solves the denoising problem for piece-wise smooth reflectors. In this expression, theta is the thresholding operator with threshold? = sigmaradicalbig2loge N, N number of samples and sigma the standard deviation of the noise [see for detail 6, 11, 3]. So how can we correctfor the coloring of the noise knowing that the Covariance is almost diagonal? The answer is simple, simply weight the thresholdwith the square root of the diagonal, i.e.? = etaradicalG with G = diag {Cov??} (6)ans eta an additional control parameter (de)-emphasizing the thresholding, setting the confidence interval (e.g. eta = 3 correspondsto a 95 % confidence interval). Applying this threshold, ? = theta? (?) to the coefficients of the image yields, after inverse-Curvelettransforming, the estimate for the migrated image. Comparison with the original noisy image in Fig. 2 shows a remarkable reductionof the noise. However, the amplitudes have not been recovered accurately because ? still contains the normal operator, ? approxequal ? .Amplitude restoration by constraint optimizationNow that we have successfully removed the noise from the migrated imaged, it is time to restore the amplitudes and impose thesparseness constraint. With the latter, we hope to (i) remove the imaging [by approximately inverting for the normal operator bythe pseudo-inverse of the diagonal G also proposed in 7], estimation and side-band artifacts; (ii) enhance the continuity along theimaged reflectors. To accomplish these goals, we propose the following strategy, which is close to the ideas presented by [4] exceptthat we include imaging operators. Solve the following constrained optimization problem? : minmJ(m) s.t. | ? ? - ?| <=< ? or approximately ? : minmJ(m) s.t. |G? - ?| <=< ?, (7)where, in the approximate case, explicit use has been made of the approximate property that ? approxequal G. This property showsthat G acts as the symbol of the normal operator in the Curvelet domain. We solve this constrained optimization by using anaugmented Lagrangian method involving a Steepest Decent Method [see for details 12]. In the approximate case, there is besidesthe computational benefits, the additional advantage that access to the modelling/demigration operator is no longer necessary. Infact, when provided with G, we are able to approximately turn a vanilla migration into a least-squares migration. Moreover, thefudge factor on estimated coefficients allows for the minimization of the additional penalty function.Remains how to get the G. In case K is a fixed (read non-velocity model dependent) homogeneous operator, such as the Radontransform, the G can be calculated analytically [as shown by 6, 3] and some of these results could potentially carry over to seismicimaging. However, acquisition and illumination imprints remain a problem and we therefore use [see also 7, 8] a Monte-Carlosampling technique to estimate the diagonal of the Covariance operator fromdiag{G} = 1MMsummationdisplayi=1?2i (8)with ?i Curvelet-transformed migrated images of independent realizations of pseudo-random white Gaussian noise. Typically,about 50 realizations are taken.ExampleTo illustrate our proposed method, we included a synthetic example of a common-offset Kirchoff migration for a constant velocitymodel. To show the capabilities of our method, we present our results broadband. We choose 512 sources, time samples, depth andlateral positions. To remove artifacts and improve the continuity along reflectors we used a L1-norm, J(m) = |m|1. As we cansee from the example shown in Fig. 2, thresholding followed by optimization yields a substantially improved image, where most ofthe coherent energy has been removed. Conversely, least-squares migration by Conjugate gradients concentrates the energy of thenoise towards regions that are not well insonified and fails to restore the amplitudes. Our method, on the other hand, removes mostof the artifacts and restores the amplitudes for the larger angles.Conclusions and discussionThe methodology presented is an example of divide and conquer. First, we got a reasonable estimate for the image by thresholdingthen we dealt with amplitude restoration and artifact removal. We build on the premise that one stands a much better chance tosolve an inverse problem by projecting data into the appropriate domain. The combination of migration with the Curvelet transformprovides such a domain, where the energy is collapsed onto a limited number of coefficients that correspond to coherent featuresalong which Curvelets align. Since Curvelets are local in space and spatial frequency, thresholding can be used on the image.By correcting the threshold with the diagonal of the Covariance, the noise is whitened and thresholding successfully removed thebulk of the coherent noise while preserving the reflectors. Imaging amplitudes and continuity were restored by approximatelyinverting the normal operator with its diagonal. Thresholded coefficients were used to start the constrained optimization whichgreatly improved the image quality by imposing the L1-sparseness constraint. Results for the approximate inversion of the normaloperator were slightly better, which can be understood because the diagonal acts as the Curvelet-equivalent for the symbol of thenormal operator. Artifacts generated by the imaging operator can correspond to a non-PsiDO behavior, which may be thought off asoff-diagonal erroneous components. Current research, is focusing on extension to pre-stack imaging; inclusion of different normsand the removal of other coherent noise-sources (multiples etc.) on the image space and the construction of imaging operatorsdirectly in the Curvelet domain.AcknowledgmentsThe authors would like to thank Emmanuel Cand? and David Donoho for making their Digital Curvelet Transforms via Unequally-spaced Fourier Transforms available (presented at the ONR Meeting, University of Minnesota, May, 2003). We also would like tothank Mauricio Sacchi for his migration code. This work was in part financially supported by a NSERC Discovery Grant.References[1] E. J. Cand` and L. Demanet. Curvelets and fourier integral operators. 2002. to appear Comptes-Rendus de l?Academie desSciences, Paris, Serie I.[2] E. J. Cand` and D. Donoho. New tight frames of curvelets and optimal representations of objects with smooth singularities.Technical report, Stanford, 2002. submitted.[3] E. J. Cand` and D. L. Donoho. Recovering Edges in Ill-posed Problems: Optimality of Curvelet Frames. Ann. Statist.,30:784?842, 2000.[4] E. J. Cand` and F. Guo. New multiscale transforms, minimum total variation synthesis: Applications to edge-preservingimage reconstruction. Signal Processing, pages 1519?1543, 2002.Figure 2: Synthetic common-offset Kirchoff-migration example. Top row: Noisy migrated and least-square migratedimage with 10 interations of CG. Notice the coloring of the noise and the failure of CG to correct for the normaloperator. Bottom row: Thresholded image (left) and constrained optimized least-squares image (right). Notice,the significant improvement in the signal-to-noise ratio by the thresholding. The constraint optimization (cf. Eq. 7)removes the artifacts and restores the least-squares amplitudes, leading to a drastically improved image.[5] M. Do and M. Vetterli. Beyond wavelets, chapter Contourlets. Academic Press, 2002.[6] D. Donoho. Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition. App. and Comp. HarmonicAnalysis, 2, 1995.[7] F. J. Herrmann. Multifractional splines: application to seismic imaging. In A. F. L. E. Michael A. Unser, Akram Aldroubi,editor, Proceedings of SPIE Technical Conference on Wavelets: Applications in Signal and Image Processing X, volume 5207,pages 240?258. SPIE, 2003.[8] F. J. Herrmann. Optimal seismic imaging with curvelets. In Expanded Abstracts, Tulsa, 2003. Soc. Expl. Geophys.[9] J. Hu, G. T. Schuster, and P. Valasek. Poststack migration deconvolution. Geophysics, 66(3):939?952, 2001.[10] H. Kuhl and M. D. Sacchi. Least-squares wave-equation migration for avp/ava inversion. Geophysics, 68(1):262?273, 2003.[11] S. G. Mallat. A wavelet tour of signal processing. Academic Press, 1997.[12] J. Nocedal and S. J. Wright. Numerical optimization. Springer, 1999.[13] D. O. Trad. Interpolation and multiple attenuation with migration operators. Geophysics, 68(6):2043?2054, 2003.[14] C. Vogel. Computational Methods for Inverse Problems. SIAM, 2002.


Citation Scheme:


Usage Statistics

Country Views Downloads
China 16 43
United States 8 0
Germany 3 0
France 2 0
Canada 1 0
Japan 1 0
City Views Downloads
Beijing 14 5
Unknown 6 0
Ashburn 4 0
Shenzhen 2 38
Houston 1 0
Ottawa 1 0
Sunnyvale 1 0
University Park 1 0
Tokyo 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