UBC Faculty Research and Publications

Curvelet imaging and processing : an overview Herrmann, Felix J. 2004

Item Metadata


52383-Herrmann_F_Curvelet_overview.doc [ 1.27MB ]
JSON: 52383-1.0107384.json
JSON-LD: 52383-1.0107384-ld.json
RDF/XML (Pretty): 52383-1.0107384-rdf.xml
RDF/JSON: 52383-1.0107384-rdf.json
Turtle: 52383-1.0107384-turtle.txt
N-Triples: 52383-1.0107384-rdf-ntriples.txt
Original Record: 52383-1.0107384-source.json
Full Text

Full Text

Curvelet imaging and processing: an overview Felix J. Herrmann Department of Earth and Ocean Sciences, University of British Columbia, Canada  Abstract  In this paper an overview is given on the application of directional basis functions, known under the name Curvelets/Contourlets, to various aspects of seismic processing and imaging. Key conceps in the approach are the use of (i) directional basis functions that localize in both domains (e.g. space and angle); (ii) non-linear estimation, which corresponds to localized muting on the coefficients, possibly supplemented by constrained optimization (iii) invariance of the basis functions under the imaging operators. We will discuss applications that include multiple and ground roll removal; sparseness-constrained least-squares migration and the computation of 4-D difference cubes.  Introduction Everybody would agree that a seismologists dream would be the existence of a basis-function decomposition which localizes in space and angle, is very sparse and well behaved under imaging operators. Well, this dream has come true with the recent introduction of directional wavelets, known under the names Curvelets or Contourlets [5,2,1,4]. As their names suggest, these basis functions are designed to represent images with edges on piece-wise smooth curves very sparsely. In our language, this corresponds to a decomposition that is nearly optimal (read as sparse as you can get) for seismic reflection events as well as for reflectors in the earth. Unfortunately, this property by itself is not enough since we are dealing with imaging operators, which may render the basis-functions less effective. For instance, discrete wavelets transforms have failed to really make a impact in migration because their lack of directional selectivity not only makes them less suitable for migration but they also fail to effectively represent the earth’s reflectors and seismic data. Curvelets (from now on I will use the name Curvelets to refer to both Curvelets and Contourlets), on the other hand, remain more or less invariant under imaging operators while they obtain nearly optimal sparse representations for seismic data as for well images. Our efforts are built on the premise that one is much better suited to solve linear imaging/filtering problems if one is able to very sparsely represent both data and denoised image [7,8]. Therefore, an important part of our research is and will be devoted towards the application of Curvelets to (i) coherent denoising problems such as multiple, ground roll elimination and the robust computation of 4D-difference cubes (refer to other contributions of the author to the Proceedings of this conference) and to (ii) imaging, where the sparseness and invariance of Curvelets are used to improve the noise-robustness and computational efficiency of sparseness-constrained least-squares migration. Our approach differs in several respects from previous attempts to use optimal basis-function decompositions in seismic processing and imaging. First, we do not solely focus on the compression of data or operators. Instead, we concentrate on the design and implementation of basis functions which (i) sparsely represent data and model; and (ii) which are preserved under (de)-migration. This combination of features allows us to reformulate the seismic imaging and inversion problems in terms of non-linear operations on the coefficients of the basis-function decomposition. By replacing linear inversion techniques by non-linear estimation techniques based on thresholding, we limit the smoothing due to regularization and hence we preserve the edges. Therefore, the success of our approach depends primarily on the degree of sparseness we can achieve to represent both data and image. In this abstract, we will present the basic methodology and we mention a number of applications. For more detail, the reader is referred to a series of accompanying abstracts to be presented at this Conference.  Methodology The inverse problem Virtually any linear problems in seismic processing and imaging can be seen as a special case of the following generic problem: how to obtain m form data d in the presence of (coherent) noise n:  d=Km+n,  (1)  in which, for imaging, K is the de-migration operator or the identity for noise-elimination by adaptive subtraction problem. The solution of the inverse problem has the following form [see e.g. [9,11]]  (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 find the appropriate domain to solve this inverse problem? In other words how can exploit the locality of the basis functions in the above setting? It appears from the work of [6], [9] and others that, for a certain class of data and/or images, one can obtain almost optimal denoising results, i.e. near optimal SNR for denoised data, by projecting noisy data onto a basis-function representation that is optimal for that particular class of models. In that case, diagonal decision operator (read simple muting of coefficient) can be constructed that minimizes the energy difference between the estimate and the true model, i.e. Grea t Exp lo r a t i o n s  – Canada and Beyond  1  findD: =D(d) such that minD|-m|22.  (3)  For basis functions that are local, one can show that simple thresholding on the coefficients [9,6] suffices to solve the above denoising problem, i.e.  (4) In this expression, B, B-1 refer to the basis-function expansion and its (pseudo)-inverse. θ is a soft/hard thresholding operator with a ν threshold that for orthonormal basis functions equals [9,6]  (5) with σ the standard deviation of the noise and N the number of data samples. Comparison of equations (2) and (4) establishes the connection between the threshold and the control parameter ν. Question of course now is: can we extend these results, proven to be valid for K=I and n white Gaussian noise, to include the effects of (i) an operator (e.g. K a de-migration operator) and/or (ii) of colored additive noise (e.g. coherent ground roll and multiples)? In other words, can we decolor predicted coherent noise such that the above thresholding procedure becomes effective again? Before answering this question let us first be more specific with respect to the choice of the appropriate basis functions for seismic data, model and coherent noise. To wet our appetite, I included in Fig. 3 an example of thresholding, which clearly demonstrates the potential of the approach applied to the coefficients of Curvelet transforms performed on constant-angle-sorted Common Angle Gathers in combination with a wavelet transform in the angle direction. Even though, we ignored in this example the influence of coloring of noise by migration we greatly improved the stack. The basis functions Curvelets as proposed by [5] and [2,1,4] respectively, constitute a relatively new family of non-separable wavelet bases that are designed to effectively represent functions with singularities (read wavefronts) that lie on piece-wise smooth curves. For these type of functions, nonseparable basis functions obtain nearly optimal non-linear approximation rates (NAR). For comparison, the m-term expansion of the sorted 2D Fourier transform obtains for the object x [3,4,5] the following NAR (measured in the L2-norm),  |x-m|22∝m-1/2  while separable wavelets (ordinary 2-D discrete wavelet transforms) obtain  |x-m|22∝m-1.  The improvement for wavelets has been one of the main reasons of their success in image processing. Clearly, the further improvement in NRA by Curvelets, which obtain  |x-m|22∝m-2(log m)3,  will likely open up a whole new series of successful applications. Besides the logarithmic term, the above NRA is equivalent to the best possible approximation rate (e.g. by a triangular meshing given the location of the edges/reflection events) attainable for this type functions. In Donoho’s words the true miracle of wavelets is that they find the singularities (read wavefronts/reflection events) at virtually no additional cost. While this statement holds for piece-wise smooth 1-D functions, it does not hold for curved events. Curvelets, on the other hand, remedy the wavelets as well the Fourier shortcomings. For example, a simple spike gives rise to a full Fourier spectrum. So how do Curvelets obtain such a high non-linear approximation rate? without being all inclusive [see for details [2,1,4,5]], the answer to this question lies in the fact that Curvelets are  • • • • • •  multi-scale, i.e. they live in different dyadic corona (see Figure 1) in the FK-domain. multi-directional, i.e. they live on wedges within these corona (see Figure 1). anisotropic, i.e. they obey the following scaling law width∝length2. directional selective with . local both in (x,t) and (k,f). almost orthogonal, they are tight frames with a moderate redundancy. Contourlets implement the pseudo-inverse in closed-form while Curvelets provide the transform and its adjoint, yielding a pseudo-inverse computed by iterative Conjugate Gradients.  In Figure 1, a number of Curvelet properties are detailed [adapted from [3,4,5]]. Figure 1 describes the Curvelet partitioning of the frequencywavenumber plane. One Curvelet lives in a wedge and becomes more directional selective and anisotropic for the higher frequencies (while Grea t Exp lo r a t i o n s  – Canada and Beyond  2  moving away from the origin). Curvelets are localized in both the space (or (x,t)) and spatial (KF)-domains and have, as consequence of their partitioning, the tendency to align themselves with curves/wavefronts (see Figure 1). As such they are more flexible then a representation yielded by e.g. high-resolution Radon [as described by e.g. [10]], because they are local and able to follow any piece-wise smooth curve. Only Curvelets that align with reflectors yield large inner products and this explains their success in solving the denoising problem for models that contain events on curves. The noise is canceled because the basis functions adapt themselves locally to the dip and hence optimizing their overlap with the data integrating out the incoherent noise. Now what about coherent noise such as predicted multiples?  Figure 1: Left: Curvelet partitioning of the frequency plane [modified from [3]]. Right: Comparison of non-linear approximation rates Curvelets and Wavelets [modified from [5]].  Figure 2: Properties of the Curvelet Transform under imaging. Curvelet and its demigration. Clearly, the Curvelet remains Curvelet-like, which can be explained from its excellent frequency-localization. Far right panel illustrates why Curvelet coefficients are large when they align and small when they are not.  Thresholding, constrained optimization and sparseness constraints Given the Curvelets and the proposed thresholding procedure how can we operate in an environment where there is coherent additive noise? Well in cases where we have reasonable accurate predictions for the noise, we can simply weight the above threshold (cf. Eq. 5) with the Curvelet coefficients of the predicted noise, i.e.  µ∼η |Bn|.  (6)  Given the predicted noise, we are now able to adaptively decide whether a certain event belongs to signal or to the noise. The n contains the predicted noise and η represents an additional control parameter which sets the confidence interval (e.g. 95 % for η=3) (de)-emphasizing the thresholding. Even though, thresholding takes care of the bulk of the separation of noise and signal, there remains the desire to (i) correct for an operator, if present, to restore the amplitudes; (ii) impose additional prior information on the model, e.g. certain sparseness constraints. With the latter, we hope to remove possible estimation and side-band artifacts, enhancing the continuity along the imaged reflectors. To accomplish these goals, we formulate the following constrained optimization problem [4]  : minmJ(m) s.t. |Ã-|≤µ ,  (7)  where our initial estimate by thresholding, , is updated, subject to constraints, as to minimize the penalty function. The symbol is used to refer to the coefficients and operators in the Curvelet domain.  Applications The methodology presented in this paper has been applied to several areas in seismic processing and imaging. The applications can roughly be divided into adaptive subtraction (K=I hence Ã=I) for Grea t Exp lo r a t i o n s  – Canada and Beyond  3  • multiple elimination, where predicted multiples and multiple-free data are represented by the coherent noise n and model m, respectively. • ground roll elimination, where ground roll-free data is the model m and ground roll the noise n. • computation 4-D difference cubes, where one vintage dataset is the noise of the other and vice versa. ∗  and sparseness-constrained least-squares migration where d is migrated data, A=K K, the normal operator and n colored noise (by migration). See for details, other contributions by the authors to the Proceedings of this Conference.  Conclusions Why do Curvelets seem to be the appropriate choice to formulate problems in seismic processing and imaging.The answer to this question is relatively simple. Curvelets provide an almost orthogonal decomposition w.r.t. to basis functions that have far superior (compared to e.g. Fourier, Radon, Wavelet) localization and sparseness properties then other transforms. In effect, they can for the high frequencies be seen as local Radon transforms with perfect reconstruction. Their locality and directional selectivity renders local thresholding/muting very effective, while their almost diagonalization of imaging operators potentially makes these basis functions a good choice to solve imaging problems, including the construction of fast migration operators.  Acknowledgments The authors would like to thank Emmanuel Candés 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 to thank Mauricio Sacchi for his migration code. This work was in part financially supported by a NSERC Discovery Grant.  Figure 3: Real data example denoised by a simple thresholding on the coefficients after Contourlet and then Wavelet transforming the image for multiple angles.Top left: original stack. Top right: denoised stack (within the window). Bottom middle: predicted noise. Notice the improved SNR and removal of migration artifacts. References  [1]  E. J. Candès and D. Donoho. New tight frames of curvelets and optimal representations of objects with smooth singularities. Technical report, Stanford, 2002. submitted. [2] E. J. Candès and D. L. Donoho. Curvelets – a surprisingly effective nonadaptive representation for objects with edges. Curves and Surfaces. Vanderbilt University Press, 2000. [3] E. J. Candès and D. L. Donoho. Recovering Edges in Ill-posed Problems: Optimality of Curvelet Frames. Ann. Statist., 30: 784–842, 2000. [4] E. J. Candès and F. Guo. New multiscale transforms, minimum total variation synthesis: Applications to edge-preserving image reconstruction. Signal Processing, pages 1519–1543, 2002. [5] M. Do and M. Vetterli. Beyond wavelets, chapter Contourlets. Academic Press, 2002. [6] D. L. Donoho and I. M. Johnstone. Minimax estimation via wavelet shrinkage. Annals of Statistics, 26(3): 879–921, 1998. [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] S. G. Mallat. A wavelet tour of signal processing. Academic Press, 1997. [10] D. Trad, T. Ulrych, and M. Sacchi. Latest views of the sparse radon transform. Geophysics, 68(1): 386–399, 2003. [11] C. Vogel. Computational Methods for Inverse Problems. SIAM, 2002.  Grea t Exp lo r a t i o n s  – Canada and Beyond  4  


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