UBC Faculty Research and Publications

Optimization strategies for sparseness- and continuity- enhanced imaging : Theory Herrmann, Felix J.; Moghaddam, Peyman P.; Kirlin, Rodney L. 2005

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

Item Metadata


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

Full Text

4251 OPTIMIZATION STRATEGIES FOR SPARSENESS- AND CONTINUITY- ENHANCED IMAGING: THEORY FELIX J.HERRMANN∗, PEYMAN P. MOGHADDAM AND RODNEY KIRLIN Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, BC, Canada Abstract Two complementary solution strategies to the least-squares migration problem with sparseness- & continuity constraints are pro- posed. The applied formalism explores the sparseness of curvelets on the reflectivity and their invariance under the demigration- migration operator. Sparseness is enhanced by (approximately) minimizing a (weighted) `1-norm on the curvelet coefficients. Continuity along imaged reflectors is brought out by minimizing the anisotropic diffusion or total variation norm which penalizes variations along and in between reflectors. A brief sketch of the theory is provided as well as a number of synthetic examples. Technical details on the implementation of the optimization strategies are deferred to an accompanying paper: implementation Introduction Least-squares migration and migration deconvolution have been topics that received a recent flare of interest [7, 8]. 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 unregularized least-squares tends to fit noise and smear the energy. In addition, artifacts may be created due to imperfections in the model and possible null space of the normal operator [11]. Regularization by minimizing an energy functional on the reflectivity can alleviate some of these problems, but may go at the expense of resolution. Non-linear functionals such as `1 minimization partly deal with the resolution problem but ignore bandwidth-limitation and continuity along the reflectors [12]. Independent of above efforts, attempts have been made to enhance the continuity along imaged reflectors by applying anisotropic diffusion to the image [4]. The beauty of this approach is that it brings out the continuity along the reflectors. However, the way this method is applied now leaves room for improvement regarding (i) the loss of resolution; (ii) the non-integral and non-data constrained aspects, i.e. this method is not constrained by the data which may lead to unnatural results and ’overfiltering’. In this paper, we make a first attempt to bring these techniques together under the umbrella of optimization theory and modern image processing with basis functions such as curvelet frames [3, 2]. Our approach is designed to (i) deal with substantial amounts of noise (SNR ≤ 0); (ii) use the optimal (sparse & local) representation properties of curvelets for reflectivity; (iii) exploit the near diagonalization of the normal operator by curvelets [2]; and (iv) use non-linear estimation, norm minimization and optimization techniques to enhance the continuity along reflectors [5]. Optimization strategies for seismic imaging After linearization the forward model has the following form d = Km+ n, (1) where K is a demigration operator given by the adjoint of the migration operator; m the model wih the reflectivity and n white Gaussian noise with standard deviation σn (colored noise can be accounted for). Irrespective of the type of migration operator (our discussion is independent of the type of migration operator and we allow for post-stack Kirchoff as well as ’wave-equation’ operators), two complementary optimization strategies are being investigated in our group. These strategies are designed to exploit the sparseness & invariance properties of curvelet frames in conjunction with the enhancement of the overall continuity along reflectors. More specifically, the first method [6, 5] preconditions the migration operator, yielding a reformulation of the normal equations into a standard denoising problem F∗d = ≈Idz}|{ F∗Fx+ F∗n (2) y = x+  (3) with ∗ the adjoint; F· = KC∗Γ−1· the curvelet-frame preconditioned migration operator with Γ2· = diag(CK∗KC∗)· ≈ CK∗KC∗· by virtue Theorem 1.1 of [2], which states that Green’s functions are nearly diagonalized by curvelets; C, C∗ the curvelet transform and its transpose; x the preconditioned model related to the reflectivity according m = CΓ−1x and  close to white Gaussian noise (by virtue of the preconditioning). Applying a soft thresholding to Eq. 2 with a threshold proportional to the standard deviation σn of the noise on the data, gives an estimate for the preconditioned model with some sparseness [see for details 01000 2000 3000 z [m ] 0 100 200 300 400 500 x [m] (c) Marmousi model Detailed velocity model x [m] y [m ] 1000 2000 3000 4000 5000 6000 7000 8000 9000 500 1000 1500 2000 2500 3000 Figure 1: Marmousi model and its normal directions. (a) Reflectivity of the Marmousi model; for the detailed velocity model. (b) Directions where the gradients (white ↑’s) is maximum for the smoothed Marmousi model. 5]. For orthogonal bases, soft thresholding [9] solves x̂ = argmin x 1 2 ‖y − x‖22 + ‖x‖1 (4) with ‖c‖1 the `1-norm. Since curvelets are redundant frames and only approximately diagonalize the normal operator, the above thresholding is only approximate which means that x is not optimally sparse. After thresholding, there remains the synthesis of the reflectivity from x̂ which involves the inversion of the preconditioning operator. During this synthesis step, a wavefront-continuity enhancing constraint is imposed on the reflectivity in the form of the following anisotropic norm Ja(m) = ‖Λ1/2∇m‖p (5) with Λ a location-dependent operator (see Fig. 1 for an example) which rotates the action of the gradient towards the approximate normal to the reflector [see 5, where this norm is used]. For p = 2, Eq. 5 corresponds to anisotropic diffusion [4], minimizing L2- norm, and for p = 1 this penalty functional corresponds to anisotropic Total Variation (TV minimizes the `1-norm of the gradient), generalizing isotropic TV. The former penalizes discontinuities in the amplitudes along the reflectors whereas the latter allows for discontinuities along the reflectors. The synthesis step involves the solution of the following constrained optimization problem m̂ = argmin m Ja(m) s.t. |ΓCm− x̂|µ ≤ eµ ∀µ ∈M. (6) In words, this optimization procedure tailors each curvelet coefficients with index µ ∈ M (M is the curvelet index set) such that the continuity along the reflectors is enhanced. The tolerance vector e ensures that optimized coefficients do not wander off too far from the coefficients yielded by the denoised migrated image. This method has been implemented for Λ fixed and we are currently researching how to make the estimation of the spatial varying rotation matrix part of the optimization. Advantages of the above divide-and-conquer approach are (i) the clearly identifiable steps of first denoising and then reconstructing with continuity enhancement; (ii) flexible constraints that connect the continuity enhancement with the denoised migrated data. Disadvantages are (i) the increased problem size (curvelet frames are redundant); (ii) the approximate nature of the thresholding procedure; (iii) the necessity to compute the Γ, which can be generalized to non-diagonal Lanczos approximations [see the ensuing paper and 10, 5]; (iv) required knowledge of the tolerances e that depend on the noise level of the migrated image and (v) the computational expense in solving the constrained optimization problem. The second optimization strategy aims to invert Eq. 1 using m̂ = argmin m J(m) s.t. ‖d−Km‖2 = Nσn (7) with N the number of data samples. As opposed to the first optimization strategy where a divide-and-conquer approach was followed, we propose to solve Eq. 7 by inverting the scattering operator while simultaneously minimizing a weighted sparseness constraint Jc(m) = ‖Cm‖1,w (8) on the curvelet coefficients and the continuity constraint of Eq. 5 on the model with J(m) = αJc(m) + Ja(m) and α+ β = 1. (9) The positive weighting vector w is designed to penalize the low and high frequency ends since these are not present in the data. Since curvelets are multiscale, this weighting simply corresponds to increasing the weight for the coarsest and finest scales. The 01000 2000 3000 z [m ] 0 2000 4000 6000 8000 x [m] (d) Noisy least-squares image 0 1000 2000 3000 z [m ] 0 2000 4000 6000 8000 x [m] (c) Denoised with detailed prior Figure 2: Conventional non-regularized least-squares imaging versus optimized imaging using the first strategy. (a) Least-squares migrated image; (b) Optimized image using the penalty functional using information on the normal directions depicted in Fig. 1. Note the striking improvement. coefficients α and β weight the relative importance of the sparseness on the curvelet coefficients versus enhancement of continuity along the reflectors. Details on the implementation of this optimization are given in the ensuing paper that deals with the implemen- tation and practical considerations. Advantages of this second method are (i) the reduced model size; (ii) the integrated approach; (iii) availability of fast solvers, which use the approximate inverse of the normal operator. Disadvantages are the more complicated penalty functional and the computation of approximate inverse of the Hessian/normal operator. Examples Our two optimization strategies are illustrated by an example for the Marmousi model. Fig. 1 shows the normal directions (white ↑’s) computed from the smoothed velocity model. In practice, similar information on the smoothed velocity model is used to compute the migration operator. Without loss of generality, we use a post-stack Kirchoff modeling-migration operator for a constant velocity model. Synthetic data is generated according to 1 for a single geophone and 1800 sources. The scattering operator K is of the common-offset Kirchoff type for a constant velocity model with c̄ = 3500ms−1. The reflectivity model is given by the derivative of the Marmousi model [1] (see Fig. 1) (a). The image size is 512 × 512 with a depth-sample interval of 5.8m and a horizontal-sample interval of 18.0m. The source locations range from s ∈ [−10800, 13540]m and the geophone record 512 time samples with a sample interval ∆t = 2ms. Both source and receiver signatures are ignored, yielding a broad-band response. Noise was added yielding a signal-to-noise ratio of SNR = 0dB, which corresponds to σ = 1 for a unit signal standard deviation. Fig’s 2-3, show examples of both optimization strategies for the above synthetic experiment. Both strategies greatly improve the image quality. The results for the second strategy are better which can be explained by the more rapid convergence of the Newton method underlying this strategy. For details refer to the ensuing paper. Discussion and conclusions We presented a new imaging approach which aims to employ a delegate balance between sparseness in the curvelet domain on the one hand and continuity in the physical domain on the other. We presented two alternative strategies both of which conduct a `1-norm minimization on the curvelet coefficients in combination with the minimization of an anisotropic diffusion/Total-Variation norm on the reflectivity but differ with respect to the employed optimization techniques. The first strategy is built on a convex optimization technique enforcing the inequality in Eq. 6, whereas the second strategy entails a non-convex optimization procedure. Irrespective of the optimization strategies, explicit use is made that curvelet coefficients are large whenever curvelets align with parts of the reflector. The `1-norm boosts these coefficients and hence brings out the reflectivity at the expense of the noise. The subsequent anisotropic norm minimization enforces lateral consistency by ’chaining” the sparse curvelets together along the reflectors. This approach can be seen as a first step towards defining a framework for non-linear regularization in seismic imaging. Acknowledgments The author would like to thank the authors of the Digital Curvelet Transform (Candes, Donoho, Demanet and Ying). We also would like to thank dr. Sacchi for his migration code. This work was carried out as part of the SINBAD project with financial support, secured through ITF (the Industry Technology Facilitator), from the following organizations: BG Group, BP, ExxonMobil and SHELL. Additional funding came from the NSERC Discovery Grants 22R81254 and 0004049. References [1] A. Bourgeois, M. Bourget, P. Lailly, M. Poulet, P. Ricarte, and R. Versteeg. The Marmousi experience, volume 5-9, chapter Marmousi data and model. EAGE, 1991. x[m] Denoised using integrated non−convex method z[m ] 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 0 500 1000 1500 2000 2500 3000 Figure 3: Results for non-convex optimimization. [2] E. J. Candès and L. Demanet. The curvelet representation of wave propagators is optimally sparse. 2004. Preprint: http://www.acm.caltech.edu/ demanet. [3] 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. [4] G. C. Fehmers and C. Hocker. Fast structural interpretation with structure-oriented filtering. Geophysics, 68(4), 2003. [5] F. Herrmann and P. Moghaddam. Sparseness- and continuity-constrained seismic imaging with curvelet frames. 2005. Submitted for publication. [6] F. J. Herrmann and P. Moghaddam. Curvelet-domain least-squares migration with sparseness constraints. In Expanded Abstracts. EAGE, 2004. [7] J. Hu, G. T. Schuster, and P. Valasek. Poststack migration deconvolution. Geophysics, 66(3):939–952, 2001. [8] H. Kuhl and M. D. Sacchi. Least-squares wave-equation migration for avp/ava inversion. Geophysics, 68(1):262–273, 2003. [9] S. G. Mallat. A wavelet tour of signal processing. Academic Press, 1997. [10] P. Moghaddam and F. J. Herrmann. Migration preconditioning with curvelets. In Expanded Abstracts, Tulsa, 2004. Soc. Expl. Geophys. [11] J. E. Rickett. Illumination-based normalization for wave-equation depth migration. Geophysics, 68, 2003. [12] Y. Wang. Sparseness-constrained least-squares inversion: Application to seismic wave reconstruction. In Geophysics, vol- ume 68, pages 1633–1638. Soc. of Expl. Geophys., 2003.


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