UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Incoherent noise suppression and deconvolution using curvelet-domain sparsity Vishal, Kumar 2009

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

Item Metadata

Download

Media
24-ubc_2009_fall_kumar_vishal.pdf [ 9.81MB ]
Metadata
JSON: 24-1.0052886.json
JSON-LD: 24-1.0052886-ld.json
RDF/XML (Pretty): 24-1.0052886-rdf.xml
RDF/JSON: 24-1.0052886-rdf.json
Turtle: 24-1.0052886-turtle.txt
N-Triples: 24-1.0052886-rdf-ntriples.txt
Original Record: 24-1.0052886-source.json
Full Text
24-1.0052886-fulltext.txt
Citation
24-1.0052886.ris

Full Text

Incoherent noise suppression and deconvolution using curvelet-domain sparsity by Vishal Kumar B.Sc., Indian Institute of Technology, Kharagpur, 2006 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Geophysics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June, 2009 c Vishal Kumar 2009  Abstract Curvelets are a recently introduced transform domain that belongs to a family of multiscale and also multidirectional data expansions. As such, curvelets can be applied to resolution of the issues of complicated seismic wavefronts. We make use of this multiscale, multidirectional and hence sparsifying ability of the curvelet transform to suppress incoherent noise from crustal data where the signal-to-noise ratio is low and to develop an improved deconvolution procedure. Incoherent noise present in seismic reflection data corrupts the quality of the signal and can often lead to misinterpretation. The curvelet domain lends itself particularly well for denoising because coherent seismic energy maps to a relatively small number of significant curvelet coefficients while incoherent energy is spread more or less evenly amongst all curvelet coefficients. Following standard processing of crustal reflection data, we apply our curvelet denoising algorithm to deep reflection data. In terms of enhancing the coherent energy and removing incoherent noise, curvelets perform better than the F-X prediction method. We also use the curvelet transform to exploit the continuity along reflectors for cases in which the assumption of spiky reflectivity may not hold. We show that such type of seismic reflectivity is sparse in the curvelet-domain. This curvelet-domain compression of reflectivity opens new perspectives towards solving classical problems in seismic processing, including the deconvolution problem. We present a formulation that seeks curvelet-domain sparsity for non-spiky reflectivity. Comparing the results with those obtained from sparse spike deconvolution, curvelets perform better than the latter by recovering the frequency components, which get degraded by convolution operator and noise.  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Abstract  List of Tables  Preface  Acknowledgments Dedication  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii  Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . 1.1 Theme . . . . . . . . . . . . . . . . . . . . . . . 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . 1.4 Theoretical background . . . . . . . . . . . . . 1.4.1 Sparsity-promoting inversion background 1.4.2 Sparsifying transforms . . . . . . . . . . 1.4.3 Sparse representation of seismic data . . 1.5 Outline . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  1 1 2 2 2 2 4 4 11  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  12  2 Incoherent noise suppression with curvelet-domain sparsitypromotion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Elementwise hard and soft thresholding . . . . . . . .  13 13 14 14 iii  Table of Contents . . . . . .  15 16 17 22 22 24  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  29  2.3 2.4 2.5 2.6 2.7  2.2.2 One-norm . . . . . . . . . . . Parameter selection for synthetic data Application to synthetic data . . . . . Parameter selection for real data . . . Application to deep reflection data . . Results and discussion . . . . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . .  41 41 41 41 42  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  44  . . . . .  . . . . .  . . . . .  . . . . .  . . . . . .  4 Conclusions . . . . . . . . . . . . . . . . . . . 4.1 Curvelet-domain sparsity . . . . . . . . . . 4.2 Application to incoherent noise suppression 4.3 Application to deconvolution . . . . . . . . 4.4 Open and future research . . . . . . . . . .  . . . . .  . . . . .  . . . . . .  40  . . . . .  . . . . .  . . . . . .  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . .  . . . . . .  31 31 32 33 33  . . . . .  . . . . .  . . . . . .  . . . . .  3 Deconvolution with curvelet-domain sparsity 3.1 Background . . . . . . . . . . . . . . . . . . . 3.2 Problem formulation . . . . . . . . . . . . . . . 3.3 Parameter selection . . . . . . . . . . . . . . . 3.4 Application and results . . . . . . . . . . . . .  . . . . .  . . . . . .  . . . . .  . . . . .  . . . . .  iv  List of Tables 2.1  2.2 3.1  3.2  The cooling method with iterative soft thresholding to solve Eq. 3.5. L is the number of inner iterations, K is the number of outer iterations, k is the counter for outer iterations, λ1 · · · λK are the values of λ in decreasing order, other symbols are as explained in the text. . . . . . . . . . . . . . . . . SNR comparison for pre-stack denoising . . . . . . . . . . . .  17 22  Pseudo code to solve the optimization problem shown in Eq. 4.4. K is maximum number of iterations, k is the counter for iterations, other quantities are as stated in the text. . . . . . . 33 SNR comparison for deconvolution methods . . . . . . . . . . 35  v  List of Figures 1.1  1.2  1.3  (a) A few curvelets in t-x, (b) Same curvelets in F-K domain. The colored arrows show the positions of some curvelets in time-space and their representation in the frequency-wavenumber domain. Each curvelet function is spatially localized because its amplitude rapidly decays to zero outside a certain region. The curvelets are localized in a wedge in the F-K domain. The orientation of a curvelet in the physical domain is perpendicular to its wedge in the F-K domain. (c) Curvelet (top) and F-K (bottom) decomposition of two dipping events (conflicting dips). Note the localization of curvelet functions and global nature of F-K functions used to decompose the image. 6 (Adapted from Neelamani et al. (2008)) . . . . . . . . . . . . Discrete curvelet partitioning of the 2-D Fourier plane into second dyadic coronae and sub-partitioning of the coronae into angular wedges. The tiling corresponds to 5 scales, 8 angles at the 2nd coarsest level and curvelets at the finest (fifth) scale. The angles double every other scale (Adapted from Hennenfent and Herrmann (2006)). . . . . . . . . . . . . 7 Curvelet decomposition of a shot gather at different frequencyband (scale) and angles (dip). Four scales are used for the decomposition. The centre (coarsest scale) shows the DC and low frequency components of the shot gather. The 2nd coarsest scale has 16 angles. The number of angles is doubled to 32 at the third scale and stays the same (32) for the fourth scale. Note the portions of shot gather captured at various angles and scales. . . . . . . . . . . . . . . . . . . . . . . . . . 8  vi  List of Figures 1.4  1.5  2.1  2.2  Decay of coefficients of synthetic models mentioned in the text. The models are taken to various transform domains and the coefficients are sorted in their descending order. (a) The shot gather, (b) Poststack model, (c) Decay rates of transform coefficients for shot gather, (d) Decay rates of transform coefficients for poststack image. The decay rate of coefficients is proportional to the sparsity of the model in the transform domain. The curvelet coefficients (solid pink) decay faster than wavelet coefficients (dashed red) and Fourier coefficients (dashed blue). Hence curvelets are the most appropriate sparsifying transform for these two models. Since curvelet transform is redundant in nature; i.e., it produces more coefficients than the data size, the x-axis is chosen as percentages of coefficients rather than absolute number of coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Partial reconstruction of the two models (Fig. 2.4(a) and Fig. 2.4(b)) in different transform domains. The model is taken to the transform domain and reconstructed from 1 % of the amplitude-largest coefficients. (a) Fourier reconstruction of shot gather, (b) Fourier reconstruction of poststack image, (c) Wavelet reconstruction of shot gather, (c) Wavelet reconstruction of poststack image, (d) Curvelet reconstruction of shot gather, (e) Curvelet reconstruction of poststack image. As can be observed, the curvelet reconstruction is the most accurate approximation in terms of resemblance of the reconstructed images with the true images. . . . . . . . . . . . . . (a) A seismic signal; dashed lines show the threshold level (λ), (b) Signal after hard thresholding (λ=2), (c) Signal after soft thresholding (λ=2). Both hard and soft thresholding sets entries less than the threshold to zero. In addition, soft thresholding shrinks all remaining entries by the threshold. This moves the signal toward zero. . . . . . . . . . . . . . . . (a) Comparison of signal-to-noise ratio (SNR) vs. threshold value for white noise for hard thresholding and soft thresholding of curvelet coefficients, (b) Comparison of signal-tonoise ratio (SNR) vs. threshold value for colored noise for hard thresholding and soft thresholding of curvelet coefficients. The arrows point to the peak SNR values and the corresponding threshold value (λ) is chosen. . . . . . . . . . .  9  10  15  18 vii  List of Figures 2.3 2.4  2.5  2.6  2.7  2.8  (a) Noise-free data (model), (b) Model with white noise, (c) Model with colored noise. . . . . . . . . . . . . . . . . . . . . Data corrupted with white noise: (a) One-norm approach, (b) Noise removed by one-norm approach, (c) Hard threshold approach, (d) Noise removed by hard threshold approach, (e) Soft threshold approach, (f) Noise removed by soft threshold approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Data corrupted with colored noise: (a) One-norm approach, (b) Noise removed by one-norm approach, (c) Hard threshold approach, (d) Noise removed by hard threshold approach, (e) Soft threshold approach, (f) Noise removed by soft threshold approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tectonic map of the acquisition area and location of SNORCLE near vertical incidence (NVI) reflection Line 1. Sediments of the Phanerozoic Western Canada Sedimentary Basin overlie the Precambrian domains west of the long dashed line. Short dashed black lines show political boundaries. CS Coronation Supergroup, SD - Sleepy Dragon, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC - British Columbia, NWT - Northwest Territories, YK - Yukon. Inset (top-left) shows location of map within Canada. Inset (green box) shows the area for the data used. . . . . . . . . . . . . . . . Shot gather from Line 1 segment highlighted in Fig. 3.6. (a) Noisy data without AGC, (b) Noisy data with AGC, (c) Result of applying curvelet denoising (inversion approach) on the shot gather with AGC, (d) Noise removed by curvelet denoising. (e) The flowchart with important processing steps for pre-stack denoising. The curvelet denoising is applied before the bandpass filtering. . . . . . . . . . . . . . . . . . . . . . . (a) The original stacked seismic section contaminated with random noise. Arrows indicate some of the features discussed in the text. (b) The flowchart with important processing steps for poststack denoising. The curvelet denoising is applied as an alternative to F-X prediction. (c) Result of applying curvelet denoising (inversion approach) on the same section. (d) Noise removed by curvelet denoising, (e) Result of applying F-X prediction on the same section, (f) Noise removed by F-X prediction. . . . . . . . . . . . . . . . . . . . . . . . . .  19  20  21  23  26  27  viii  List of Figures 2.9  3.1  3.2  3.3  3.4  (a) Stacked section after pre-stack curvelet denoising, (b) Stacked section after conventional processing (without F-X prediction), (c) Stacked section after post-stack curvelet denoising. The poststack curvelet denoising is observed to be better than the other two. . . . . . . . . . . . . . . . . . . . . (a) True reflectivity model, (b) Noisy data obtained by convolving the reflectivity model with Ricker wavelet and addition of noise (σ=50), (c) Estimated reflectivity with curvelet deconvolution, (d) Estimated reflectivity with sparse spike deconvolution. The box shows the region of interest that is enlarged in Fig. 4.2. . . . . . . . . . . . . . . . . . . . . . . . (a) Zoom-in plot of true reflectivity model, (b) Zoom-in plot of noisy data, (c) Zoom-in plot of estimated reflectivity with curvelet deconvolution, (d) Zoom-in plot of estimated reflectivity with sparse spike deconvolution. Note the resemblance of curvelet-estimated reflectivity with the true reflectivity. . . (a) True reflectivity model, (b) Noisy data obtained by convolving the reflectivity model with Ricker wavelet and addition of noise (σ=.01), (c) Estimated reflectivity with curvelet deconvolution, (d) Estimated reflectivity with sparse spike deconvolution. The box shows the region of interest (reservoir area) for which enlarged versions of the sections are displayed in Fig. 4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Zoom-in plot of true reflectivity model, (b) Zoom-in plot of noisy data, (c) Zoom-in plot of estimated reflectivity with curvelet deconvolution, (d) Zoom-in plot of estimated reflectivity with sparse spike deconvolution. Note the resemblance of curvelet-estimated reflectivity with the true reflectivity. . .  28  36  37  38  39  ix  Preface The results in this thesis was prepared with Madagascar; a reproducible research software package available at rsf.sf.net and Claritas; a standard package for seismic processing. The reproducibility allows direct access to the experiments performed in this work allowing transparency and dissemination of knowledge. The programs required to reproduce the reported results are Madagascar programs written in C/C++, Matlab R , or Python. The results in this thesis are entirely reproducible with the exception of real data examples in Chapter 2.  x  Acknowledgments I would first like to thank my supervisors Felix Herrmann and Ron Clowes, for their support since my first days at UBC. Without their vision and courage none of this (and much, much more) would have been possible. I would also like to thank the SLIM team (both past and present) for their encouragement, commitment and spirit. It has been an honor to be part of this dynamic group and to witness the amazing process of development that it has undergone. I would like to thank Jounada Oueity for his help and support with our combined research. I am very fortunate to have Michael Bostock as part of my supervisory committee. The knowledge, insight and inspiration that he has provided throughout my time at UBC have been invaluable. I would like to thank Geoscience BC, Society of Exploration Geophysicists and Canadian Society of Exploration Geophysicists for providing annual scholarship. Without the support and love of my parents and the rest of my family, I certainly would never have accomplished any of this. They have my unending thanks and love. Finally, I would like to thank my wife Khushboo Kumari for her support and understanding that enabled me to produce results at my work. All the work for this thesis was financially supported by the NSERC Discovery Grants 22R81254 (FH) and 22R87707(RC) and CRD Grant DNOISE 334810-05 (FH). It was carried out as part of the project with support, secured through ITF, from the following organizations: BG Group, BP, Chevron, ExxonMobil and Shell.  xi  To my family.  xii  Statement of Co-Authorship A version of Chapter 2 will be submitted for publication. I prepared the manuscript with inputs from Jounada Oueity, Ron Clowes and Felix Herrmann. A version of Chapter 3 is a published extended abstract for which the reference is: V. Kumar and F. J. Herrmann, 2008, Deconvolution with curveletdomain sparsity, In Expanded Abstracts, Society of Exploration Geophysicists, Tulsa. The manuscript was jointly written with Felix Herrmann. He led the more theoretical sections, and I led the more applied ones.  xiii  Chapter 1  Introduction Processing of seismic reflection data faces the challenge of noise contamination (Olhovich, 1964) and loss of frequency components (Yilmaz, 2001). The noise contamination and missing frequency in seismic reflection data corrupt the quality of the signal and can often lead to misinterpretation.  1.1  Theme  The main theme of this thesis is a transform-based approach to the seismic denoising and deconvolution problems. Two key features motivate our approach, namely: • High dimensionality: Seismic data are typically five dimensional: time, two spatial coordinates for the sources, and two spatial coordinates for the receivers (for a 3D survey). • Strong geometrical structure: Seismic data are a spatial-temporal sampling of the reflected wavefield that correspond to different interactions of the incident wavefield with inhomogeneities in the Earth’s subsurface. The wavefields have wavefront-like structure with spatial continuity that needs to be exploited. To make the most of these features, our approach uses the curvelet transform (Cand`es and Donoho, 2000), which is data-independent, multiscale, and multidirectional. The basic elements of this transform, the curvelets, are localized in the frequency domain and are of rapid decay in space. Spatial curvelets look like localized plane waves. Because of these properties, curvelets are very efficient at representing curve-like events (e.g., wavefronts, seismic images). The transform concentrates the energy signal amongst a few significant coefficients. In other words, only a few curvelets are needed to represent the complexity of real seismic data. We use this sparsity to solve the problem of noise suppression (also referred to as denoising) and frequency enhancement (also referred to as deconvolution). 1  Chapter 1. Introduction  1.2  Objectives  The objectives of this thesis are twofold: • To utilize the curvelet domain as a sparse representation of seismic data/models that exploits high dimensionality and spatial continuity. • To formulate practical sparsity-promoting algorithms, which address seismic denoising and deconvolution problems, using inversion theory.  1.3  Outline  In Chapter 1, we give a theoretical background regarding the sparsitypromoting inversion approach, various sparsifying transforms and curvelets. In Chapter 2, we present an overview of the curvelet-transform-based denoising algorithms and their application to synthetic data. Three curveletbased methods are tested and the most appropriate (with greater signal-tonoise ratio) is applied to the real deep crustal data. For the real data, results of our curvelet-based algorithm are compared with standard F-X prediction (also known as F-X deconvolution) results. The parameters are carefully selected such that we don’t harm the coherent energy of the signal. Chapter 3 deals with the application of the curvelet-based inversion approach to the deconvolution problem and its application. The two-dimensional structure of reflectivity is exploited by curvelets. Results are compared with standard trace-by-trace sparse spike deconvolution method. Chapter 4 has conclusions and highlights the scope for future work.  1.4  Theoretical background  The fundamental concept behind the methods discussed in this thesis is based on sparsity of seismic signal in a transform domain that facilitates the use of sparsity-promoting methods. For this reason we need to understand the underlying theory behind sparsity-promoting inversion.  1.4.1  Sparsity-promoting inversion background  From the perspective of inversion, observed data can be expressed as: y = Km + n,  (1.1)  2  Chapter 1. Introduction where y is the known data, K is the linearized modeling operator, m is the unknown model and n is additive zero-centered Gaussian noise. The noise can be white (full-band) or colored (band-limited). Given K and y, our objective is to estimate m. In a transform domain, the model can be approximated by the superposition of basic elements of the transform (basis or frames1 ) with appropriate weighting and thus the model can be written as: m= λ i ϕi , (1.2) i∈I  where λi are the weighting of each basic elements ϕi . The transform coefficients λi are inner products of the basic elements and the model at a specific location and orientation. In equation form, this can be represented by: m=  m, ϕi ϕi .  (1.3)  m = S† Sm = S† x,  (1.4)  i∈I  In matrix-vector notation:  where S is the transform, S† is the inverse transform and x is the vector of transform coefficients. The columns of S† are the basic elements (discrete) and λi are the elements of x. For this work, we will seek a transform domain that provides sparse representation of the model; i.e., the elements of x should have few significant values (or fast decay). Assuming, we have a sparsifying transform for the model; we can use transform domain sparsity as a priori information to formulate our inversion problem: ˜ = arg min ||x||1 x x  s.t. ||y − KS† x||2 ≤ ε,  (1.5)  ˜ represents the estimate of the transform coefficient vector x, S† is where x the inverse transform operator and is proportional to the noise level. By solving Eq. 1.5, we find the sparsest set of transform coefficients that explain ˜ = the data within the noise level. The final estimated model is given by m † ˜ . For this work, we are interested in estimating two-dimensional modSx els and for this reason we need to understand the various two-dimensional sparsifying transforms. 1  As opposed to orthonormal basis, redundant frame expansions decompose a length K signal into a frame expansion with M > K elements. Consequently, the composition matrix (S† ) is rectangular with the number of columns exceeding the number of rows.  3  Chapter 1. Introduction  1.4.2  Sparsifying transforms  The Fourier Transform is by far the most important used in seismology. Fourier’s theory states that a given signal can be synthesized as a summation of sinusoidal waves of various amplitudes, frequencies and phases. The limitation of Fourier transform is that it needs a large number of coefficients to account for any jump in the signal and hence not sparse for any signal that is discontinuous. The solution to handling discontinuities (jumps) comes in the form of wavelets that are local in nature and hence can approximate a jump with few wavelet coefficients. Although, the wavelet transform is local (multiscale) in nature, it is poorly directional in higher dimensions. This is because the basis functions of 2D/3D wavelet are isotropic in nature (has no sense of direction). Since seismic wavefields contain wavefronts moving in all directions, wavelets require a larger number of elements to represent the seismic data. More recently, curvelets have been introduced to resolve the issues of directionality of complicated wavefronts as they belong to a family of multiscale, and also multidirectional data expansions. Curvelets obey a parabolic scaling principle with length2 width, with the frequency of the oscillations across the curvelet increasing with increasing scales. Fig. 1.1 shows a few curvelets in the physical and Fourier domains. The comparison of curvelet and F-K decomposition of two dipping events is also shown in Fig 1.1. The curvelets differ from F-K in the sense that the curvelets are able to decompose conflicting dips into localized functions of different frequencywavenumber bands, while F-K functions are global in nature. The increasingly anisotropic nature of curvelets can be attributed to the tiling of the frequency plane into a collection of angular wedges positioned at different scales with the number of angular wedges doubling at every second scale as seen in Fig. 1.2. These wedges in the frequency plane are strictly localized while the curvelets in the physical domain have a rapid spatial decay. Fig. 1.3 shows the decomposition of a shot gather at different scales and angles in the curvelet domain (also referred to as a mosaic plot).  1.4.3  Sparse representation of seismic data  A transform’s ability to efficiently represent a seismic signal can be judged by examining the decay curve of magnitude-sorted coefficients. In other words, the sparsifying transform is able to concentrate most of the signal information within a few large coefficients. A decay curve with a rapid decay to near zero will represent a more parsimonious representation of the seismic 4  Chapter 1. Introduction image as the transform concentrates most of the signal information within few large coefficients. Fig. 1.4 shows the decay of transform coefficients for a synthetic shot gather and a synthetic poststack model. The shot gather (Fig. 1.4(a)) is taken from the work of Herrmann et al. (2007). The poststack model (Fig. 1.4(b)) is a subset of the Marmousi model (Versteeg, 2008), which is a function of depth and offset. For the purpose of our experiments, we assume that the depth axis is the time axis with sampling rate of 4 ms. As we see from Fig. 1.4, the curvelet coefficients have the fastest decay of magnitudes and hence the two models are relatively sparse in the curvelet domain. Another way to measure sparsity is by partially reconstructing a signal with few significant transform coefficients. The signal is taken into the transform domain, the transform coefficients are sorted in descending order of their amplitudes and then 1 % of the largest amplitude sorted coefficients are kept and the remainder are discarded (put to zero). The transform coefficients are then inverse transformed to the physical domain and the reconstructed image is observed. The transform that closely approximates the signal in the process of reconstruction is considered to be best-suited sparsifying transform. Fig. 1.5 shows the partial reconstruction of two models in Fourier, wavelet and curvelet domain. As we observe, the curvelet reconstruction shows a better approximation of the initial test models. Thus, both the test models are sparse in curvelet domain and provides motivation for the development of effective sparsity-promoting algorithms in subsequent chapters of this thesis. Considering the above-mentioned examples, we choose curvelets as our sparsifying transform. For this work, we choose the numerically tight FDCT via wrapping as our curvelet transform (Candes et al., 2005). For this transform, the pseudo-inverse (denoted by the symbol †) equals the adjoint and we have m = C† x = CH x implying CH C = I . Hence our inversion problem (Eq. 1.5) simplifies to: P :  ˜ = arg minx ||x||1 x ˜ = CH x ˜. m  s.t. ||y − KCH x||2 ≤ ε,  (1.6)  For the sake of simple notation, we will use A = KCH in the subsequent chapters of this thesis. Thus for a given modeling operator K, we can formulate our problem in the preceding manner for promoting curvelet-domain sparsity of models. The selection of parameters will be discussed in detail in subsequent chapters of this thesis. 5  Chapter 1. Introduction  (b)  (a)  (c)  Figure 1.1: (a) A few curvelets in t-x, (b) Same curvelets in F-K domain. The colored arrows show the positions of some curvelets in time-space and their representation in the frequency-wavenumber domain. Each curvelet function is spatially localized because its amplitude rapidly decays to zero outside a certain region. The curvelets are localized in a wedge in the F-K domain. The orientation of a curvelet in the physical domain is perpendicular to its wedge in the F-K domain. (c) Curvelet (top) and F-K (bottom) decomposition of two dipping events (conflicting dips). Note the localization of curvelet functions and global nature of F-K functions used to decompose the image. (Adapted from Neelamani et al. (2008))  6  Chapter 1. Introduction  k2 2j/2  angular wedge  2j  k1  Figure 1.2: Discrete curvelet partitioning of the 2-D Fourier plane into second dyadic coronae and sub-partitioning of the coronae into angular wedges. The tiling corresponds to 5 scales, 8 angles at the 2nd coarsest level and curvelets at the finest (fifth) scale. The angles double every other scale (Adapted from Hennenfent and Herrmann (2006)).  7  Chapter 1. Introduction  Figure 1.3: Curvelet decomposition of a shot gather at different frequencyband (scale) and angles (dip). Four scales are used for the decomposition. The centre (coarsest scale) shows the DC and low frequency components of the shot gather. The 2nd coarsest scale has 16 angles. The number of angles is doubled to 32 at the third scale and stays the same (32) for the fourth scale. Note the portions of shot gather captured at various angles and scales.  8  Chapter 1. Introduction  (a)  (b)  (c)  (d)  Figure 1.4: Decay of coefficients of synthetic models mentioned in the text. The models are taken to various transform domains and the coefficients are sorted in their descending order. (a) The shot gather, (b) Poststack model, (c) Decay rates of transform coefficients for shot gather, (d) Decay rates of transform coefficients for poststack image. The decay rate of coefficients is proportional to the sparsity of the model in the transform domain. The curvelet coefficients (solid pink) decay faster than wavelet coefficients (dashed red) and Fourier coefficients (dashed blue). Hence curvelets are the most appropriate sparsifying transform for these two models. Since curvelet transform is redundant in nature; i.e., it produces more coefficients than the data size, the x-axis is chosen as percentages of coefficients rather than absolute number of coefficients.  9  Chapter 1. Introduction  (a)  (b)  (c)  (d)  (e)  (f)  Figure 1.5: Partial reconstruction of the two models (Fig. 2.4(a) and Fig. 2.4(b)) in different transform domains. The model is taken to the transform domain and reconstructed from 1 % of the amplitude-largest coefficients. (a) Fourier reconstruction of shot gather, (b) Fourier reconstruction of poststack image, (c) Wavelet reconstruction of shot gather, (c) Wavelet reconstruction of poststack image, (d) Curvelet reconstruction of shot gather, (e) Curvelet reconstruction of poststack image. As can be observed, the curvelet reconstruction is the most accurate approximation in terms of resemblance of the reconstructed images with the true images. 10  Chapter 1. Introduction  1.5  Outline  In Chapter 2, we present an overview of the curvelet-transform-based denoising algorithms and their application to synthetic data. Three curveletbased methods are tested and the most appropriate (with greater signal-tonoise ratio) is applied to the real deep crustal data. For the real data, results of our curvelet-based algorithm are compared with standard F-X prediction (also known as F-X deconvolution) results. The parameters are carefully selected such that we don’t harm the coherent energy of the signal. Chapter 3 deals with the application of the curvelet-based inversion approach to the deconvolution problem and its application. The two-dimensional structure of reflectivity is exploited by curvelets. Results are compared with standard trace-by-trace sparse spike deconvolution method. Chapter 4 has conclusions and highlights the scope for future work.  11  Bibliography Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899. Cand`es, E. J. and D. L. Donoho, 2000, Curvelets – a surprisingly effective nonadaptive representation for objects with edges. Curves and Surfaces. Vanderbilt University Press. Hennenfent, G. and F. J. Herrmann, 2006, Seismic denoising with nonuniformly sampled curvelets: Computing in Science and Engineering, 8, 16–25. Herrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Nonlinear primary-multiple separation with directional curvelet frames: Geophysical Journal International, 170, 781–799. Neelamani, R., A. I. Baumstein, D. G. Gillard, M. T. Hadidi, and W. L. Soroka, 2008, Coherent and random noise attenuation using the curvelet transform: The Leading Edge, 27, 240–248. Olhovich, V. A., 1964, The causes of noise in seismic reflection and refraction work: Geophysics, 29, 1015–1030. Versteeg, R., 2008, The marmousi experience: Velocity model determination on a synthetic complex data set: The Leading Edge, 13, 927936. Yilmaz, O., 2001, Seismic data analysis. Society of Exploration Geophysicists.  12  Chapter 2  Incoherent noise suppression with curvelet-domain sparsity-promotion 2.1  Background  In seismic data, recorded wavefronts (i.e., reflections) arise from the interaction of the incident wavefield with inhomogeneities in the Earth’s subsurface. The wavefronts can become contaminated with various types of noise during acquisition or due to processing problems (Olhovich, 1964). Thus, separation of signal and noise is an important issue in seismic data processing, particularly in crustal data where the signal-to-noise ratio is low. By noise, we refer to the incoherent (random) noise that is present in the data. Because noise attenuation is a key step in seismic processing, many methods have been developed to suppress such incoherent noise (Donoho, 1995; Zhang and Ulrych, 2003; Jones and Levy, 1987; Ulrych et al., 1999). Predictive deconvolution is one of the most common techniques in seismic prospecting to attenuate noise. It is also known as F-X deconvolution or F-X prediction (Abma and Claerbout, 1995). The prediction filter works in a windowed sense, predicting one dip at a time, and tends to attenuate all other secondary dips in that window. Thus F-X prediction methods are not good for curved events and for events with conflicting dips (e.g. caustics, pinch-outs). More recently, curvelets have been introduced to resolve the issues of directionality of complicated wavefronts as they belong to a family of multiscale, and also multidirectional data expansions (Candes et al., 2005). Each curvelet is associated with a position, frequency and an angle. In seismic data, the underlying geological signal differs from most noises in terms of at 1  A version of this chapter will be submitted for publication. V. Kumar, J. Oueity, F. J. Herrmann and R. M. Clowes (2009) Curvelet denoising: application to crustal reflection data.  13  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion least one of the following features: angle, frequency, or location. For example, a linear dipping event will map to the coefficients of specific angles in the curvelet domain. Consequently, the signal maps to some specific large coefficients whereas noise maps everywhere as small coefficients after curvelet transformation. This makes the curvelet transform an appropriate choice for detecting wave fronts and suppressing noise as noise and signal naturally separate in the curvelet domain (Neelamani et al., 2008; Hennenfent and Herrmannn, 2006). Our objective is to demonstrate the effectiveness of the curvelet transform in removing incoherent noise without harming the signal. We compare three curvelet-based approaches to suppress random noise with application to synthetic data and deep reflection data. The real data were recorded along LITHOPROBE’s SNorCLE Line 1 in the Northwest Territories, Canada (Cook et al., 1999).  2.2  Problem formulation  The forward problem of seismic denoising can be written as: y = m + n,  (2.1)  where y is the known noisy data, m is the unknown noise free data (the model) and n is the zero-centered Gaussian noise. The noise can be white (full band) or colored (band-limited) depending on the previous processing applied on the data. Our objective is to recover m. The curvelet-based methods for noise suppression are discussed next.  2.2.1  Elementwise hard and soft thresholding  The most common method for noise suppression in a transform domain is hard and soft thresholding of transform coefficients (Donoho, 1995). Both ˜ after hard of them are element-wise operation. The estimated model m thresholding is defined as:  ˜ = CH Hλ (Cy), where  m Hλ : (2.2) Hλ (xi ) = xi , if xi ≥ λ   Hλ (xi ) = 0, otherwise  The symbol Hλ represents the hard thresholding function, C is the forward curvelet transform operator and CH is the inverse curvelet transform operator. The symbol H denotes conjugate transpose that is equivalent to the 14  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion ˜ after inverse for this choice of curvelet transform. The estimated model m soft thresholding is defined as: ˜ = CH Sλ (Cy), where m Sλ (xi ) = sgn(xi ) · max(0, |xi | − |λ|),  Sλ :  (2.3)  where Sλ (x) is the soft thresholding operator and λ is the threshold value. Fig.2.1 shows the application of hard and soft thresholding applied on a signal with λ = 2.  (a)  (b)  (c)  Figure 2.1: (a) A seismic signal; dashed lines show the threshold level (λ), (b) Signal after hard thresholding (λ=2), (c) Signal after soft thresholding (λ=2). Both hard and soft thresholding sets entries less than the threshold to zero. In addition, soft thresholding shrinks all remaining entries by the threshold. This moves the signal toward zero.  15  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  2.2.2  One-norm  The denoising problem with curvelet-domain sparsity can be cast into the following constrained optimization problem (Elad et al., 2005): P :  ˜ = arg minx ||x||1 x ˜ = CH x ˜, m  s.t. ||y − Ax||2 ≤ ε,  (2.4)  ˜ represents the estimated curvelet where x is the curvelet transform vector, x ˜ is the estimated model, is proportional to transform coefficient vector, m the noise level. A = CH is the curvelet synthesis operator with a restriction r that allows only those coefficients to enter the solution that survive the best hard threshold. The reason for using the restricted curvelet transform is that we already know from the best hard thresholding (explained later) about the positions of those coefficients that contain no useful signal and we would not like those coefficients to enter our solution while solving Eq. 2.4. Honoring this restriction, we find the sparsest set of curvelet coefficients (by minimizing the one-norm) that explains the data within the noise level. The constrained optimization problem (Eq. 2.4) is solved by a series of the following unconstrained optimization problems (Herrmann and Hennenfent, 2008; Elad et al., 2005; Figueiredo and Nowak, 2003): 1 ˜ = arg min ||y − Ax||22 + λ||x||1 , x x 2  (2.5)  where λ is a regularization parameter that determines the trade-off between data consistency and the sparsity in the curvelet domain. Eq. 2.5 can be solved by iterative soft thresholding (Daubechies et al., 2005). The solution is found by updating x with x ← Sλ (x + AH y − Ax ), where Sλ is the soft thresholding operator defined in Eq. 2.3. We solve a series of such problems (Eq. 2.5) starting with high λ and decreasing the value of λ until ||y − Ax||2 ≤ , which corresponds to the solution of our optimization problem. The algorithm is summarized in Table 2.1.  2.3  Parameter selection for synthetic data  For simplicity of our experiments, the real-valued curvelet transform is used. The curvelet transform has 5 scales, 16 angles at the 2nd coarsest level and curvelets at the finest (fifth) scale. To choose the best threshold level for hard and soft thresholding, we did an experiment to plot a graph of signal-tonoise ratio (SNR) vs. threshold for both approaches. The SNR is calculated 16  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion Choose: L, K Initialize: k = 1, AH y ∞ > λ1 > · · · > λK , x=0 while y − Ax 2 > and k ≤ K do for l = 1 to L x = Sλk x + AH (y − Ax) end for k =k+1 end while m = CH x  Table 2.1: The cooling method with iterative soft thresholding to solve Eq. 3.5. L is the number of inner iterations, K is the number of outer iterations, k is the counter for outer iterations, λ1 · · · λK are the values of λ in decreasing order, other symbols are as explained in the text. by the following formula: SNR = 20 ∗ log10  ||m||2 , ||m − m||2  (2.6)  where m is the true model. As we see from Figs. 2.2, the SNR shows a gradual increase and then starts to decrease. This is expected as the initial threshold level discards the noisy component of the data resulting in an increase of the SNR value. After a certain point (the arrows in Fig. 2.2) the threshold level starts harming the signal components and thus we see a decrease in SNR value. Based on this experiment, the threshold value is chosen for the peak value of SNR as shown by the arrows in Fig. 2.2. Since we know the noise, exact value of noise level ( = ||n||2 ) is used for the one-norm denoising. We used two outer iterations (K=2) and one inner iteration (L=1) for the algorithm (Table 2.1). We choose the value of K and L after intensive experiments in order to get the best possible SNR value. We observe that increasing the number of iterations cause a decrease in SNR value because iterative soft thresholding started introducing bais in the amplitudes of estimated curvelet coefficients. This bias corresponds to harming the coherent components. Since our criterion is to avoid coherent damage, we limit our number of iterations. 17  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion 16  16 Hard Soft  Hard Soft  12  12  10  10 SNR  14  SNR  14  8  8  6  6  4  4  2 0  20  40  60  80  100 120 threshold value  140  160  180  200  2 0  20  40  60  (a)  80  100 120 threshold value  140  160  180  200  (b)  Figure 2.2: (a) Comparison of signal-to-noise ratio (SNR) vs. threshold value for white noise for hard thresholding and soft thresholding of curvelet coefficients, (b) Comparison of signal-to-noise ratio (SNR) vs. threshold value for colored noise for hard thresholding and soft thresholding of curvelet coefficients. The arrows point to the peak SNR values and the corresponding threshold value (λ) is chosen.  2.4  Application to synthetic data  A synthetic shot gather is used for our experiments (same shot gather as described in Chapter 2). The model (Fig. 2.3(a)) lives in the frequency range 5-60 Hz. Noisy data (Fig. 2.3(b)) are obtained by addition of zero-centred white Gaussian noise (σ=50) to the model. We also created a colored noise realization by doing a band-pass filter (5-60 Hz) on the same white noise realization. Fig. 2.3(c) shows the noisy data for additive colored noise. All three curvelet denoising approaches are applied to the noisy data and the signal-to-noise ratio (SNR) is calculated for the denoised data. The estimated models and the differences are shown in Fig. 2.4 and Fig. 2.5. For the one-norm denoising, the difference for white noise (Fig. 2.4(b)) has almost no coherent energy while the difference for colored noise (Fig. 2.5(b)) has some coherent energy. In terms of the difference, both one-norm (Fig. 2.4(b), 2.5(b)) and hard thresholding (Fig. 2.4(d), 2.5(d)) show that they have negligible impact on the signal energy. The difference image of soft thresholding (Fig. 2.4(f), 2.5(f)) shows a significant amount of coherent energy. Because soft thresholding causes biasness in the amplitudes of curvelet coefficients, the noise removed (difference) contains a lot of coherent energy. This is a well-known effect and may not be important for image denoising. However this is important to us as we do not want to harm any coherent energy. 18  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion In terms of SNR (Table 2.2), one-norm denoising has the highest SNR compared to hard and soft thresholding of curvelet coefficients (for both white and colored noise). As we are dealing with redundant transform, onenorm puts a constraint on the coefficients and thus stabilizes the solution. Thus, one-norm gives the highest SNR. Table 2.2: SNR comparison for Noise Data One-norm White 3.44 14.69 Color 7.67 15.44  (a)  pre-stack denoising Hard Soft 14.44 12.77 15.20 14.01  (b)  (c)  Figure 2.3: (a) Noise-free data (model), (b) Model with white noise, (c) Model with colored noise.  19  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  (a)  (b)  (c)  (d)  (e)  (f)  Figure 2.4: Data corrupted with white noise: (a) One-norm approach, (b) Noise removed by one-norm approach, (c) Hard threshold approach, (d) Noise removed by hard threshold approach, (e) Soft threshold approach, (f) Noise removed by soft threshold approach.  20  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  (a)  (b)  (c)  (d)  (e)  (f)  Figure 2.5: Data corrupted with colored noise: (a) One-norm approach, (b) Noise removed by one-norm approach, (c) Hard threshold approach, (d) Noise removed by hard threshold approach, (e) Soft threshold approach, (f) Noise removed by soft threshold approach.  21  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  2.5  Parameter selection for real data  For the case of real data, we kept the same number of scales and angles for curvelets as for synthetic data. Observing the success of one-norm over hard and soft thresholding, we apply one-norm to real data. Choosing the value of noise level ( ) is always a challenge for real data. For this work, we looked at various threshold levels for hard threshold and observed the noise removed for a series of threshold values. Out of these experiments we picked the result that caused reasonable noise attenuation with difference as random noise. We call this best hard threshold result. Based on this best hard thresholding, we defined the restricted curvelet transform CH r and . In other words, we allow only those curvelet coefficients that survive best hard thresholding to enter the solution. As we will see from the examples, choosing the parameters this way caused minimal harm to the signal.  2.6  Application to deep reflection data  Near-vertical incidence crustal seismic data were recorded along SNorCLE Line 1 as part of LITHOPROBE multidisciplinary studies in the PaleoproterozoicArchean domains of the Northwest Territories, Canada (Fig. 2.6). Line 1 is 725 km long and follows the highway from east of Yellowknife to west of Fort Simpson. For the seismic data acquisition, four to five vibroseis trucks were used with a 90 m shot point spacing, 60 m receiver spacing and a total of 404 channels. The sweep frequencies ranged between 10 and 80 Hz with a record length of 32 s (Cook et al., 1999). Fig. 2.7(a) shows that raw shot gather is contaminated with ground roll and random noise, which mask many features, especially deeper reflectors. As observed, the shot gathers are characterized by a gradual decrease in reflectivity with time. The decrease is due to the weak reflections from deeper levels. Some weak reflections are observed at 11.0 s that can be related to the crust-mantle transition. We first balanced the amplitudes of the shot gathers to deal with high amplitude groundroll and bad traces. The balancing is done with an automatic gain control (AGC with 500 ms window, see Fig. 2.7(b)) and subsequently ran the data through the curvelet denoising (Fig. 2.7(e)). Fig. 2.7(c) and Fig. 2.7(d) show the result of applying curvelet denoising to the data of Fig. 2.7(b) and the noise that is removed, respectively. Fig. 2.8(a) shows the stacked section after standard processing (Fig. 2.8(b)) for a 60-km-long segment. The clarity of the recorded seismic data is highly 22  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  39  Cordillera  Slave  Wopmay  1100  CS  1104  Deformed Ancestral North America  1101 Fort Simpson 1107 1112  1109  YK BC 1113  128  Seismic Profile  Extent of Phanerozoic cover 120  Seismic Figure 2.6: Tectonic mapProfile of the acquisition area andoflocation of SNORCLE Extent Phanerozoic cover near vertical incidence (NVI) reflection Line 1. Sediments of the PhaneroShot Points zoic Western Canada Sedimentary Basin overlie the Precambrian domains west of the long dashed line. Short dashed black lines show political boundaries. CS - Coronation Supergroup, SD - Sleepy Dragon, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC - British Columbia, NWT - Northwest Territories, YK - Yukon. Inset (top-left) shows location of map within Canada. Inset (green box) shows the area for the data used.  Fig. 1. Oueity and Clowes Paleoproterozoic slab subduction  23  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion degraded by the low signal-to-noise ratio that precludes detailed interpretation of reflectivity images. The stacked seismic profile shows a decrease in reflectivity between the crust and upper mantle; the crust-mantle transition is identified at 11.0 s (∼34 km). On the southwestern end of the segment shown, the crust-mantle transition becomes more diffuse due to high background noise (arrow at 11 s). Some northeast-dipping events can also be identified on the southwestern end at 4.0 and 6.5 s but can only be followed over very short distances (arrow at 7 s Fig. 2.8(a)). We applied curvelet denoising to the pre-stack and the post-stack data to compare our results with the conventional processing. Processing to stack followed standard procedures used in Lithoprobe studies (Cook et al., 1999) with curvelet denoising before the bandpass filter for the pre-stack data as shown in Fig. 2.7(e). Curvelet denoising is applied on post-stack data as an alternative to F-X prediction as shown in Fig. 2.8(b). Fig. 2.8(c) shows the result from denoising with curvelet processing while Fig. 2.8(e) shows results from F-X prediction. We do quality control of our results by observing the differences (noise removed) between the original data and the denoised data for each of the denoising applications. The noise removed is shown in (Fig. 2.8(d) for curvelet denoising, Fig. 2.8(f) for F-X prediction)  2.7  Results and discussion  In the case of pre-stack data, the results demonstrate that curvelet denoising successfully suppressed both the groundroll energy and random noise, thereby enhancing deep reflections. The Moho is clearly visible in our results at around 11 s. Individual reflectors throughout the crust are enhanced and more clearly visible (Fig. 2.7(c)). The difference plot between the original data and the curvelet result (Fig. 2.7(d)) shows that the noise removed by curvelet processing contains groundroll, corrupted traces and random noise. In the case of post-stack data, the improvement in image resolution is clear when comparing the two data sets before and after the curvelet denoising (Fig. 2.8). The curvelet denoised seismic image (Fig. 2.8(c)) shows a highly reflective crust with significant dipping reflectivity, a relatively flat Moho with an offset of about 1s at about CDP 1100 and a transparent upper mantle. The lower crust contains a series of discrete, east-dipping reflections within a zone that is 2.0 s to 3.0 s thick (about 6.5 km to 10 km) and flattens into the Moho. Below about 9.5 s at the southwestern end of the section, the reflectivity becomes sub-horizontal to about the position of Moho offset. These events are not visible on the original stacked 24  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion image. The Moho is characterized by a sharp, narrow band of reflections (200 - 300 ms) that are piece-wise continuous, except for the offset, over the 60-km-long segment. The difference plot between the original data and the curvelet result (Fig. 2.8(d)) shows that the noise removed by curvelet processing is random in nature. Thus, the technique achieves excellent noise attenuation without removing any coherent events. The F-X predicted image (Fig. 2.8(e)) shows most of the features observed on the curvelet denoised image, but the northeast-dipping and sub-horizontal reflectivity is not clearly expressed. The noise removed by F-X prediction (Fig. 2.8(f)) also contains some coherent reflectivity. Fig. 2.9(a) shows the stack produced by the denoised shot gather. The stacked image did produce a better stack compared to the stack created by conventional processing (Fig. 2.9(b)) but the overall quality of the stack is not as good as what we got after post stack curvelet denoising (Fig. 2.9(c)). The possible reason for this is that stacking is a powerful operation that itself is used to attenuate noise and errors creep in during subsequent processing steps. Nevertheless, the pre-stack denoising produces clearer reflection hyperbolae that are advantageous for pre-stack processing (e.g. velocity analysis, deconvolution, Normal Moveout correction).  25  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  (a)  (b)  (c)  (d) Crooked line geometry Static correction Curvelet Denoising Bandpass filter Deconvolution NMO  Stack  (e)  Figure 2.7: Shot gather from Line 1 segment highlighted in Fig. 3.6. (a) Noisy data without AGC, (b) Noisy data with AGC, (c) Result of applying curvelet denoising (inversion approach) on the shot gather with AGC, (d) Noise removed by curvelet denoising. (e) The flowchart with important processing steps for pre-stack denoising. The curvelet denoising is applied before the bandpass filtering. 26  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion SW  Crooked line geometry  NE  Static correction Bandpass filter Deconvolution NMO Stack  F-X Prediction  (a) SW  (b) NE  SW  (c) SW  NE  (d) NE  (e)  Curvelet Denoising  SW  NE  (f)  Figure 2.8: (a) The original stacked seismic section contaminated with random noise. Arrows indicate some of the features discussed in the text. (b) The flowchart with important processing steps for poststack denoising. The curvelet denoising is applied as an alternative to F-X prediction. (c) Result of applying curvelet denoising (inversion approach) on the same section. (d) Noise removed by curvelet denoising, (e) Result of applying F-X prediction on the same section, (f) Noise removed by F-X prediction. 27  Chapter 2. Incoherent noise suppression with curvelet-domain sparsity-promotion  SW  NE  (a)  (b) SW  NE  (c)  Figure 2.9: (a) Stacked section after pre-stack curvelet denoising, (b) Stacked section after conventional processing (without F-X prediction), (c) Stacked section after post-stack curvelet denoising. The poststack curvelet denoising is observed to be better than the other two.  28  Bibliography Abma, R. and J. Claerbout, 1995, Lateral prediction for noise attenuation by t-x and f-x techniques: Geophysics, 60, 1887–1896. Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899. Cook, F., A. van der Velden, K. Hall, and B. Roberts, 1999, Frozen subduction in canada’s northwest territories: lithoprobe deep lithospheric reflection profiling of the western canadian shield: Tectonics, 18, 124. Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint: Comm. Pure Appl. Math., 57, 1413–1457. Donoho, D., 1995, De-noising by soft thresholding: IEEE Trans. Inform. Theory, 41, 613–627. Elad, M., J. L. Starck, P. Querre, and D. L. Donoho, 2005, Simulataneous cartoon and texture image inpainting using morphological component analysis (MCA): Appl. Comput. harmon. Anal., 19, 340–358. Figueiredo, M. and R. Nowak, 2003, An EM algorithm for wavelet-based image restoration: IEEE Trans. Image Processing, 12, 906–916. Hennenfent, G. and F. J. Herrmannn, 2006, Seismic denoising with nonuniformly sampled curvelets: Computing in Science and Engineering, 8, 16–25. Herrmann, F. J. and G. Hennenfent, 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, 233–248. Jones, I. and S. Levy, 1987, Signal-to-noise enhancement in multichannel seismic data via the karhunen-love transform: Geophysical Prospecting, 35, 12–32. Neelamani, R., A. I. Baumstein, D. G. Gillard, M. T. Hadidi, and W. L. Soroka, 2008, Coherent and random noise attenuation using the curvelet transform: The Leading Edge, 27, 240–248.  29  Bibliography Olhovich, V. A., 1964, The causes of noise in seismic reflection and refraction work: Geophysics, 29, 1015–1030. Ulrych, T., M. Sacchi, and J. Graul, 1999, Signal and noise separation, art and science: Geophysics, 64, 1648– 1656. Zhang, R. and T. J. Ulrych, 2003, Physical wavelet frame denoising: Geophysics, 68, 225–231.  30  Chapter 3  Deconvolution with curvelet-domain sparsity 3.1  Background  Seismic deconvolution is one of the most widely researched seismic signal processing tools. The primary goal of seismic deconvolution is to remove the characteristics of the source wavelet from the recorded seismic time series (data), so that one is ideally left with only the reflection coefficients (reflectivity). The reflection coefficients identify and quantify the impedance mismatches between different geological layers that are of great interest to the geophysicist. The forward problem of seismic deconvolution can be written as: y = s m + n, (3.1) where y is the data, s is the source wavelet, m is the reflectivity model and n is zero-centered white Gaussian noise. The symbol denotes linear convolution. The only quantity we know is the data and we have to estimate the reflectivity and source signature (wavelet). In other words, we have one equation and two unknowns. In practice, we make assumptions that help us to estimate both the reflectivity and source wavelet from data (Jurkevics and Wiggins, 1984). For this work, we assume that we have an estimate of source wavelet and our goal is to estimate the reflectivity given the data and the source wavelet. The source wavelet is also assumed to be stationary, i.e., it does not change with time or offset. In matrix-vector notation the forward problem (Eq. 3.1) can be written as: y = Km + n,  (3.2)  where K is the convolution operator (Toeplitz matrix). The columns of K contain the shifted versions of s. Given K and y, we need to find m. 1  A version of this chapter has been published. V. Kumar and F. J. Herrmann (2008) Deconvolution with curvelet-domain sparsity, In Expanded Abstracts, Society of Exploration Geophysicists, Tulsa.  31  Chapter 3. Deconvolution with curvelet-domain sparsity Since the early 1980’s, researchers have cast this problem as a 1 -norm minimization (Taylor et al., 1979; Oldenburg et al., 1981, 1988; Yarlagadda et al., 1985), where the reflectivity is assumed to comprise spikes. The method is often referred to as sparse-spike deconvolution or spiky deconvolution. The optimization problem solved for spiky deconvolution is given by: ˜ = arg min ||m||1 s.t. ||y − Km||2 ≤ ε, m (3.3) m  ˜ represents the estimate of model vector and is proportional to the where m noise level. In recent work, Herrmann (2005) showed that the assumption of spiky reflectivity is too limited to describe seismic reflectivity. Our objective is to use the curvelet transform in solving the classical deconvolution problem to invert for a non-spiky reflectivity. We show that curvelet sparsity can be used to estimate the non-spiky reflectivity and hence the Earth model. We start with an introduction to our algorithm followed by application to synthetic data for a shot gather and a poststack model.  3.2  Problem formulation  The deconvolution problem with curvelet-domain sparsity can be cast into the following constrained optimization problem: ˜ = arg min ||x||1 x x  s.t. ||y − Ax||2 ≤ ε,  (3.4)  ˜ repwhere A = KCH , CH is the curvelet transform synthesis operator, x resents the estimated transform coefficient vector, and is proportional to the noise level. By solving Eq. 3.4, we find the sparsest set of transform coefficients that explains the data within the noise level (Herrmann and Hen˜ = CH x ˜ . The nenfent, 2008). The final estimated reflectivity is given by m constrained optimization problem (Eq. 3.4) is solved by a series of the following unconstrained optimization problem (Herrmann and Hennenfent, 2008; Hennenfent and Herrmannn, 2006; van den Berg and Friedlander, 2008): 1 ˜ = arg min ||y − Ax||22 + λ||x||1 , x x 2  (3.5)  where λ is a regularization parameter that determines the trade-off between data consistency and the sparsity. To speed up the convergence, we solve these problems (Eq. 3.5) starting with a high λ and decreasing its value until ||y − Ax||2 ≤ , which corresponds to the solution of our optimization problem. The lowering of λ is done in a controlled way so that we reach the 32  Chapter 3. Deconvolution with curvelet-domain sparsity optimum λ rapidly using the spectral projected-gradient algorithm (SPG 1 ) (van den Berg and Friedlander, 2008; Hennenfent et al., 2008). Details about the algorithm are beyond the scope of this thesis but can be found in van den Berg and Friedlander (2008). The algorithm to solve Eq. 3.4 is summarized in Table 3.1. By putting A = K, Eq. 3.4 simplifies to the problem of spiky deconvolution (Eq. 3.3). Choose: K, λ1 > · · · > λK Initialize: k = 1, while y − Ax 2 > and k ≤ K do ˜ = arg minx 12 ||y − Ax||22 + λ||x||1 x k =k+1 end while m = CH x  Table 3.1: Pseudo code to solve the optimization problem shown in Eq. 4.4. K is maximum number of iterations, k is the counter for iterations, other quantities are as stated in the text.  3.3  Parameter selection  The curvelet transform used for this work has 5 scales, 16 angles at the 2nd coarsest level and curvelets at the finest (fifth) scale. In the case of additive white Gaussian noise with standard deviation σ, the square norm of error ||n||22√is a chi-square random variable with mean σ 2 N and standard deviation σ 2 2N , where N is the total number of data points (Cand´es et al., 2005). For this work, we assume that the probability of ||n||22 exceeding its 2 mean by two standard deviations is small. √ The maximum ||n||2 within two 2 standard deviations is given by σ (N + 2 2N ). Thus, we solve Eq. 3.4 with √ 2 = σ 2 (N + 2 2N ).  3.4  Application and results  The reflectivity models used for this work are the same models as described in Chapter 1. Noisy data are obtained by convolving a Ricker wavelet (central frequency=25 Hz) with the reflectivity models followed by the addition of zero-centred Gaussian noise. We apply our curvelet deconvolution to the 33  Chapter 3. Deconvolution with curvelet-domain sparsity noisy data to estimate the reflectivity. For comparison, we also include results obtained by sparse spike deconvolution on the same data. For both approaches, we keep the same (noise-level estimate) to enable a fair comparison. We also calculate signal-to-noise ratio (SNR) for the estimated model given by the expression: SNR = 20 ∗ log10  ||m||2 , ||m − m||2  (3.6)  where m is the true reflectivity model. First, we start with the prestack reflectivity model as shown in Fig. 3.1(a). Fig. 3.1(b) shows the data that are degraded by the convolution operator and the noise. Fig. 3.1(c) and Fig. 3.1(d) show the estimated reflectivity obtained by both algorithms. Fig. 3.2 shows an enhancement for the region within the box of Fig. 3.1. The reflectivity model (Fig. 3.2(a)) is not made up of spikes (does not have zero-order discontinuity). Fig. 3.1(b) shows the reflectivity model after convolution and with added noise. We observe that the coherent reflectivity is degraded. The curvelet estimated reflectivity (Fig. 3.2(c)) is very similar to the true reflectivity (Fig. 3.2(a)). The spiky estimate (Fig. 3.2(d)) consists of sharp spikes that is inconsistent with the true reflectivity (Fig. 3.2(a)). We also apply the algorithm for the poststack model as shown in Fig. 3.3(a). Usually the poststack model has sharp boundaries but to show that our algorithm is suited for fractional order singularities, we did a fractional (half) integration of the model in the frequency domain to obtain a non-spiky model. The non-spiky reflectivity is still sparse in the curvelet domain because of its smooth variation in the offset direction. Fig. 3.3(b) shows the data that were input to the two deconvolution procedures. Fig. 3.3(c) and Fig. 3.3(d) shows the estimated reflectivity obtained by curvelet deconvolution and spiky deconvolution applied to the noisy data shown in Fig. 3.3(b). Fig. 3.4 shows a zoom-in wiggle plot of the reservoir area of Fig. 3.3. As observed, the model (Fig. 3.4(a)) is non-spiky in nature; convolution and additive noise degrades the coherent reflectivity as seen in Fig. 3.4(b). The curvelet-estimated reflectivity (Fig. 3.4(c)) has a close similarity to the true reflectivity (Fig. 3.4(a)). In particular, the curvelet result (Fig. 3.4(c)) improves the resolution of the pinch-out at 1.1 seconds, which was degraded by the convolution operator and noise. The spiky estimate (Fig. 3.4(d)) consists of sharp spikes that do not resemble the true reflectivity (Fig. 3.4(a)). The SNR values are summarized in Table. 3.2. In terms of SNR, we observe that the curvelet-regularized deconvolution has high SNR compared to sparse spike deconvolution for both models. 34  Chapter 3. Deconvolution with curvelet-domain sparsity Table 3.2: SNR comparison for deconvolution methods Data Curvelet decon Spiky decon Poststack -3.32 12.01 8.05 Pre-stack -3.22 14.09 8.27 By exploiting the continuity along reflectors, our curvelet algorithm enhances the degraded frequency components. Sparse spike deconvolution is a trace-by-trace operation, and doesn’t account for the two-dimensional structure of the reflectivity model. On the other hand, curvelets, because of their localization and similarity with curved events, exploit the structure in both time and spatial coordinate.  35  Chapter 3. Deconvolution with curvelet-domain sparsity  (a)  (b)  (c)  (d)  Figure 3.1: (a) True reflectivity model, (b) Noisy data obtained by convolving the reflectivity model with Ricker wavelet and addition of noise (σ=50), (c) Estimated reflectivity with curvelet deconvolution, (d) Estimated reflectivity with sparse spike deconvolution. The box shows the region of interest that is enlarged in Fig. 4.2.  36  Chapter 3. Deconvolution with curvelet-domain sparsity  (a)  (b)  (c)  (d)  Figure 3.2: (a) Zoom-in plot of true reflectivity model, (b) Zoom-in plot of noisy data, (c) Zoom-in plot of estimated reflectivity with curvelet deconvolution, (d) Zoom-in plot of estimated reflectivity with sparse spike deconvolution. Note the resemblance of curvelet-estimated reflectivity with the true reflectivity.  37  Chapter 3. Deconvolution with curvelet-domain sparsity  (a)  (b)  (c)  (d)  Figure 3.3: (a) True reflectivity model, (b) Noisy data obtained by convolving the reflectivity model with Ricker wavelet and addition of noise (σ=.01), (c) Estimated reflectivity with curvelet deconvolution, (d) Estimated reflectivity with sparse spike deconvolution. The box shows the region of interest (reservoir area) for which enlarged versions of the sections are displayed in Fig. 4.4.  38  Chapter 3. Deconvolution with curvelet-domain sparsity  (a)  (b)  (c)  (d)  Figure 3.4: (a) Zoom-in plot of true reflectivity model, (b) Zoom-in plot of noisy data, (c) Zoom-in plot of estimated reflectivity with curvelet deconvolution, (d) Zoom-in plot of estimated reflectivity with sparse spike deconvolution. Note the resemblance of curvelet-estimated reflectivity with the true reflectivity.  39  Bibliography Cand´es, E., J. Romberg, and T. Tao, 2005, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Hennenfent, G. and F. J. Herrmannn, 2006, Seismic denoising with nonuniformly sampled curvelets: Computing in Science and Engineering, 8, 16–25. Hennenfent, G., E. van den Berg, M. P. Friedlander, and F. J. Herrmann, 2008, New insights into one-norm solvers from the pareto curve: Geophysics, 73, A23–A26. Herrmann, F. J., 2005, Seismic deconvolution by atomic decomposition: a parametric approach with sparseness constraints: Integr. Computer-Aided Eng., 12, 69–91. Herrmann, F. J. and G. Hennenfent, 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, 233–248. Jurkevics, A. and R. Wiggins, 1984, A critique of seismic deconvolution methods: Geophysics, 49, 2109–2116. Oldenburg, D. W., S. Levy, and K. P. Whittall, 1981, Wavelet estimation and deconvolution: Geophysics, 46, 1528–1542. Oldenburg, D. W., T. Scheuer, and S. Levy, 1988, Recovery of the acoustic impedance from reflection seismograms: Inversion of geophysical data, 245– 264. Taylor, H. L., S. C. Banks, and J. F. McCoy, 1979, Deconvolution with l1 -norm: Geophysics, 44, 39–52. van den Berg, E. and M. P. Friedlander, 2008, Probing the pareto frontier for basis pursuit solutions: SIAM Journal on Scientific Computing, 31, 890–912. Yarlagadda, R., J. B. Bednar, and T. L. Watt, 1985, Fast algorithms for lp deconvolution: IEEE Trans. on Acoustics, Speech, and Signal Processing, 174–182. 40  Chapter 4  Conclusions 4.1  Curvelet-domain sparsity  The purpose of this thesis is to explore the ability of sparsity-promoting methods as applied to noise suppression and deconvolution. As shown in examples, curvelets provide a sparse representation of seismic data. This sparsity in the curvelet domain provides a scope for new sparsity-promoting methods that can be utilized and applied to seismic data, such as the two algorithms discussed in this thesis. As we see from the presented examples, curvelet-based-sparsity-promotion helps us to overcome the limitations of traditional methods for denoising and deconvolution.  4.2  Application to incoherent noise suppression  We showed how the ability of curvelets to detect wavefronts can be used to suppress incoherent noise. The inversion approach (one-norm denoising) with curvelet-domain sparsity not only recovers the features that are degraded by noise, but also gives the highest SNR for both white and colored noise without compromising the signal. We demonstrate that incoherent noise in deep crustal reflection data can be suppressed using curvelet processing procedures and that the results are better than those from more conventional methods, including F-X prediction. The algorithm is able to remove the incoherent noise and improve the resolution of the seismic image (especially near the Moho) both in pre-stack and post-stack data without compromising the coherent reflectivity.  4.3  Application to deconvolution  We demonstrate that non-spiky reflectivity can be effectively recovered using curvelet-regularized deconvolution and that the results show a better frequency enhancement than those from more conventional methods, including sparse spike deconvolution. The assumption of spiky reflectivity is 41  Chapter 4. Conclusions too limited and may not be applicable in all cases; leading to deteriorated performance of the traditional algorithms. Our algorithm is able to enhance the frequency components that get degraded by noise and a convolution operator. Thus, sparsity of such type of reflectivity in the curvelet-domain is a strong characteristic that is exploited within our deconvolution algorithm.  4.4  Open and future research  The introduction of curvelet processing adds two new tools to the growing family of curvelet-based approaches that includes primary-multiple separation, groundroll removal, seismic data regularization, imaging and inversion. The use of the curvelet domain for seismic processing will most likely continue to expand as new techniques and methods are developed to take advantage of the sparsity. The curvelet domain also allows access to specific scales, angles and locations that have yet to be fully utilized. Another avenue for future work is the influx of new transforms. While the curvelet transform is the most parsimonious of the transforms discussed, many other similar systems are being developed. Shearlets (Labate et al., 2005) are a good example of this. New developments may also lead to a transform that is sparser for seismic data than the curvelet transform. Careful depth and angular weighting within the curvelet domain may further enhance the denoising procedure. The deconvolution algorithm can be incorporated to do simultaneous migration and deconvolution. Both the algorithms can be extended to three dimensions using 3D-curvelets so that continuity of structure in three dimensions can be exploited. One of the most challenging parts of the methods that were described center on parameter selection (especially for the real data). There is currently no fully automated way of selecting optimal parameters for the methods described in this thesis. While fully automated parameter selection may be a bit ambitious, schemes that automatically adjust the parameters for the algorithms described may be possible. A good example of automatically setting the threshold level is the Generalized Cross-Validation (GCV) criterion (Jansen et al., 1997). Another challenge posed by the curvelets is the curvelet vector length that is more than the number of data points. Curvelets are around eight times redundant for two-dimensional data and around twenty-four times redundant in three-dimensions. While two-dimensional data generally are small enough that vector storage is not a problem, problems could be encountered for large three-dimensional data cubes. Thus a parallel version of 42  Chapter 4. Conclusions curvelet algorithm is favorable to handle three-dimensional and even fourdimensional (time-lapse) seismic data.  43  Bibliography Abma, R. and J. Claerbout, 1995, Lateral prediction for noise attenuation by t-x and f-x techniques: Geophysics, 60, 1887–1896. Cand´es, E., J. Romberg, and T. Tao, 2005a, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Candes, E. J., L. Demanet, D. L. Donoho, and L. Ying, 2005b, Fast discrete curvelet transforms: SIAM Multiscale Model. Simul., 5, 861–899. Cand`es, E. J. and D. L. Donoho, 2000, Curvelets – a surprisingly effective nonadaptive representation for objects with edges. Curves and Surfaces. Vanderbilt University Press. Cook, F., A. van der Velden, K. Hall, and B. Roberts, 1999, Frozen subduction in canada’s northwest territories: lithoprobe deep lithospheric reflection profiling of the western canadian shield: Tectonics, 18, 124. Daubechies, I., M. Defrise, and C. de Mol, 2005, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint: Comm. Pure Appl. Math., 57, 1413–1457. Donoho, D., 1995, De-noising by soft thresholding: IEEE Trans. Inform. Theory, 41, 613–627. Elad, M., J. L. Starck, P. Querre, and D. L. Donoho, 2005, Simulataneous cartoon and texture image inpainting using morphological component analysis (MCA): Appl. Comput. harmon. Anal., 19, 340–358. Figueiredo, M. and R. Nowak, 2003, An EM algorithm for wavelet-based image restoration: IEEE Trans. Image Processing, 12, 906–916. Hennenfent, G. and F. J. Herrmann, 2006, Seismic denoising with nonuniformly sampled curvelets: Computing in Science and Engineering, 8, 16–25. Hennenfent, G., E. van den Berg, M. P. Friedlander, and F. J. Herrmann, 2008, New insights into one-norm solvers from the pareto curve: Geophysics, 73, A23–A26.  44  Bibliography Herrmann, F. J., 2005, Seismic deconvolution by atomic decomposition: a parametric approach with sparseness constraints: Integr. Computer-Aided Eng., 12, 69–91. Herrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Nonlinear primary-multiple separation with directional curvelet frames: Geophysical Journal International, 170, 781–799. Herrmann, F. J. and G. Hennenfent, 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, 233–248. Jansen, M., M. Malfait, and A. Bultheel, 1997, Generalized cross-validation for wavelet thresholding: Signal Processing, 56, 33–44. Jones, I. and S. Levy, 1987, Signal-to-noise enhancement in multichannel seismic data via the karhunen-lo(c)ve transform: Geophysical Prospecting, 35, 12–32. Jurkevics, A. and R. Wiggins, 1984, A critique of seismic deconvolution methods: Geophysics, 49, 2109–2116. Labate, D., W.-Q. Lim, J. Kutyniok, and G. Weiss, 2005, Sparse multidimensional represenations using shearlets: Proc. SPIE Wavelets XI, 254– 262. Neelamani, R., A. I. Baumstein, D. G. Gillard, M. T. Hadidi, and W. L. Soroka, 2008, Coherent and random noise attenuation using the curvelet transform: The Leading Edge, 27, 240–248. Oldenburg, D. W., S. Levy, and K. P. Whittall, 1981, Wavelet estimation and deconvolution: Geophysics, 46, 1528–1542. Oldenburg, D. W., T. Scheuer, and S. Levy, 1988, Recovery of the acoustic impedance from reflection seismograms: Inversion of geophysical data, 245– 264. Olhovich, V. A., 1964, The causes of noise in seismic reflection and refraction work: Geophysics, 29, 1015–1030. Taylor, H. L., S. C. Banks, and J. F. McCoy, 1979, Deconvolution with l1 -norm: Geophysics, 44, 39–52. Ulrych, T., M. Sacchi, and J. Graul, 1999, Signal and noise separation, art and science: Geophysics, 64, 1648– 1656. van den Berg, E. and M. P. Friedlander, 2008, Probing the pareto frontier for basis pursuit solutions: SIAM Journal on Scientific Computing, 31, 890–912. 45  Bibliography Versteeg, R., 2008, The marmousi experience: Velocity model determination on a synthetic complex data set: The Leading Edge, 13, 927–936. Yarlagadda, R., J. B. Bednar, and T. L. Watt, 1985, Fast algorithms for lp deconvolution: IEEE Trans. on Acoustics, Speech, and Signal Processing, 174–182. Yilmaz, O., 2001, Seismic data analysis. Society of Exploration Geophysicists. Zhang, R. and T. J. Ulrych, 2003, Physical wavelet frame denoising: Geophysics, 68, 225–231.  46  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052886/manifest

Comment

Related Items