UBC Faculty Research and Publications

Seismic deconvolution revisited with curvelet frames Hennenfent, Gilles; Herrmann, Felix J.; Neelamani, Ramesh 2005

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

Item Metadata


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

Full Text

4255  SEISMIC DECONVOLUTION REVISITED WITH CURVELET FRAMES 1  GILLES HENNENFENT1 , FELIX HERRMANN1 and RAMESH NEELAMANI2 University of British Columbia, Department of Earth and Ocean Sciences, 6339 Stores Road, Vancouver, BC V6T 1Z4, Canada 2 ExxonMobil Upstream Research Company, 3319 Mercer St., Houston, TX 77027-6019, USA  Summary We propose an efficient iterative curvelet-regularized deconvolution algorithm that exploits continuity along reflectors in seismic images. Curvelets are a new multiscale transform that provides sparse representations for images (such as seismic images) that comprise smooth objects separated by piece-wise smooth discontinuities. Our technique combines conjugate gradient-based convolution operator inversion with noise regularization that is performed using non-linear curvelet coefficient shrinkage (thresholding). The shrinkage operation leverages the sparsity of curvelets representations. Simulations demonstrate that our algorithm provides improved resolution compared to the traditional Wiener-based deconvolution approach. Introduction In this paper, we address the classical discrete-time deconvolution problem. The forward problem is d = Km + n,  (1)  where d is the data, K the convolution operator (cyclic matrix), m the seismic reflectivity and n zero-mean additive white Gaussian noise with variance σ 2 . The assumption of white noise can be relaxed to colored noise as long as its covariance matrix is near diagonal in curvelet frames. The convolution kernel of (1) is defined as N −1  Gn−i mi + ni .  di =  (2)  n=0  The inverse problem is to determine m given d and K. Notice that we do not consider the problems of blind deconvolution and source signature estimation. In the Fourier domain, we have ˆ ˆ d(ω) = G(ω) m(ω) ˆ +n ˆ (ω).  (3)  A naive approach to obtain an estimate for m ˜ would be m ˜ = Ad = m + An,  (4)  where A· = K−1 · or A· = K† ·, the pseudo-inverse for a non-invertible K. Unfortunately, the discrete deconvolution problem is ill-conditioned as is evidenced by the size of the conditioned number of K cond(K) =  ˆ max G , ˆ min G  (5)  ˆ min is the smallest value of the magnitude of the Fourier transform and G ˆ max the biggest. In that case, the where G variance of An is large making m ˜ an unsatisfactory deconvolution estimate. Since the 60’s, deconvolution filters based on the Wiener theory are widely used. Wiener filters minimize the meansquared-error (MSE) [14]. As early as the mid-70’s to early 80’s, this minimal energy approach was generalized towards an 1 -norm minimization, designed to bring the spikyness of the reflectivity [3, 9, 10]. In the image processing EAGE 67th Conference & Exhibition — Madrid, Spain, 13 - 16 June 2005  and wavelet community, similar ideas emerged in the early 90’s where deconvolution problems are solved by also exploiting sparseness of the model but now on certain bases that are designed to be sparse on the model (see e.g. [13, 8]). Our approach borrows from both the minimal structure concepts introduced in geophysics and from above recent ideas in image processing, where the sparseness of certain frame expansions is used to regularize deconvolution [13, 4]. Since the Earth can be considered to consist of reflectors on piece-wise smooth curves, recently developed curvelet frames are the appropriate choice. We begin by providing a brief overview of curvelet frames. We then present our curvelet-regularized deconvolution algorithm and discuss its practical implementation. An illustrative synthetic example follows. Curvelet Frames Curvelet properties The curvelet transform is a new member in the family of Computational Harmonic Analysis tools [2, 1]. Tight curvelet frames were initially designed to provide optimally sparse representations for objects that are smooth except along smooth curve-like discontinuities (e.g. image with edges). Curvelets are obtained by partitioning the 2D Fourier plane into dyadic coronae and sub-partitioning those into angular wedges. Curvelets obeys a parabolic scaling law – at scale 2−j , each element has an effective support of length 2−j/2 and width 2−j . The resulting curvelet frames are multiscale, multi-directional, highly anisotropic, and localized both in space and frequency. Curvelets can perform decomposition and reconstruction of any function much like an expansion in an orthonormal basis, but in contrast, curvelet representations are moderately redundant. On the practical side, there exists fast O(N log N ) algorithms that allow for a decomposition of n-by-m images f with n ∗ m = N . These algorithms are numerically tight and have an explicit construction for the adjoint that equals the pseudo-inverse C∗ · = C† · and C∗ C = I,  Figure 1: Five curvelets at different scales and orientations. (6)  where C represents the curvelet transform. As a reminder, in an inner product space, say H, the adjoint operator T∗ of a given operator T is defined by f , Tg = T∗ f , g ∀ f , g ∈ H. (7) Hence, f can be written as f , ϕγ ϕγ = C† Cf . def  f=  (8)  γ ∈Γ Curvelet Sparse Representations Consider a function f of two variables that is piece-wise twice differentiable. Further, let the discontinuity curves separating the smooth pieces of f also be piece-wise twice differentiable. From the non-linear perspective, the optimal approximation rate equals [5] 2 −2 f − fO for k → ∞ (9) k 2 ∝k where f O k is the partial reconstruction of f using the k largest terms in the basis. Curvelet frames achieve [2] f − fC k  2 2  ≤ C k −2 (log k)3 for k → ∞  (10) 2  In words, we quote from [2] “there is no basis in which coefficients of an object with an arbitrary C singularity would decay faster than in a curvelet frame”. For comparison, the decay rate for Fourier coefficients is O(k −1/2 ) and for wavelet coefficients O(k −1 ). Curvelet shrinkage-based signal estimation Consider (1) with K = I (i.e. denoising problem). In a basis of wavelets, a simple shrinkage (i.e. soft-thresholding) of the wavelet coefficients of the noisy data minimizes a quadratic distance to the data penalized by a 1 -norm [7]. The variational problem is 1 (11) m ˜ = arg min d − m 22 + λ||Wm||1 , m 2 where W represents the wavelet transform. In curvelet frames, soft-thresholding solves only approximately (11) but already provides a good estimate for the original signal [12].  Curvelets for seismic deconvolution Main features in seismic images correspond to unconformities in the Earth’s subsurface which tend to follow piecewise smooth curves, allowing for pinch-outs and faults. Hence the curvelet representation of seismic images can be expected to be sparse. We use this additional piece of information to fill the null space of the operator and make the inversion more stable. The forward problem (1) is reformulated as d = Fx + n, ∗  (12) ∗  where x is a vector of curvelet coefficients such that m = C x, and F = KC . The variational problem now becomes [4] 1 x ˜ = arg min d − Fx 22 + λ||x||1 . (13) x 2 In words, we want to deconvolve the data and stabilize the process by imposing a sparsity constraint on the curvelet representation of the model. Intuitively speaking, (13) solves the deconvolution problem ”curvelet-wise” but the exact solution is not trivial due to the non-differentiability of the 1 -norm at the origin. A popular choice is the iterative reweighted least-squares (IRLS) method that is known to give a good approximation to the 1 -norm [6, 11] to solve (13). We propose a method inspired by the work in [4]. We combine the Conjugate Gradient (CG) method with soft thresholding to obtain an approximate solution (13). First, the CG-method is applied to the normal equations F∗ d = F∗ Fx + F∗ n. Typically, the noise that gets amplified by the small singular values of F∗ F is regularized by stopping CG prematurely; the stopping criterion is proportional to the noise level. Our approach regularizes the noise by soft thresholding the curvelet coefficients after a limited number of CG iterations and then restarting CG with the thresholded estimate. Besides removing the noise, soft thresholding enhances the spikyness (minimizes the 1 -norm of the coefficients). This iterative CG regularization approximately solves (13) while preserving the frequency information associated with the small singular values. The threshold level is taken large at the beginning and made smaller towards the end. The algorithm is given below: 1. Initialize threshold Γ0 , first guess x0 , set number of CG iterations between thresholdings J, and number of iterations Lmax . 2. For i = 1 : Lmax • x ¯i = CG(˜ xi−1 , J), • Soft-threshold the coefficients x ¯i with the threshold Γi and obtain x ˜i . The reflectivity estimate is given by the inverse transform m ˜ = C∗ x ˜. Results The noisy data are obtained by convolving the Ricker wavelet with the Marmousi seismic reflectivity and adding white Gaussian noise (SNR ≈ 14 dB). There are 512 traces with 512 samples each. As a reference, Fig. 2(a) shows the seismic reflectivity estimate obtained using the Wiener filter, i.e. ˆ H(ω) =  ˆ G(ω) 2 + Pn (ω) ˆ |G(ω)| Pm (ω)  (14)  ˆ where G(ω) is the Fourier transform of the seismic source signature, Pm (ω) the power spectrum of the model, and Pn (ω) the power spectrum of the noise. In practice, Pm (ω) is substituted by Pd (ω) the power spectrum of the data as a first approximation. Fig. 2(b) shows the deconvolution estimate obtained using the proposed algorithm. Unlike the Wiener filter, our method is able to reconstruct frequency components which have been degraded by the noise by exploiting the sparsity of the reflectivity’s curvelet representation. As a result, the frequency content of the deconvolution estimate is improved (Fig. 2(c)). In particular, we point out that the target area is better resolved. Moreover, the ringing effect is also reduced but some artifacts, related to the residual noise and Gibbs-like effects, are still present in the curvelet deconvolved image. Conclusions We developed and demonstrated a new iterative curvelet-regularized deconvolution algorithm that exploits continuity along reflectors in seismic images. The motivation for the curvelet approach stems from the realization that classical methods like Wiener filtering are not able to preserve frequency contents degraded by noise. Sparseness in the curvelet domain, on the other hand, brings out the signal components mostly influenced by the noise. Our algorithm can be applied to a wide range of applications as long as 1) the model is known to be otherwise smooth having discontinuities EAGE 67th Conference & Exhibition — Madrid, Spain, 13 - 16 June 2005  0  50  100  150  200  250  300  350  400  450  500 −1  (a)  (b)  −0.5  0  0.5  1  1.5  2  2.5  3  (c)  Figure 2: (a) Reflectivity estimate using the Wiener filter, (b) reflectivity estimate using iterative curvelet-regularized deconvolution algorithm, and (c) from left to right: trace #360 in the true model, in our estimate, and in the Wiener filter estimate. along piece-wise smooth curves and 2) the covariance matrix of the noise is near diagonal in curvelet frames. So far, our regularization focussed on enhancing the sparseness of the curvelet coefficients. The next step of our research will be to exploit apparent continuity along reflectors which, based on our experience with continuity-enhanced imaging, will further reduce the previously mentioned artifacts. Acknowledgments The authors would like to thank the authors of the Digital Curvelet Transform (Candes, Donoho, Demanet and Ying). 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, SHELL. Additional funding came from the NSERC Discovery Grants 22R81254. References [1] E. Candes and D. Donoho. Curvelets: A surprisingly effective nonadaptive representation of objects with edges. Curves and Surfaces, 1999. L. L. Schumaker et al. (eds), Vanderbilt University Press, Nashville, TN. [2] E. Candes and D. Donoho. New tight frames of curvelets and optimal representations of objects with C 2 singularities. Technical report, Caltech, 2002. [3] J. Claerbout and F. Muir. Robust modeling with erratic data. Geophysics, 38(05):826–844, 1973. [4] I. Daubechies, M. Defrise, and C. de Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constrains. Comm. Pure Appl. Math., pages 1413–1457, 2005. [5] D. Donoho. Sparse components of images and optimal atomic decomposition. Constr. Approx, 17:353–382, 2001. [6] A. Gersztenkorn, J. B. Bednar, and L. R. Lines. Robust iterative inversion for the one-dimensional acoustic wave equation. Geophysics, 51(2), February 1986. [7] S. Mallat. A Wavelet Tour of Signal Processing. Academic Press, 1999. [8] R. Neelamani, H. Choi, and R. Baraniuk. ForWaRD: Fourier-Wavelet Regularized Deconvolution for IllConditioned Systems. IEEE Transactions on Signal Processing, 52(2), February 2004. [9] D. W. Oldenburg, S. Levy, and K. P. Whittall. Wavelet estimation and deconvolution. Geophysics, 46(11):1528– 1542, 1981. [10] M. D. Sacchi, D. R. Velis, and A. H. Cominguez. Minimum entropy deconvolution with frequency-domain constraints. Geophysics, 59(06):938–945, 1994. [11] J. Scales and A. Gersztenkorn. Robust methods in inverse theory. Inverse problems, 4:1071–1091, 1988. [12] J. Starck, E. Candes, and D. Donoho. The curvelet transform for image denoising. IEEE Transactions on Image Processing, 11:670–684, 2000. [13] J.-L. Starck, M. Nguyen, and F. Murtagh. Wavelets and curvelets for image deconvolution: a combined approach. Signal Processing, 83(10):2279–2283, 2003. [14] N. Wiener. The extrapolation, interpolation and smoothing of stationary time series, with engineering applications. Wiley, 1949.  


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