UBC Faculty Research and Publications

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

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

Item Metadata


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

Full Text

2881  CURVELET-DOMAIN LEAST-SQUARES MIGRATION WITH SPARSENESS CONSTRAINTS FELIX J. HERRMANN AND PEYMAN MOGHADDAM Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, BC, Canada  Abstract A non-linear edge-preserving solution to the least-squares migration problem with sparseness constraints is introduced. The applied formalism explores Curvelets as basis functions that, by virtue of their sparseness and locality, not only allow for a reduction of the dimensionality of the imaging problem but which also naturally lead to a non-linear solution with significantly improved signalto-noise ratio. Additional conditions on the image are imposed by solving a constrained optimization problem on the estimated Curvelet 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.  Introduction Least-squares migration and migration deconvolution has been a topic that received a recent flare of interest [9, 10]. This interest is for a good reason because inverting for the normal operator (the demigration-migration operator) restores many of the amplitude artifacts related to acquisition and illumination imprints. However, the downside to this approach is that least-squares tends to smear the energy leading to a loss of resolution. By imposing certain sparseness constraints on the imaged reflectivity, progress has been 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 of noise; (ii) use the optimal (sparse & local) representation properties of Curvelets and their almost diagonalization of the imaging operators; (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], in which 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 appropriate choice for seismic imaging and processing. We proceed by introducing the non-linear estimation procedure with thresholding in the image space. We show that Monte-Carlo sampling techniques can be used to compute a correction for the threshold that incorporates the coloring of the noise due to migration. We conclude by introducing a constrained optimization approach aimed to (approximately, via the diagonal/symbol of the normal operator in the Curvelet domain) invert for the normal operator while imposing 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 migration operator.  The seismic imaging problem In all generality, the seismic imaging problem can after linearization and high-frequency approximation be cast into the following form for the forward model describing our data d = Km + n, (1) where K is the (Kirchoff) (de)-migration/scattering operator, m the model with the reflectivity and n white Gaussian noise. The pertaining inverse problem has the following general form [see e.g. 11, 14] m ˆ :  1 min d − Km m 2  2 2  + µ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 find the appropriate domain to solve this inverse problem? By hitting the data with the migration operator (the adjoint of the scattering operator denoted by ∗ ), followed by sandwiching the normal operator between basis-function (de)-compositions, collapses the energy onto a limited number of coefficients. Question is how to recover these coefficients from the now colored noise, ΨDO  FIO ∗  ∗  Bu = B K K B Bm + B K∗ n  ˜m or u ˜=A ˜ +n ˜.  (3)  At this point, the cruciality of the right choice for the basis functions becomes apparent. Not only do we wish to represent the model  200  0  400  200  0  400 200  0  400 −2  10  demigrated migrated normal  −4  0.5  10  0.5  −6  10  50  1.0  −8  10  1.0  −10  10  −12  10  1.5  −14  10  1.5  100 −16  10  −18  10  2.0  A Curvelet  0  10  2.0  Demigrated Curvelet  1  2  10  10  3  10  4  10  5  10  6  10  Curvelet in FK-domain  Figure 1:  Properties of the Curvelet Transform under imaging. Curvelet and its demigration (first 2 panels). Clearly, the Curvelet remains Curvelet-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 approximation rate for the model and an almost diagonalization of the (de)-migration and normal operators. Mathematically these operators correspond to Fourier Integral (FIO))and Pseudo Differential Operators (ΨDO), respectively. Only under these conditions, we can hope for a significant improvement in the image quality by solving m ˆ :  1 ˜m min u ˜−A ˜ m 2  2 2  + µJ(m) or approximately  m ˆ :  1 ˜ m min u ˜ − diag{A} ˜ m 2  2 2  + µJ(m).  (4)  Curvelets Curvelets as proposed by [2, 4], constitute a relatively new family of non-separable wavelet bases that are designed to effectively represent seismic data with reflectors that generally lie on piece-wise smooth curves (m ∈ C 2 ). Since it has also been shown that these basis functions are well behaved under the action of the (de)-migration and normal operators [1] Curvelets - as opposed to Wavelets - are the right choice. The well-behavedness corresponds to the property that these basis function remain approximately invariant under the operators as we can see from Fig. 1, where not only Curvelets remain Curvelet-like but where also a rapid decay for the coefficients of images of Curvelets are observed. Evidently, these properties make Curvelets suitable to be used in imaging since they not only represent the model very well (they obtain near optimal non-linear approximation rates), but they also almost diagonalize the normal operator and hence the Covariance of to the colored noise on the image, Covn˜ n˜ = E [˜ nn ˜ ∗ ]. With the appropriate correction for the threshold, we are in the position to relatively easily separate noise from model. How do Curvelets obtain such a high non-linear approximation rate? Without being all inclusive [see for details 2, 4, 5], the answer to this question lies 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 ∝ length2 . 1 . directional selective with # orientations ∝ √scale  • local both in (x, y)) and (kx , ky ). • almost orthogonal, they are tight frames.  Migration denoising by thresholding After applying migration, the estimate for the u (denoted by ) involves the solution of a denoising problem in the presence of colored Gaussian noise. For white Gaussian noise, thresholding on the coefficients solves inverse problems of the type given in Eq. 2 by setting K = I. When using Curvelets, u ˆ = B−1 θµ (Bu) (5) solves the denoising problem for piece-wise smooth reflectors. In this expression, θ is the thresholding operator with threshold µ = σ 2 loge N , N number of samples and σ the standard deviation of the noise [see for detail 6, 11, 3]. So how can we correct for the coloring of the noise knowing that the Covariance is almost diagonal? The answer is simple, simply weight the threshold with the square root of the diagonal, i.e. √ µ = η Γ with Γ = diag {Covn˜ n˜ } (6) ans η an additional control parameter (de)-emphasizing the thresholding, setting the confidence interval (e.g. η = 3 corresponds ˆ to a 95 % confidence interval). Applying this threshold, u ˜ = θµ (˜ u) to the coefficients of the image yields, after inverse-Curvelet transforming, the estimate for the migrated image. Comparison with the original noisy image in Fig. 2 shows a remarkable reduction ˆ ˜ of the noise. However, the amplitudes have not been recovered accurately because u ˆ still contains the normal operator, u ˜ ≈ Am.  Amplitude restoration by constraint optimization Now that we have successfully removed the noise from the migrated imaged, it is time to restore the amplitudes and impose the sparseness constraint. With the latter, we hope to (i) remove the imaging [by approximately inverting for the normal operator by  the pseudo-inverse of the diagonal Γ also proposed in 7], estimation and side-band artifacts; (ii) enhance the continuity along the imaged reflectors. To accomplish these goals, we propose the following strategy, which is close to the ideas presented by [4] except that we include imaging operators. Solve the following constrained optimization problem ˆ : m  min J(m) s.t. m  ˆ ˜m |A ˜ −u ˜| ≤ µ  or approximately  ˆ : m  min J(m) s.t. m  ˆ |Γm ˜ −u ˜ | ≤ µ,  (7)  ˜ ≈ Γ. This property shows where, in the approximate case, explicit use has been made of the approximate property that A that Γ acts as the symbol of the normal operator in the Curvelet domain. We solve this constrained optimization by using an augmented Lagrangian method involving a Steepest Decent Method [see for details 12]. In the approximate case, there is besides the computational benefits, the additional advantage that access to the modelling/demigration operator is no longer necessary. In fact, when provided with Γ, we are able to approximately turn a vanilla migration into a least-squares migration. Moreover, the fudge factor on estimated coefficients allows for the minimization of the additional penalty function. Remains how to get the Γ. In case K is a fixed (read non-velocity model dependent) homogeneous operator, such as the Radon transform, the Γ can be calculated analytically [as shown by 6, 3] and some of these results could potentially carry over to seismic imaging. However, acquisition and illumination imprints remain a problem and we therefore use [see also 7, 8] a Monte-Carlo sampling technique to estimate the diagonal of the Covariance operator from diag{Γ} =  1 M  M  n ˜ 2i  (8)  i=1  with n ˜ i Curvelet-transformed migrated images of independent realizations of pseudo-random white Gaussian noise. Typically, about 50 realizations are taken.  Example To illustrate our proposed method, we included a synthetic example of a common-offset Kirchoff migration for a constant velocity model. To show the capabilities of our method, we present our results broadband. We choose 512 sources, time samples, depth and lateral positions. To remove artifacts and improve the continuity along reflectors we used a L1 -norm, J(m) = |m|1 . As we can see from the example shown in Fig. 2, thresholding followed by optimization yields a substantially improved image, where most of the coherent energy has been removed. Conversely, least-squares migration by Conjugate gradients concentrates the energy of the noise towards regions that are not well insonified and fails to restore the amplitudes. Our method, on the other hand, removes most of the artifacts and restores the amplitudes for the larger angles.  Conclusions and discussion The methodology presented is an example of divide and conquer. First, we got a reasonable estimate for the image by thresholding then we dealt with amplitude restoration and artifact removal. We build on the premise that one stands a much better chance to solve an inverse problem by projecting data into the appropriate domain. The combination of migration with the Curvelet transform provides such a domain, where the energy is collapsed onto a limited number of coefficients that correspond to coherent features along 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 the bulk of the coherent noise while preserving the reflectors. Imaging amplitudes and continuity were restored by approximately inverting the normal operator with its diagonal. Thresholded coefficients were used to start the constrained optimization which greatly improved the image quality by imposing the L1 -sparseness constraint. Results for the approximate inversion of the normal operator were slightly better, which can be understood because the diagonal acts as the Curvelet-equivalent for the symbol of the normal operator. Artifacts generated by the imaging operator can correspond to a non-ΨDO behavior, which may be thought off as off-diagonal erroneous components. Current research, is focusing on extension to pre-stack imaging; inclusion of different norms and the removal of other coherent noise-sources (multiples etc.) on the image space and the construction of imaging operators directly in the Curvelet domain.  Acknowledgments The authors would like to thank Emmanuel Cand´es and David Donoho for making their Digital Curvelet Transforms via Unequallyspaced 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.  References [1] E. J. Cand`es and L. Demanet. Curvelets and fourier integral operators. 2002. to appear Comptes-Rendus de l’Academie des Sciences, Paris, Serie I. [2] E. J. Cand`es 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`es and D. L. Donoho. Recovering Edges in Ill-posed Problems: Optimality of Curvelet Frames. Ann. Statist., 30:784–842, 2000. [4] E. J. Cand`es and F. Guo. New multiscale transforms, minimum total variation synthesis: Applications to edge-preserving image reconstruction. Signal Processing, pages 1519–1543, 2002.  0  200  400  0.5  0.5  1.0  1.0  1.5  1.5  2.0  2.0  Noisy Image 0  200  200  0  400  Least-squares migrated Image 400  0.5  0.5  1.0  1.0  1.5  1.5  2.0  2.0  Denoised after Thresholding  200  0  400  Constrained Optimization  Figure 2: Synthetic common-offset Kirchoff-migration example. Top row: Noisy migrated and least-square migrated image with 10 interations of CG. Notice the coloring of the noise and the failure of CG to correct for the normal operator. 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. Harmonic Analysis, 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:


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