Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Trad, Daniel
Doctor of Philosophy - PhD
University of British Columbia
The Radon transform (RT) has many desirable properties which make it particularly useful for multiple removal and interpolation of seismic reflections. However many practical difficulties arise as a consequence of poor sampling and limited aperture in the offset dimension. Furthermore the ever increasing volumes of seismic data make computing time a key issue in any practical implementation. The standard implementations of the Radon transform used in seismic processing fulfill very well the requirement of a fast transform but do not allow proper handling of problems associated with sampling and aperture. Many of these difficulties can be partially solved by using inverse theory to compute the transform subject to a sparseness constraint. The main thrust and contribution of this thesis lies in the exploration of the sparseness constraint and the design and implementation of the RT. Although this technique reduces sampling and aperture problems, it also considerably increases the computation time. Some of the possible solutions to decrease computation time discussed in this thesis are conjugate gradient methods, irregularly sampled model space and efficient operations with sparse matrices. Real data often have very complicated structure that makes sparseness criteria difficult to implement. Noise, non hyperbolic events and amplitude variation with offset conspire against sparseness. Therefore, the success of this method depends on the applied algorithms and on a delicate balance between sparseness and fit of data. Several real data examples of multiple removal and interpolation show the success of the proposed algorithms to achieve this balance. Another aspect of this thesis is the design of new applications for the Radon transform in seismic processing. Many new applications can be developed by designing new types of Radon operators. Some new applications implemented in this thesis are elliptical and hybrid linear-hyperbolic Radon transforms that should prove to be useful for slant-stack processing and ground roll removal respectively. As a final topic, this thesis presents an application of the RT to the enhancement of a prestack time migration method. Multiple attenuation, artifact filtering and automatic velocity model corrections result from the merging of the R T with equivalent offset migration.
49287596 bytes
application/pdf
IMPLEMENTATIONS AND APPLICATIONS OF T H E SPARSE RADON TRANSFORM By Daniel Trad Geofisico, Universidad Nacional de San Juan, Argentina, 1994 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S Department of E A R T H A N D O C E A N S C I E N C E S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A September 2001 © Daniel Trad, 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences The University of British Columbia 6339 Stores Road Vancouver, Canada V 6 T 1Z4 Date: October a * 2.001 A b s t r a c t The Radon transform (RT) has many desirable properties which make it particularly use-ful for multiple removal and interpolation of seismic reflections. However many practical difficulties arise as a consequence of poor sampling and limited aperture in the offset dimension. Furthermore the ever increasing volumes of seismic data make computing time a key issue in any practical implementation. The standard implementations of the Radon transform used in seismic processing fulfill very well the requirement of a fast transform but do not allow proper handling of problems associated with sampling and aperture. Many of these difficulties can be partially solved by using inverse theory to compute the transform subject to a sparseness constraint. The main thrust and contribution of this thesis lies in the exploration of the sparseness constraint and the design and implementation of the R T . Although this technique reduces sampling and aperture problems, it also considerably increases the computation time. Some of the possible solutions to decrease computation time discussed in this thesis are conjugate gradient methods, irregularly sampled model space and efficient operations with sparse matrices. Real data often have very complicated structure that makes sparseness criteria dif-ficult to implement. Noise, non hyperbolic events and amplitude variation with offset conspire against sparseness. Therefore, the success of this method depends on the ap-plied algorithms and on a delicate balance between sparseness and fit of data. Several real data examples of multiple removal and interpolation show the success of the proposed algorithms to achieve this balance. ii Another aspect of this thesis is the design of new applications for the Radon transform in seismic processing. Many new applications can be developed by designing new types of Radon operators. Some new applications implemented in this thesis are elliptical and hybrid linear-hyperbolic Radon transforms that should prove to be useful for slant-stack processing and ground roll removal respectively. As a final topic, this thesis presents an application of the R T to the enhancement of a prestack time migration method. Multiple attenuation, artifact filtering and auto-matic velocity model corrections result from the merging of the R T with equivalent offset migration. iii T a b l e o f C o n t e n t s A b s t r a c t i i L i s t o f F i g u r e s v i i i A c k n o w l e d g e m e n t x v i i i 1 I n t r o d u c t i o n 1 1.1 The Radon transform in different fields 2 1.2 Invertibility of The Radon transform in continuous domains 5 1.3 The RT in the discrete world 8 1.4 Different Radon basis functions 9 1.4.1 Linear RT 10 1.4.2 Parabolic and Pseudo-hyperbolic RT 11 1.4.3 Hyperbolic RT 12 1.5 From processing to inversion: the sparse RT 13 1.6 Time versus frequency domain 18 1.7 Scope and contribution of this thesis 19 2 I n v e r s i o n for sparse m o d e l s : a l g o r i t h m s a n d e x a m p l e s 21 2.1 Some basics about inverse theory 21 2.1.1 The Bayesian approach 23 2.2 Hyperparameters 25 2.3 Methods for sparse inversion 27 iv 2.3.1 The Iteratively Reweighted Least Squares method (IRLS) 27 2.3.2 The conjugate gradient algorithm 31 2.4 Example 1: Sparse Deconvolution 34 2.4.1 Deconvolution in terms of inversion 35 2.4.2 S V D and G S V D 37 2.4.3 Preconditioned C G 47 2.5 Example 2: Denoising by sparse hybrid inversion 48 2.5.1 Problem 50 2.5.2 The iterative re-weighted least square method for the hybrid inversion 51 2.5.3 Log Barrier interior point method 52 2.5.4 Results 54 2.6 Example 3: Re-sampling with the Fourier transform (FT) 64 2.7 Conclusions 65 3 T i m e invariant R T 69 3.1 Synthetic data examples 70 3.2 Non sparse R T 71 3.3 Sparse R T 75 3.3.1 Weighted conjugate gradients least squares ( W T C G L S ) 76 3.3.2 C G - F F T R T 80 3.3.3 De-aliased R T 85 3.3.4 Hybrid time-frequency domain R T 86 3.3.5 Other methods 91 3.4 Real data examples for marine and land data 94 3.5 Conclusions 108 v 4 T i m e variant R a d o n transforms 110 4.1 P R T v s H R T 110 4.2 Time variant Radon transforms 112 4.2.1 The Hyperbolic Radon Transform 112 4.2.2 Elliptical Radon Transform 113 4.3 Implementation 114 4.3.1 C G solution for the linear problem 114 4.3.2 Fast sparse implementation for the operators 115 4.3.3 Velocity grid 118 4.4 Examples 119 4.5 Discussion 124 4.6 Conclusions 126 5 A hybr id L inear -Hyperbol ic Radon transform 135 5.1 The hybrid operator 135 5.2 In quest of a zero nullspace 138 5.3 Inversion for sparse hybrid models 140 5.4 Examples 141 5.5 Discussion and conclusions 143 6 Al i a s ing 152 6.1 Aliasing in different flavors 152 6.2 Different views of aliasing 154 6.3 Operator aliasing and inversion 158 6.4 Irregular sampling 160 6.5 Alias attenuation by the time domain sparse R T 164 6.6 Examples 165 vi 6.7 Conclusions 168 7 Radon transform applied for migration 174 7.1 E O M 175 7.2 The E O M - R T method 178 7.2.1 Denoising 179 7.2.2 Focusing . . 181 7.3 Implementation 181 7.3.1 Equivalent Offset Migration and CSPs 182 7.3.2 Radon Transform 183 7.4 Examples 186 7.5 Discussion 190 7.6 Conclusions 191 8 Discussion and Conclusion 200 References 203 Appendices 210 A Table of Symbols 210 vii List of Figures 1.1 The Radon transform is a line integral along a particular kind of curve. . 3 1.2 Polar R T used in medical tomography, (a) Data (b) Radon transform . . 4 1.3 Radon transform operator in the frequency domain 8 1.4 The Fourier projection theorem (from Stewart, 1991) 11 1.5 Differences between smooth and sparse R T 18 1.6 Time invariant vs time variant R T 19 2.1 (a) Ricker wavelet, (b) reflectivity (c) convolution of the wavelet with reflectivity (clean data), (d) data plus noise 35 2.2 (a) True reflectivity (b) Truncated S V D solution for noise-free data, (c) T S V D for noisy data 38 2.3 (a) Exact data (b) Noisy data 40 2.4 Four of the first singular vectors for the model space of the given kernel . 41 2.5 Generalized singular values for the convolution kernel 42 2.6 Four generalized singular vectors for the convolution kernel and the sparse constraint 43 2.7 (a) The solution constructed by using the first 27 singular vectors (no prior information), (b) The solution obtained by adding the last 9 generalized singular vectors (second iteration of IRLS). (c) The same but after the third iteration of IRLS 44 2.8 First four singular vectors for the modified convolutional kernel AL-1. . 45 viii 2.9 Evolution of the solution built by using the modified kernel A L - 1 as the iteration proceeds and L is more precisely computed 46 2.10 The solution built by using Preconditioned C G . (a) Initial model, (b) Solution by solving the expanded system of equations 47 2.11 The solution built by using left preconditioned C G . (a) First iteration, L = I because there is not prior information yet. (b) fifth iteration, (c) Tenth iteration 48 2.12 Basis functions and models, (a) Signal basis, (b) True model, (c) Noise basis, (d) Model for noise 55 2.13 (a) Signal, (b) Noise, (c) Data 56 2.14 (a) Zero order regularization (b) IRLS with correct model weights (perfect prior), (c) IRLS with perfect model weights for the noise, (d) IRLS with model weights obtained from the data 57 2.15 Results for IRLS with W T C G L S (blocky noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and pre-dicted noise (line) 58 2.16 Results for Log Barrier method (blocky noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and pre-dicted noise (line) 59 2.17 Results for nonlinear C G method (blocky noise test) with preconditioning. Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line) 60 ix 2.18 Basis functions and models, (a) Signal basis function, (b) True model, (c) One of the basis functions for the coherent noise, (d) Model for noise . . 61 2.19 (a) Signal, (b) Noise, (c) Data 62 2.20 Results for IRLS with W T C G L S (harmonic noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line) 63 2.21 Results for Log Barrier method (harmonic noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and pre-dicted noise (line) 63 2.22 Fourier transform of irregularly sampled data, a) F T of regularly sampled data, b) D F T of irregularly sampled data, c) Least squares F T of the same data, d) Sparse F T of the same data 67 2.23 Fourier transform and irregular data, a) true signal with regular sampling (line) and predicted by F T (dots), b) irregularly sampled data (line) and predicted data by D F T (dots), c) data (line) and predicted signal by L S F T (dots), d) data (line) and predicted signal by sparse F T (dots) 68 3.1 (a) Synthetic gather, (b) Non-Sparse R T . (c) Spectrum of (b). (d) Recov-ered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a) 72 3.2 (a) Synthetic data, (b) Synthetic gather with gaps, (c) Non-Sparse R T . (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a) 73 x 3.3 (a) Synthetic gather, (b) Sparse R T obtained with the W T C G L S method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a) 77 3.4 (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse R T obtained with the W T C G L S method, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a) 78 3.5 (a) Synthetic gather, (b) Sparse R T obtained with the C G - F F T method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a) 81 3.6 (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse R T obtained with the C G - F F T method, (d) Spectrum of (c). (e) interpolated data from (b) . (f) Residuals obtained by subtracting (e) from (a) 82 3.7 (a) Synthetic gather, (b) Sparse R T obtained with the de-aliased R T method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a) 87 3.8 (a) Synthetic gather, (b) Sparse R T obtained with the time-frequency domain C G method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a) 88 3.9 (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse R T obtained with the time-frequency domain C G method, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a) 89 3.10 (a) Synthetic data, (b) Synthetic gather with gaps, (c) Irregular R T obtained with a frequency dependent q sampling interval, (d) Spectrum of (c) . (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a) 92 xi 3.11 (a) Marine C D P . (b) Primaries obtained by subtraction, (c) Multiples obtained from the inverse R T of (e). (d) Non-Sparse R T . (e) R T after mute of primaries. Bottom - Semblance analysis for: (f) data, (g) primaries, and (h) multiples 95 3.12 (a) Marine C D P . (b) Primaries obtained by subtraction, (c) Multiples obtained from inverse R T of (e). (d) Sparse R T - W T C G L S . (e) R T after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples 96 3.13 (a) Marine C D P . (b) Primaries obtained by subtraction, (c) Multiples obtained from inverse R T of (e). (d) Sparse R T - T i m e frequency domain, (e) R T after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples 97 3.14 (a) The Mississipi Canyon data after N M O correction and stack 100 3.15 (a) The Mississipi Canyon data after multiple removal with P R T ( W T C G L S ) , N M O correction and stack 101 3.16 (a) Marine C M P . (b) Non sparse R T obtained with Levinson recursion. (c) Interpolated data obtained from (b). (d) Residuals 102 3.17 (a) Marine C M P . (b) Sparse R T obtained with the W T C G L S method, (c) Interpolated data obtained from (b). (d) Residuals 103 3.18 (a) C M P from land data, (b) Non sparse R T . (c) Sparse R T obtained with the W T C G L S method, (d) Sparse R T obtained with the time-frequency domain method 105 3.19 (a) C M P from land data, (b-c-d) Interpolation by using b) non sparse R T . c) Sparse R T obtained by the W T C G L S method, d) Sparse R T obtained by the time-frequency domain method 106 xn 3.20 (a) C M P from land data after eliminating some traces, (b-c-d) Interpola-tion by using b) non sparse R T . c) Sparse R T obtained by the W T C G L S method, d) Sparse R T obtained by the time-frequency domain method. . 107 4.1 Algorithm P C G N R 116 4.2 Sparse matrix operations. A two column matrix connects the data and model spaces 117 4.3 Sparse operator 118 4.4 (a) Synthetic C M P gather with a near offset gap and very limited aperture. (b) Sparse hyperbolic R T . Horizontal axis represents slowness at zero time. (c) Gather produced by the Radon model, (d) Residuals between complete synthetic gather (not shown) and recovered after R T 127 4.5 (a) The same synthetic C M P gather with a near offset gap and very limited aperture, (b) Nonsparse parabolic R T computed in the frequency domain (after N M O ) . Horizontal axis represents residual moveout per offset square, (c) Gather produced by the Radon model, (d) Residuals between complete synthetic gather (not shown) and recovered after R T 128 4.6 (a) Complete simulated gather, (b) Gather a) with missing traces and truncation, (c) f-k spectrum from a), (d) Slant stack from the complete gather, (e) Slant stack from incomplete gather b). (f) Residuals (d-e) (g) Velocity model obtained from e) by inverse elliptical R T . (h) Forward E R T from g). (i) Residuals (d-h) 129 4.7 (a) C D P from marine data with three artificial gaps, (b) Hyperbolic Radon Transform with irregular Radon space, (c) Original and interpolated traces using the model in b). (d) Difference between original (complete set) and synthetic traces (Gather generated from b) 130 xiii 4.8 (a-b-c) Near offsets: Gap within -460 to 0 m (a) Original data, (b) Inter-polated with sparse P R T . (c) The same with H R T . (d-e-f) Middle offsets: Gap between -2600 to -2140 m (d) Original data, (e) Interpolated with sparse P R T . (f) The same with H R T . (g-h-i) Large offsets: Gap between -4740 to -4276 m (g) Original data, (h) Interpolated with sparse P R T . (i) The same with H R T 131 4.9 (a) C M P gather at location 1000. (b) Linear R T of (a), (c) C M P gather at location 1000 with gaps and one every two traces, d) Linear R T of (c). (e) Model obtained from (d) using elliptical R T . (f) Linear R T of (c) obtained from (e) 132 4.10 (a) Marine C D P . (b) Primaries obtained by substraction. (c) Multiples obtained from inverse R T of (e). (d) Sparse H R T . (e) R T after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples 133 4.11 (a) The Mississipi Canyon data after multiple removal with H R T , N M O correction and stack 134 5.1 Hybrid inversion with sparseness constraints 139 5.2 (a) Shot gather, (b) Recovered data after Parabolic R T . (c) Recovered data after L R T . (d) Recovered gather after L R T + P R T 145 5.3 (a) Shot gather, (b) Residuals after Parabolic R T . (c) Residuals after L R T . (d) Residuals after L R T + P R T 146 5.4 (a) Synthetic gather with hyperbolic and linear events, (b) Data in the f — kx domain, (c) / — k Filtered data, (d) Filtered / — kx domain. . . 147 xiv 5.5 (a) Synthetic gather with hyperbolic and linear events, (b) Sparse linear-hyperbolic R T . Horizontal axis shows trace numbers. Slopes for the left hand side (linear space) and curvatures for the right hand side increase to the right side, (c) Predicted noise obtained by muting the parabolic space. (d) Signal obtained by subtracting c) from a) 148 5.6 (a) Shot gather, (b) Combined Linear + Pseudo Hyperbolic Radon do-main, for negative offsets (traces 1-160), and for positive offsets (traces 161-320). (c) Recovered gather from linear space only, (d) Signal obtained by difference between (a) and (b) 149 5.7 (a) Synthetic gather with hyperbolic and linear events, (b) Non-Sparse Linear-hyperbolic R T . (c) Spectrum of (b). (d) Recovered data from hy-perbolic model, (e) Recovered data from linear model 150 5.8 (a) Synthetic gather with hyperbolic and linear events, (b) Sparse Linear-hyperbolic R T . (c) Spectrum of (b). (d) Recovered data from hyperbolic model, (e) Recovered data from linear model 151 6.1 (a) Regular sampling function, (b) Fourier transform of sampling function, (c) Irregular sampling function, (d) Fourier Transform of Sampling Function . . . 154 6.2 Condition number as a function of frequency, a) L R T . Note the jump at / = 40Hz where aliasing arises - b) P R T 160 6.3 Smearing filter LHL (linear r — p ) a) / = 10Hz, regular offset no aliasing, b) / = 25Hz, regular offset, aliasing, (c) / = 25Hz, irregular offset, no aliasing 161 6.4 Resolution and smearing filters (linear R T ) a) resolution filter from smear-ing filter b) (regular offset) - c) resolution filter from smearing filter d) (irregular offset) 162 xv 6.5 Smearing filter LHL (parabolic R T ) a) / = 10Hz regular offset - b) / = 25Hz regular offset -c) / = 25Hz irregular offset 163 6.6 (a) Regularly sampled time series with fx = 0.35Hz, f2 = 0.7Hz, (b) Irregularly sampled time series, (c) Fourier transform of time series in (a), (d) Fourier transform of time series in (b) 165 6.7 Two events sampled with (a) regular and (b) irregular sampling, (c) Radon domain (f-p) from (a), showing the first aliased event (left hand side) for the first high slope event and the first aliased event (right hand side) for the second lower slope event, (d) Radon domain (f-p) from (b), with attenuated alias events 170 6.8 Filtering low velocity events in the Radon domain beyond alias from regu-lar sampling data does not remove completely the event but it does almost completely if the data are irregularly sampled (a) Radon domain from regular sampled data (b) Radon domain from irregular sampled data (c) Filtered from (a) and (d) Filtered from (b) 171 6.9 (a) Synthetic aliased data, (b) f-k spectrum of the aliased data (c) Non sparse L R T . (d) Spectrum of (c). (e) Recovered data from (c) after mut-ing slopes smaller than 0.002s/m. (f) Sparse L R T obtained with a time domain method, (g) Spectrum of (f). (h) Recovered data from (f) after the same muting 172 6.10 (a) Synthetic aliased data (parabolas), (b) f-k spectrum of the aliased data (c) Non sparse P R T . (d) Spectrum of (c). (e) Recovered data from (c) after muting events with Radon parameter smaller than 3 .10 - 7 s/m2. (f) Sparse P R T obtained with a time domain method, (g) Spectrum of (f). (h) Recovered data from (f) after the same muting 173 xvi 7.1 The prestack volume of seismic data and CSPs mapping 175 7.2 Relation between two way traveltime equation and the equivalent offset traveltime equation 177 7.3 Difference between mapping to one trace (Kirchoff) and several traces (RT).180 7.4 The E O M + Radon transform 182 7.5 Computing velocity corrections on the output Radon space 185 7.6 (a) Model, (b) E O M using correct R M S velocities, (c) E O M using wrong velocities, (d) Radon E O M using wrong velocities 193 7.7 Signal to noise ratio =1 (a) E O M (b) Radon E O M 194 7.8 Semblance analysis on clean gathers located at 300m, 500m and 9000m. first row: C M P gathers, second row C S P gathers, third row, C S P gathers + R T (c) 195 7.9 Semblance analysis on noisy gathers located at 300m, 500m and 9000m. first row: C M P gathers, second row C S P gathers, third row, C S P gathers + R T (c) 196 7.10 (a) Stack, (b) E O M , (c) R A D O N - E O M , (d) Radon (in frequency) -E O M 197 7.11 (a) C M P gather at location 1000. (b) C S P gather at location 1000. (c) C S P gather after filtering horizontal artifacts in the R T domain and sub-tracting multiples, (d-e-f) Semblance analysis from (a-b-c) respectively. . 198 7.12 (a) C S P gather at location 1000 obtained without the C M P gather 1000. (b) The same gather after multiple subtraction, (c) The same gather after multiple subtraction and filtering artifacts, (d-e-f) Semblance analysis from (a-b-c) respectively 199 x v i i Acknowledgement I would consider myself a very poor human being if, after four years of hard work, I had no one to thank. Fortunately, this is not the case. This thesis would not have been possible without the large or small contributions of many people. First of all, I want to thank Tadeusz Ulrych for having me on his team, but more for being my friend. He gave me encouragement and guidance. He trusted me when I did not trust myself. Little by little he covered the holes in my knowledge about signal processing and inversion. Doug Oldenburg, Michael Bostock, Bruce Buffet and other excellent professors taught me enough to start me on my journey. This thesis contains many of the concepts that I learned from them. Mauricio Sacchi, and Eldad Haber showed me the path among the clouds, gave me ideas and helped until the sky was clear. The universe is too complicated to look at it with a single pair of eyes. I am grateful to my friends, James Irving, Sam Kaplan, Kris Innanen, Albert Zhang, Shougan Song and many more who often opened my mind through discussions about the wonderful world of science, and most important, they gave me friendship in difficult times. To Jandyr Travassos, who taught me the "abc" of science, and showed me that there are no limits, only obstacles, and that obstacles can be overcome with courage and persistence. On a more fundamental level, I give my sincere acknowledgment to my parents, who gave me the support, comprehension, love, and who taught me everything they know, enabling me to face life's challenges. They gave me wings to fly, but waited patiently until they grew, never asking anything from me, other than to be a good human being. I thank in a special way Betina, who was my wife, my best friend, the column that sustained my smile, and the wind that pushed my boat in the sea of life. xviii To all of them, and many more that made me feel at home in this country blessed by God, thank you... xix C h a p t e r 1 I n t r o d u c t i o n / keep six honest serving-men (They taught me all I know); Their names are What and Why and When And How and Where and Who. I send them over land and sea, I send them east and west; Rudyard Kipling This thesis addresses different implementations and applications of the Radon trans-form (RT). In particular, it deals with the sparse R T . This transform maps events with different shapes, curvatures or slopes, from images, seismic gathers or other similar spaces, to a new space, called the Radon domain or transform domain. Most R T applications in seismic processing consist of mapping mixed and overlapping events in the seismic gather to a new transform domain where they can be separated. A transformation with higher resolution, called the sparse R T , can be obtained by applying inverse techniques with model dependent regularization (Thorson and Claerbout, 1985; Sacchi and Ulrych, 1995). The R T is currently calculated in the time-offset, frequency-offset or frequency-wavenumber domains, with linear, parabolic or hyperbolic basis functions, each one with many possi-ble implementations. Problems appear as a result of discrete sampling in offset and time, discrete velocity parameterization, finite aperture in offset and velocities and missing data due in particular to low fold. Even though the sparse R T has been recognized as being superior to the standard R T , it is not extensively used in routine seismic processing. There are many reasons for 1 Chapter 1. Introduction 2 this, including large computation time, the introduction of artifacts and difficulties in setting inversion parameters. Devising robust solutions to these problems is part of the objective of this work. This thesis is strongly directed towards the solution of practical aspects related to the R T since it is in this area that the sparse R T has shown deficiencies. A large part of this thesis therefore deals with computer implementations and data examples. 1.1 The R a d o n transform i n different fields Scientists have dedicated much effort and time to the design of transformations that map data into new representations or domains. The key idea is that some features of the data are most easily identified in certain domains. The Fourier transform, for example, is ubiq-uitous in many disciplines. Its basis functions are eigenvectors of differential operators that describe many physical processes. The Wavelet transform is another mathematical tool that is receiving much attention due mainly to its localization properties. The Radon transform inherits many of the problems of other transforms. For example, its resolution decreases when the data consist of a discrete representation with limited aperture. The goal of increasing resolution is very important for the R T because its main use is to separate events that map very closely in the transformed domain. Developed by J . Radon in 1917, the Radon transform has become very popular in seis-mic data processing, image processing and geophysical and medical tomography. Radon (1917) showed that if a complete set of projections of some object's parameters were measured then the parameters of the object can be calculated from these projections. In seismic processing the R T is presented as a set of line integrals in the seismic gath-ers. The result of these integrals, i.e. the projections, are represented in a new domain called the Radon panel or transformed domain. The line integrals can follow different Chapter 1. Introduction 3 data space model space H„. m(q,r)= \d[h,t = ft^{z,q,h)¥h Figure 1.1: The Radon transform is a line integral along a particular kind of curve. curves, but usually only straight lines, parabolas or hyperbolas are used. Another way of looking at the RT is a decomposition of the data with a particular set of basis functions. Let d(h, t) represent the signal in time t and offset h, the RT can be calculated as /+0O d(h,t = fu.T(T,q,h))dh -oo (1.1) where t — fti_T(r, q, h) defines the integration path as a function of the model and data space coordinates (r,q,h). q is a parameter that defines the shape of the curve (see Figure 1.1). In tomography, the RT acquires a different role. The tomographic problem is usually posed as k = / p{x,y,z)dl'(x,y,z) (1.2) where U is a particular ray between source and receiver, p(x, y, z) is the slowness at (.x, y, z) and U is the traveltime. After discretization dl becomes Al, the length of the path for the ray i in a cell at (x, y, z). In this expression the goal is to find the slowness Chapter 1. Introduction 4 (a) (b) x 1 0 5 20 40 60 80 100 120 9 (degrees) Figure 1.2: Polar RT used in medical tomography, (a) Data (b) Radon transform for different spatial coordinates (x, y, z). In medical tomography, the linear RT is applied using polar coordinates, (Rf) (p, 9) = J J dxdyf(x, y)5{xcos9 + ysinO - p) (1.3) where / is a continuous function with compact support, 8 is the Dirac delta function, p is the perpendicular distance from one particular projection to the origin, and 9 is the angle formed by the distance vector with the vertical (an example of polar RT is shown in Figure 1.2). In Equation (1.3), the image (Rf)(p,9) is found experimentally by exposing the patient's body to X rays. The function f(x,y) in equation (1.3) is the physical image of the patient, and represents the unknown in the problem. It is computed with one of many different approaches for inverting the transform, that depends on the dimensionality of the problem. For example, for 2 dimensional problems, the image is computed by using a time derivative and the Hilbert transform, f(x,y) = f I Jr(Rf(u}19))\u]\exp(iu}(xcos9 + ysin9))du;d9. (1.4) J0,7T J — 0 0 , + 0 0 Chapter 1. Introduction 5 A general problem in all applications is insufficient data as the R T requires infinite data to obtain a fully invertible transform. The particular manner in which this prob-lem manifests itself depends on the application and different approaches are taken to solve it. In medical imaging, Rf(0) is measured over a full aperture of 360 degrees us-ing a controlled experiment with a high sampling density. In seismic data processing the issue of insufficient data arises because of limited aperture recordings and missing traces. Inversion techniques have been applied to overcome this problem. In medical tomography, excessive x-ray exposure is dangerous, therefore new approaches using the wavelet transform are being considered to obtain more localized information with less global data (Olson and DeStefano, 1994). A localized R T has been attempted as well by using a window function W\, (Rlocf)(p,6,q) = j jdxdyf(x,y)5(xcos8 + ysin6 — p)W\(—xsin9 + ycos9 — q). (1.5) 1.2 Invertibility of The Radon transform in continuous domains Let us define two Hilbert spaces, one containing all possible data vectors, and another containing all possible model vectors with inner product defined by < . , . >. Given a transformation or mapping L between the two spaces, the data space d and the model space m, the adjoint transformation L* is defined as the one satisfying {L,L*) is called an adjoint pair. For the operator defined in the equation (1.1), the adjoint transform is < d, Lm >=< L*d, m > . (1.6) (1.7) Chapter 1. Introduction 6 Because the Radon transform is not an orthogonal operator, the result of the adjoint transformation, d(h, £), is not equal to the original data set d(h, t) computed using equa-tion (1.1). If the curve defined by t = / t < _ T is time invariant, for example t = T + q<b(h), the adjoint transform for the Radon transform becomes d(h,t)= m(q,T = t-q<f>(h))dq. (1.8) 7 — 0 0 When the curves are time invariant, as in equation (1.8), the transform can be per-formed in the frequency domain using the Fourier transform shift theorem, i.e. the time shift factor becomes a multiplication with a complex exponential. The forward-adjoint transformations become /+oo d(h, OJ) exp(icjq<p(h))dh, (1.9) -oo ~ r+oo d(h, OJ) = / m(q, OJ) exp(—iojqcf)(h))dq, (1-10) J—oo Following Zhou and Greenhalgh (1994), the relation between d(h,u) and d(h,oj) can be obtained by replacing m(q,w) in the expression for d(h,u), equation (1.9), _ r+oo r+oo d(h,u)= / d(h',oj)exp(tojq((/){h) - (f){ti)))dqdti = d{h,oj) * p(h,oo). (1.11) J—oo J — oo where * denotes convolution and p(h,oo) is called the p filter. To recover the original data d(h,oj) from d(h,u), the filter p(h,u) has to be deconvolved for every frequency OJ. Note that this is a spatial and not a time deconvolution. The shape of p(h, OJ) depends on the characteristics of the data (aperture, offset sampling, etc.) and on the shape of the integration curve. For linear r-p and data with infinite aperture and dense sampling the p filter can be deconvolved from d, giving an exact inverse R T . For a finite aperture and non-dense Chapter 1. Introduction 7 sampling the p filter is band limited and must be deconvolved using some inversion algorithm utilizing regularization. It is advantageous to perform the inversion or deconvolution in the transformed Radon domain because it is in this domain that high resolution is required. With this purpose in mind, the definition of the forward RT operator and its adjoint is changed to construct another RT pair (Zhou and Greenhalgh, 1994). First, the model m(q, r) is defined as the quantity that generates the data according to the Radon integral, /+oo m(q,T = fT^t(t,q,h))dq. (1.12) -oo and the adjoint model m(q, r) is defined as the back projection of the data into the Radon space, /+oo d(h,t = ft+_T(r,q,h))dh. (1.13) -oo The RT is computed by applying some procedure to m(q, r) to obtain the model m(q, r). For time invariant curves, for example t = T + q(/)(h), these transformations can be expressed in the frequency domain, /+oo m(q,Lo) ex.p(—iuiq(/)(h))dq, (1-14) -oo /+oo d(h,w) exp(iuiq(j)(h))dh. (1-15) -oo As before, the connection between m(q, LO) and the true model m(q, co) can be obtained by replacing d(h,co) in the last equation (1.15). /+oo r+oo / m(q', w) exp(iuj(j)(h)(q — q'))dhdq' (1-16) -oo J—oo or, m(q,u) = m(q,u)*g(q,u}). (1.17) Chapter 1. Introduction 8 d=Lm "•(<»<i.„/j = exp(-iafo\q}) nitvuij) '""•"l,n„,> Figure 1.3: Radon transform operator in the frequency domain. where g(q, w) is the spreading function. This equation means that for every frequency w, to find the model that generates the data d(h, co) requires a deconvolution of the spreading function g(q,uj) from the quantity computed by the integral (1.15) (deconvolution on the Radon coordinate q). 1.3 The R T i n the discrete wor ld After discretization, the transformations considered can be performed by means of matrix multiplications (Beylkin, 1987) (Figure 1.3), rH = FHWdd = LHd, (1.18) d = F W m m = L m , (1.19) where, if q defines the Radon parameter and h the offset, = diag[Ahmin, • • •, Ahmax] and W m = diag[Aqmin, • • •, Aqmax] account for the irregular geometry ofthe array and for irregular sampling of the variable q, respectively. The matrix F has dimension (nh x nq) and elements given by Chapter 1. Introduction 9 Fij = exp(-zu;0(fc,)?j)- (1-20) The operators L and LH form an adjoint but not an inverse pair (the transformation is non-orthogonal). The non-orthogonality implies that the recovered data have been altered by the projection onto the non-orthogonal basis functions of the model space. This alteration is partially eliminated using an orthogonality filter. The notation m reminds us that something else has to be done in the adjoint model m to get m such that d = Lm. For example, for the least-squares method, the R T is computed as m = ( L w L ) - 1 L / / d (1.21) and the data are recovered with d = Lm. (1.22) When the R T is defined in terms of modeling and inversion, as is done in this thesis, the mapping from model to data space is called the forward problem, and the reverse mapping from data to model space is called the adjoint transformation. The adjoint is a first approximation to the inverse operator (Claerbout, 1992). 1.4 Different Radon basis functions A n important distinction is one between time invariant RTs, where the integration is along curves with the general expression, t = r + q<j)(h), and time variant RTs. The former is usually computed in the frequency domain, the latter in the time domain. Three different line integrals are commonly used in the frequency domain, giving three RTs: linear R T (LRT) or slant stack, parabolic R T (PRT) and pseudo hyperbolic ( P H R T ) . Time variant curves are not often used because they require computation in the time Chapter 1. Introduction 10 domain. Two varieties of time variant RTs useful for seismic data are hyperbolic and elliptical, but other curves may be interesting for special cases. Implementations in the time domain are the subject of continuing research, even for time invariant curves (Cary, 1998), because they allow for sparseness both in the time and offset domains. Direct solutions in the time domain are not feasible in general, but because the Radon operator in this domain is sparse, approximate solutions can be obtained using steepest descent (Thorson and Claerbout, 1985), conjugate gradients algorithms (Cary, 1998; Trad et al., 2001b) or Lanczos bidiagonalization (Sacchi, 1996). 1.4.1 Linear R T The linear Radon transform (LRT) , also called slant stack, can be approximately inter-preted as a plane wave decomposition of seismic data (Treitel et al., 1982). Because of this physical interpretation, there are a large variety of mathematical expressions based on the wave equation for different geometries. The most simple, and the one used in this thesis, is achieved by applying a linear moveout to the data and summing over amplitudes (Yilmaz, 1987), where p = s in#/V is the ray parameter or slowness, V is the velocity and 9 is the angle between the ray and the vertical. From a geometrical point of view, this transformation approximately maps lines and hyperbolas in data space to points and ellipses in model space, respectively. The L R T has various applications, the chief of these being predictive deconvolution because multiples preserve periodicity with offset in this domain. Furthermore, the whole sequence of seismic processing, from stack to migration, can be applied in the linear Radon domain (for laterally homogeneous media) (Dunne and Beresford, 1997). (1.23) Chapter 1. Introduction 11 K Figure 1.4: The Fourier projection theorem (from Stewart, 1991) The linear RT is usually computed using an inversion approach but, for regular sam-pling, it can be computed very efficiently using the ID and 2D FFTs. The Fourier projection theorem states that the 2D Fourier transform of an image or medium can be obtained from the ID Fourier transforms of the projections (Stewart, 1991) and, as a result, the image or medium can be reconstructed from projections. Figure 1.4 shows a scheme from Stewart(1991) of the LRT obtained by using first a 2D F F T ofthe data and then a ID inverse F F T . 1.4.2 P a r a b o l i c a n d P s e u d o - h y p e r b o l i c R T The parabolic RT (PRT) and pseudo-hyperbolic RT (PHRT) are widely used to sepa-rate primaries from multiples based on moveout differences. The PRT transformation represents the data set as a summation of parabolic basis functions. To make seismic reflections more similar to parabolas two different transformations are commonly applied to the data, normal moveout (NMO) correction (Hampson, 1986) or a time stretching Chapter 1. Introduction 12 [t —> t2 transformation (Yilmaz, 1987)]. In the first case, after the normal moveout cor-rection, those events that have a stacking velocity close to the velocity trend used in the N M O correction are mapped close to zero curvature parabolas, while those that are not, map to short hyperbolas in the model space. In the second case, hyperbolic events become parabolic so a perfect hyperbolic reflection maps to a single (point) or, more precisely, a wavelet. A different approach for decomposing hyperbolic events in time invariant basis func-tions (that can be computed in the frequency domain) is the pseudo hyperbolic transform (Foster and Mosher, 1992). The integration is performed along curves that depend on the offset h and a particular depth zre{ Each one of these methods has advantages and disadvantages, mostly depending on the data and goals. The normal moveout correction sometimes is not a good enough approx-imation to achieve high resolution (reflections map to hyperbolas), the time stretching transformation in the data distorts the wavelet for shallow times, and the pseudo hy-perbolic R T is time invariant only if the depth is fixed to one single value of interest. However, the three methods can be applied almost without change in the R T code, and they can be applied to the same data set to see which one produces better results. 1 . 4 . 3 H y p e r b o l i c R T The H R T , often termed velocity stack, can be obtained with a summation along hyper-bolas (Thorson and Claerbout, 1985) (1.24) (1.25) Chapter 1. Introduction 13 where h is the offset, m(r, V) is the Radon transform, V represents stacking velocity, and r is two way travel time. This operator (equation 1.25) is equivalent to an N M O correction with velocity V , followed by stacking. It is very large but also sparse. Its total size is (nt x nh) x (nr x riV), but it only has (nt x nh x nV) nonzero elements. Typical sizes of these parameters are nt « 1000, nh nV 100. Often a hyperbolic RT is preferred, in particular when events with very different and large moveout appear very close in time (not uncommon in deep water marine data). In this situation the N M O correction does not completely flatten some of the hyperbolic events, and the parabolic RT maps them to hyperbolas rather than points in the Radon panel. In these hyperbolic RT computed in the time domain produces better results. The HRT is not commonly used because, as already mentioned, it takes consid-erably more time to compute than time invariant transforms. Nevertheless, computation of the RT in the time domain can be speeded up considerable by the utilization of itera-tive solvers like C G with sparse matrix computations (Trad et al., 2001b) and semblance analysis to estimate a mask in the Radon space where the solution is focused (Marfurt and Nemeth, 2001). 1.5 F r o m processing to inversion: the sparse R T Processing and inversion are often difficult to separate. Usually we think of "transforms" as signal or image processing tools, while inversion defines the process of obtaining the model from the data gathered in a physical experiment. Also, inversion is usually asso-ciated with the minimization of a cost function that penalizes undesired features of the model and differences between observed and predicted data. If, however, a transform is defined in terms of an operator that, acting upon a model, produces the data, then the transform can be thought of as the process of inverting the data to obtain such a model. Chapter 1. Introduction 14 To obtain a simple picture we follow Claerbout (1992) to say that a forward Radon operator, L , produces time shifting and change of color (frequencies) on a given model, m, to produce the data, d, d = L m . (1.26) Another operator, the adjoint L H , has the ability to undo the time shifting but not to recover the original color, producing the adjoint model m a = LHd. (1.27) To recover the original frequencies we have to deconvolve the action of LHL as the adjoint model, m a , is related to the original one by m a = LHLm, (1.28) Let us pause at this last equation to think of the Radon transform as a signal processing tool. The matrix L ^ L can be arranged as a circulant matrix. Therefore L f f L m is a circular convolution of a single column of the matrix L ^ L with the model m and the Radon adjoint model is simply a filtered version of the unknown true model. Most problems associated with Radon transforms are associated with the non-invertibility of the smearing filter, L ^ L , as it has to be deconvolved from the adjoint model to recover the true unknown model. For large aperture, and fine sampling in offset and velocities, the obtained models will have (almost) no loss of resolution. The recovered model, m r , will be similar to the true model, m, m r = ( L H L ) _ 1 L H L m , (1.29) and data will be recovered completely according to d = L ( L H L ) _ 1 L H d . (1.30) Chapter 1. Introduction 15 W i t h insufficient data (which is always the case) the problem is underdetermined and an infinite number of solutions exist for each input gather. A cost function must be specified to select the more appropriate solutions among all the possible ones. From a sta-tistical point of view (Tarantola, 1987), inversion consists of finding the model that, while maximizing the probability of its predicted data (minimizing the data misfit), has a prob-ability density function in agreement (in shape and dispersion) with our prior knowledge or assumptions (Ulrych et al., 2001). This prior information selects a particular solution among the infinity of models that produce the same data. Under this philosophy, the main problem is how to include in the mathematical for-mulation the statistical parameters that define the likelihood of the data and the prior for the model. These include the model and noise covariances, C m and Cd respectively, the mean model and noise. Assuming a null initial model, the optimization problem is cast with the following cost function (Sacchi and Ulrych, 1995a) = m f f C ; 1 m + (Lm - d) f f C r f "' (Lm - d) (1.31) where C m , covariance of the model, is a diagonal matrix for uncorrelated models with elements given by the assumed variance of the model. C d , covariance of the residuals, is a diagonal matrix for uncorrelated noise with elements given by the variance of the data. The matrix C m depends in general on the model, and on the statistical parameters, shape and standard deviation, of the assumed distributions for the model parameters. In the next chapter this model covariance matrix will be the subject of further discussions, but for now it is enough to say that by choosing a prior pdf with thick tails, the cost function will penalize less heavily large elements in the model, and a sparse model rather than a smooth model will be favored. The minimization of <b gives rise to two expressions for the model m, the overdeter-mined and the underdetermined version. In most of this thesis, I have worked with the Chapter 1. Introduction 16 overdetermined version, because it has some computational advantages when using iter-ative reweighted least squares. By working with one frequency at a time, the model that minimizes this cost function (the maximum a posteriori solution in Bayesian teminology) is m = (L^C^L + C-1)-1^ -^^ (1.32) The resolution matrix, that connects the true model with the computed model, is (Tarantola, 1987) R = (L^C^L + O ^ C A (1.33) In general, the effects of the covariance matrices are, 1 - The data covariance matrix modifies the calculated model by down weighting the information given by components of the data with large variance. 2 - The model covariance matrix modifies the calculated model by favoring elements of the model associated with large variance. In other words, CTO acts as a focusing tool. The filter LHL is called the smearing filter because it represents the loss of resolution by information leaking among the elements of the model. The smearing filter can be computed as dj = YI Ahkexp(-ibjg(hk)(pi - pA) (1.34) k where g{hk) is the particular function of offset used (h for L R T , h2 for P R T ) . Note that when we formulate the Radon transform as an inverse problem, the initial geometrical problem of mapping curves to points and vice verse, becomes a mathematical problem of expressing the functions that exist in the data space with a limited set of basis functions that are spanned by the model space. If we have a parabolic event in the data space whose curvature is not exactly available in the range of possible curvatures spanned by the Radon space, this event will be represented using a combination of a large number Chapter 1. Introduction 17 of the basis functions that are available. As a consequence, aliasing may occur (Wiggins and Abma, 1990). Likewise, we can expect that complex information present in the data will be represented with a combination from elements of the model space, thus loosing the desired property of compact support. The adjoint model obtained by applying the Radon operator backward and forward on a data set is equivalent to applying a smearing filter in the true Radon model. The effects of finite aperture and sampling reflect in the impossibility of exactly deconvolving the smearing filter because of its band limited nature. A band limited filter results in loss of information, and its deconvolution requires some inversion algorithm with regularization. The purpose of regularization is to supply, somehow, the missing information in the data by assuming some characteristic of the model. Plainly explained, the regularization used in the sparse R T supplies the kind of information that leads to higher resolution in the model space. Theoretically, this information arises from assuming a smoothness in the missing data rather than in the unknown model. Figure 1.5 gives a very simple picture of the sparse R T . The data in general are incomplete because there are gaps or truncations (question marks in the Figure). The underdetermined nature of the problem requires the incorporation of some assumptions regarding the data inside the question marks. We can safely assume that the missing data is likely similar to the available data (smooth data). This assumption implies sparse models. If we seek to separate events in the model space, then a sparse model is clearly superior in terms of our objectives. A smooth model and a sparse model will honor the data where there is information, but the gaps will be treated differently. When mapping back to the data space a sparse model will preserve continuity of events and a smooth model will preserve low energy of events, i.e., it will fill the gaps with zeros. Chapter 1. Introduction 18 t t Figure 1.5: Differences between smooth and sparse R T . 1.6 T ime versus frequency domain When the curves are time invariant the R T can be solved by decomposing the large problem into several smaller problems, one per frequency. When the curves are time variant this is not possible and the problem has to be solved globally. For comparison, the forward transform for a data set of size nt x nh requires an operator of size (nt x nh) x (nr x np). Instead in the frequency domain nt/2 inverse problems each one of size nh x np must be solved. At first it seems that computing the R T in the time domain is much inferior to computing it in the frequency domain. However there are many reasons to work in the time domain. Time variant curves can better match the natural events and therefore can produce parsimonious models (see Figure 1.6). In later chapters we will see that results in the time domain are cleaner because of our ability to impose sparseness Chapter 1. Introduction 19 Figure 1.6: Time invariant vs time variant R T in time as well as curvature. In short, the model and data weights are more flexible. O n the other hand, frequency domain operators better preserve the waveforms, and above all, are much faster since the solution of several small problems is much faster than that of one large problem. 1.7 Scope and contribution of this thesis This thesis addresses many problems that appear when the sparse R T is applied to real data, and compares many (or most) implementations of the R T published in the geophysical literature. I introduce also some new implementations of the R T and discuss applications for noise attenuation, interpolation and migration. Chapter 1. Introduction 20 Original contributions from this thesis are, to the extent of my knowledge, the im-plementations for the time domain RTs, the irregularly sampled Radon space, the sparse elliptical R T , the hybrid R T using two different sets of basis functions, the implementa-tion for the sparse R T in the frequency domain using weighted conjugate gradients with general cross validation for subspaces, the considerations of aliasing in the R T for irreg-ular sampled data and the merging of R T and Equivalent Offset Migration for multiple attenuation-migration, velocity analysis and velocity corrections. Partially new topics are the elimination of aliasing by using sparse R T in the time domain and sparse decon-volution with hybrid inversion to remove coherent noise. As a result of this work I have developed several programs, written in C , preserving the Seismic Unix format from the Colorado School of mines. This format allows their immediate use for real data processing. The developed programs have been used for all the examples of this thesis, and are fully adapted to work with real data. These programs include among other general modules, standard R T , sparse R T in the frequency and combined time-frequency domains, both with a fast C G - F F T algorithm, time domain hyperbolic and elliptical RTs , hybrid R T with linear and hyperbolic basis functions and the Radon E O M migration algorithm in the frequency and time domains. C h a p t e r 2 Invers ion for sparse models : a lgor i thms a n d examples As a signal processing tool, the Radon transform is often studied in terms of filtering and mapping. Because most implementations of the R T involve inversion, however, we will deal in this chapter with regularization, cost functions and operators. I will present some basic concepts behind sparse inversion and some numerical algorithms commonly used in these problems. The complexity of the two dimensional R T transformation makes the understanding of the effect of sparse inversion on the data and model spaces rather diffi-cult. Therefore, I will first present several one dimensional examples of signal processing problems that are attacked by means of sparse inversion techniques. 2.1 S o m e basics about inverse theory Most physical problems deal with a transformation from model to data space, d = L m , (2.1) where L is an operator that contains the physics of the experiment and generates the data d from the physical model m. The geophysical problem is often the reverse of the physical experiment, i.e. given the results of some experiment we need to invert the operator to find the physical model that gave rise to the data. To invert operators, including the Radon transform, we first define a cost or objective function. This cost function is a mathematical expression for the undesired characteris-tics of the unknown model, usually an expression for the differences between the model 21 Chapter 2. Inversion for sparse models: algorithms and examples 22 response and the data, that we try to minimize. The cost function is problem depen-dent but, in general, we pursue a model that honors the data and has a minimum of information not required by the data. This statement of goals is commonly presented as minimize |Wm |£, (2-2) subject to ||Wd(d-Lm)||9 = 0d, where <bd is some estimation of the noise level in the data plus the error due to the deficiencies of the proposed model to explain the data. is a matrix of data weights, often a diagonal matrix, containing the inverse of the standard error of the data. Wm is a matrix of model weights that we can tune in order to enhance our preference in the model as, for example, resolution or smoothness. and ||.||^ are the norms that we use to measure the size of the resulting mathematical quantities (upper subscript indicates a power, lower subscript indicates the norm). By using two different symbols, p and q, I indicate that different norms can be applied for the misfit and the model. The choice of the norms and model weights is related to the desired characteristic of the solution. Some norms heavily penalize large elements in the model, leading to models with many small elements (smooth models) rather than few large components (sparse models). This is the case for the commonly used £ 2 norm. O n the contrary, the t\ norm is not sensitive to the size of the individual elements but only to the total size. For this reason, it is often used to obtain sparse models. Likewise, the model weights are powerful tools to favor a particular type of solution, but some prior knowledge about the problem is required. The model that fulfills our quest can be found by setting these measures of failure ("anti-goals") in a cost function and finding the model that is obtained by minimization, J=||Wd(Lm-d)||9 + A||Wm ||P, (2.3) Chapter 2. Inversion for sparse models: algorithms and examples 23 where A is a tradeoff hyperparameter (often a Lagrange multiplier), that assigns differing relevance to one term or another. The minimization of the cost function in equation (2.3) leads to the following expres-sion for the model, m = ( A W £ W m + L T W j W r f L ) - 1 L r W j W d d . (2.4) 2.1.1 T h e B a y e s i a n a p p r o a c h The link between prior information and inversion is well explained in Bayesian terms (Ulrych et al., 2001), where inversion consists of estimating the most probable model, i.e., the model that maximizes the probability density function (pdf) of the model given all the available information. Bayes's theorem provides a method of computing this pdf by using the likelihood of the data and any other prior information available. If a prior pdf for the model is known before the experiment then, following the experiment, the updated (pdf) becomes Pfmld) = (2.5) This equation represents the posterior information concerning the problem at hand. P ( d | m ) is the likelihood of the data, P ( m ) is the prior information in terms of the pdf of the model and P(d) is the pdf of the data. The maximum a posteriori ( M A P ) of P ( m | d ) defines the most probable model. When the pdfs are expressed in terms of negative exponentials (e.g. as in the case of the generalized Gaussian pdf), then the M A P solution also minimizes the following cost function Jp = m T C - 1 m + ( L m - d ) T C , 1 ( L m - d ) , (2.6) Chapter 2. Inversion for sparse models: algorithms and examples 24 where and C m are covariance matrices associated with the data and model spaces respectively, that are either known or assumed. Thus, the M A P solution is m = ( C " 1 + L T C d " 1 L ) - 1 L T C d - 1 d . (2.7) The data covariance matrix becomes very important when dealing with outliers and, in seismic processing, with bad traces and offset truncations. Its function is that of a preconditioner of the data space. Likewise, the model covariance matrix plays an extremely important role in the design of high resolution operators. The following two diagonal model covariance matrices have been proven to be very convenient for the estimation of sparse models: [ C m ] M = cT^ + m z 2 , (2.8) which is obtained when we assume a Cauchy pdf f(m) = ( T r c x ™ ) " 1 ^ + p - ) " 1 for the model, and [Cm]i,i = 2(jm\mi\, (2.9) which arises by assuming an exponential l\ pdf, f(m) = 2~ 1cr~ 1 exp(—\m\/o~m). A small threshold e is included to assign a constant weight to small samples. W i t h this change the model covariance matrix becomes [Cr. = < 2crm ImA if xxii > e 1 ' (2.10) 2crm|e| otherwise In the last two expressions, the parameter o m represents some assumed background dispersion measure or scale parameter for the pdf of the model, and acts as a hyper-parameter that allows a trade-off between the degree of fitting and sparseness of the solution [see discussion in Thorson and Claerbout (1985), and Sacchi and Ulrych, (1995)]. Chapter 2. Inversion for sparse models: algorithms and examples 25 By simple comparison we can associate the weighting matrices from equation (2.3) with the model and data covariance matrices of equation (2.6), (2.11) Cd 1 = W j W , . In this way, we can use a Bayesian methodology to compute the model weights. The opposite is also true, we can take some commonly used model weight operators and un-derstand the implications from a statistical point of view. For example, the l\ norm is equivalent to the use of the £ 2 norm and diagonal model weights with elements propor-tional to [W m ] , , = - L , (2.12) because IHIi = £ K I = m T W ^ W m m = | | W m m | | 2 . (2.13) i These model weights result also from the factorization of the model covariance matrix for a double exponential distribution (equation 2.9). Other common model weights are differential operators, used to obtain smooth mod-els. These operators result from the factorization of non-diagonal model covariance matrices, implying correlation between the model elements. Large variations between neighboring model parameters are penalized, therefore the model becomes smooth. 2.2 H y p e r p a r a m e t e r s The question that naturally arises is, once we know the shape of the model covariance, equations (2.8) or (2.10), or their equivalent expression in terms of weighting matrices, how do we choose the deviation or dispersion <rm? This topic is indeed difficult because Chapter 2. Inversion for sparse models: algorithms and examples 26 we do not have realizations of the model to perform sampling statistics. One method is to perform computation with different hyper-parameters to obtain the residual-model norm curve, commonly called the L curve, (see for example Hansen, 1998). However, for the seismic Radon transform problem, such a strategy is too expensive and impractical. A simple empirical method is preferred because it is not possible to afford expensive methods for the estimation of the hyperparameters and, at the same time, the choice of them is not critical for the results. This is a common characteristic of all image-inversion problems. When working in the time domain, only one value for is required for the whole gather, and the solution is not very sensitive to this parameter because the iterative algorithm is stopped much before obtaining a complete solution. In the frequency domain, one value of is required for each frequency, therefore it is important to choose a robust measure of the size or scale of the model parameters. Robust estimators are, for example, those based on quantiles of the distribution. A simple procedure is to normalize the model vector by dividing by some quantile, and then to choose the same cr^ for all frequencies. When we expect to obtain a very sparse transform because the data are very similar to the basis functions of the decomposition, then o2m can be set very small. If er^ is set too small, however, the residuals will become larger as a consequence of forcing excessive sparseness. O n the other hand, a large cr^ parameter implies a non-sparse transformation or, equivalently, a Gaussian least squares solution. When increasing o2m, the solution first resembles the standard non-stochastic Radon transform result, and then the adjoint solution. Hence, simple inspection of the data should be enough for the initial estimate. If the residuals are large, o2m should be increased. If the R T looks very similar to the standard R T , then probably we can obtain more sparseness by decreasing o2m. In my experience, a good practical approach is to set cr^ to a very small value, and then to use the number of external iterations as a tradeoff parameter (see the next section for Chapter 2. Inversion for sparse models: algorithms and examples 27 a definition of external iterations). By using this approach I have been able to compute most of the results presented in this thesis with a similar value for o2m. How do we find the dispersion of the data o\ necessary for the data covariance matrix? The dispersion of the data can be found by means of a statistical approach (Sacchi, 1995). We can also find the variance of the data at high frequencies, where the seismic signal consists of noise only. We will see later that an efficient approach is to use the number of internal iterations for an inner linear solver as a tradeoff parameter for o\. The number of iterations can been estimated by a tolerance parameter or a general cross validation routine for subspaces methods (Haber, 1997), explained in the next chapter. 2.3 M e t h o d s for sparse invers ion 2.3.1 T h e I terat ive ly R e w e i g h t e d Least Squares m e t h o d ( I R L S ) A model-dependent covariance matrix or a weight matrix leads to a non-linear optimiza-tion problem. The problem is iteratively solved by fixing the model at every iteration and solving a linear system of equations. In each one of these "external" iterations, the model weights or model covariance matrices are updated and a simple linear solver is applied. The inner linear problem can be also solved iteratively by applying "internal" iterations. This is known as the Iteratively Re-weighted Least Squares (IRLS) algorithm (Scales et. al, 1988). If the £2 norm is applied for both the misfit and model constraint, then the cost function becomes J = | | W d ( d - L m ) | | 2 + | | W m m | | 2 , (2.14) which is known as a quadratic cost function. When a different norm is desired, it is often possible to design the model and data weight operators to cast the problem in quadratic form, for example as shown before in equation 2.13. Cost functions with quadratic form Chapter 2. Inversion for sparse models: algorithms and examples 28 are simple to solve by quadratic programming. A further simplification can be applied by transforming Equation (2.14) to the stan-dard form, J=||(d-Lm)||2 + ||rh||2, (2.15) where L = WdLW"1, m = Wm , and d = Wd . (Hanke and Hansen, 1993). Such a transformation is trivial when dealing with sparseness constraints because the model weighting matrices are diagonal. For more complicated weighting functions, for example when imposing smoothness, W"1 is substituted by the pseudo-inverse. In this way, the model weights matrices are incorporated into the forward and adjoint operators. Note that the previously linear operators become nonlinear operators. The minimization of the general cost function in equation (2.14) leads to the following expression for the model, m = (W^Wm + ^ WjW.Lr^ WjW.d , (2.16) but in the standard form, equation (2.15) the model is m= (AI + LTL)_1Lrd. (2.17) Note that the system of equations that results from minimizing the cost function in (2.14) is equivalent to the following extended system of equations WdL Wdd m = 0 (2.18) whereas the system of equations resulting from the cost function in (2.15) is equivalent to L d m = AI 0 Chapter 2. Inversion for sparse models: algorithms and examples 29 If a subspace method, like a conjugate gradient algorithm is used, then equation (2.19) can be solved by setting A = 0 (no regularization) and letting the number of internal iterations play the role of regularizer. In this case, the inner problem is partially solved because the conjugate gradient algorithm is stopped before the solution is complete. A C G algorithm will produce a solution inside the span of the conjugate vectors, known as the Krylov space. These vectors are associated with the singular values of the transformation, through the Ritz polynomials (Hansen, 1998). When the algorithm is stopped before completion, only some of the Krylov vectors, those corresponding to the performed iterations, contribute to the solution. This approach will be successful if the contribution from the Krylov vectors associated with the small singular values of the transformation can be neglected in the solution. If these vectors contribute significantly, a better approach would be to obtain a complete C G solution, where the choice of the hyperparameter A will be essential for a good solution. A complete C G solution also dampens the contribution from the poor resolved singular vectors, but allows some partial contribution from them. Among other possible choices, the number of iterations performed can be controlled by a Generalized Cross Validation Criterion ( G C V ) . A non-costly G C V computation for subspace methods (Haber, 1997) has proven to be very efficient. The iterative algorithm is stopped when the G C V function reaches a minimum. The G C V function for the solution of d = L m can be computed as GCV(P) = ^ krl\;> „ (2.20) trace(l - Ck,k)2 where Ck,k are the elements of the matrix resolution matrix C = C(B) = L ( L T L + / 3 W ^ W m ) " 1 L T . (2.21) Chapter 2. Inversion for sparse models: algorithms and examples 30 If a subspace method is applied, the matrix C can be approximated by (Haber, 1997) Ek-i \d(8)k - dibs]2 GCV{ittr) = ^ / _ \ t e r ) : ] . (2.22) This last approximation does not involve any extra operation inside the C G algorithm other than a simple quotient per iteration. In summary, the sparse inversion with the Iteratively Reweighted Least Squares (IRLS) algorithm is performed in this thesis as: • Transforming the cost function from the general form to the standard form. • Applying an iterative process to minimize the cost function. The model weights are computed using the result from the previous iteration, and a linear minimization algorithm is carried out at every iteration. • Solving the inner linear problem by a subspace method, where the number of iter-ations plays the role of the regularizer. A variant of this approach is to eliminate external iterations and use model weighting functions estimated from other sources, for example, a semblance analysis or a previous C D P gather. When solving non-linear problems with the IRLS method, a linear search is generally required to check convergence. However, for the kind of problems presented in this thesis, convergence is (at least in theory) always achieved. For the double exponential distribution the model constraint is equivalent to the use of the l\ norm . Therefore, the cost function is a combination of an l2 measure of misfit and an h measure of model structure. Alliney and Ruzinsky (1994) prove the convergence of the IRLS method for this kind of problem. For the Cauchy distribution, a proof of convergence is given by Sacchi (1995). In my experience, the problem is locally convergent, and that either the Chapter 2. Inversion for sparse models: algorithms and examples 31 adjoint model or the least squares solution with uniform model weighting appear to be good starting models. The cost function is computed after every external iteration to make sure that the new solution gives a smaller cost function. Between 3 to 5 iterations are generally enough, and every new iteration just adds some sparseness to the solution while trying to honor the data. However, too much sparseness can produce undesirable results, and the optimal minimum of the cost function is often, for real data, not a good solution. It is important to note that, in this approach to the Radon transform, we are generally not interested in achieving the minimum solution for the optimization problem, but rather in improving the results obtained by the simpler linear approach. The first iteration is the zero order regularization least squares result. After that, every new iteration increases resolution by favoring a sparse model. The problem converges to a sparse solution for perfect hyperbolic data, but there is a limitation for real data on how sparse the model can be. Perfect sparseness in the R T domain only arises when working with perfect synthetic hyperbolas of constant amplitude. Therefore, for real data, we always stop before reaching the optimal solution. In summary, the number of external iterations are usually better estimated by visual inspection of the data, whereas the number of internal iterations is efficiently estimated by using an automatic stopping criterion. 2.3.2 The conjugate gradient algorithm Given a system of equations Ax « b , a solution for x can be obtained by minimizing the £2 square norm of the residuals ||b — Axil?,. If A is a symmetric positive definite matrix (SPD), then the conjugate gradient algorithm (CG) can be applied, requiring only the result of the forward operator A on x. If the matrix A is not S P D , then a least squares Chapter 2. Inversion for sparse models: algorithms and examples 32 method can be used, Ar Ax = Arb, (2.23) which again results from minimizing the l i square norm of the residuals |b — Ax||2. The conjugate gradient least squares algorithm (CGLS) is simply a reformulation of C G for the least squares problem to avoid explicit computation of ArA, and to take advantage of the matrix of coefficients being factored. It is also known as conjugate gradient for the normal residuals ( C G N R ) . Another method, the conjugate gradient algorithm for the normal equations ( C G N E ) , consists of applying a change of variables, After solving the problem for u, the solution is x = ATu. In both cases we need the forward and adjoint operations, A x and A T b respectively, where the adjoint operator is the transpose of A . The drawback of these two methods is that the condition number of As a result of the transformation to standard form explained before, the operator A can be modified by weighting operators to decrease the influence of very large residuals or to increase sparseness of the solution. In the same manner a preconditioner acting either on the data or model space can be applied to decrease the condition number of A and to arrive more quickly to a useful solution. When preconditioning is applied, the eigenvalues of the resulting system are better clustered than those of the original system, and therefore the iterative solution converges faster than without preconditioning. In C G L S , the solution is built in the Krylov subspaces, xC9 = x0 + K(ATA, ATr0) = x0 + span{ATr0, ATAATr0,(ArA)m_1ATr0}. (2.25) Every term in the solution is left multiplied by the adjoint operator A , which is unable to incorporate information from the nullspace. Hence, C G is an efficient regularizer because x = ATu, such that AATu = b. (2.24) ATA or AAT is the square of the one for A. Chapter 2. Inversion for sparse models: algorithms and examples 33 it excludes the nullspace from the solution in a fashion similar to that in obtaining a minimum norm solution. However, when the preconditioning acts on the model space, the nullspace of the transformation can be altered. The C G solution is free from the nullspace of A W " 1 but not necessarily from the nullspace of A . Thus, by preconditioning, we can add prior information to the model that lives in the nullspace of the transformation (Nichols, 1997). Different C G L S implementations arise from different ways of manipulating the weight-ing matrices to keep the symmetry required for C G methods (Saad, 1996). In this thesis I have used several variants of C G . In the frequency domain I have used mainly the right preconditioned version, d = A W - 1 W m m , (2.26) whose least squares solution is W " T A T d = W - T A T A W " 1 u . (2.27) In the time domain, I found very useful the left preconditioned conjugate gradient residual method, M _ 1 A T d = M ^ A ^ A m , (2.28) where M = W m W ^ . The asymmetry in the matrix M - 1 A T A does not affect the application of C G because the Euclidean dot product is replaced with a M-product (Saad, 1996). The systems (2.27) and (2.28) have the same eigenvalues (Shewchuk, 1994) and hence both preconditioned schemes have the same convergence. Two small advantages of the left preconditioned system are that, first, we do not require to factor the model covariance matrix and, secondly, less operations are required. Let us summarize the advantages of the C G method: • After every iteration an approximation to the true solution is obtained. Therefore a fast solution can be obtained by stopping the algorithm before completion. Chapter 2. Inversion for sparse models: algorithms and examples 34 • The solution at every iteration is approximately similar to the solution obtained by truncated singular value decomposition (SVD), using as many different singular values and vectors as iterations have been performed in the C G . Therefore, an automatic regularization is performed by limiting the number of iterations because the largest singular values contribute to the solution first. • Every element used to build the solution has a null contribution to the nullspace so the solution is regularized. The reason is that the nullspace of the transformation does not map onto the Krylov subspace. • The only required operations are by means of forward and adjoint operators. These operators can be applied as pure operators or as sparse matrix vector multiplica-tions. • Among all the iterative methods, C G is the best known and its convergence prop-erties are very well understood. 2.4 Example 1: Sparse Deconvolution This is an example of a signal processing problem addressed as an inverse problem with a sparse constraint. The purpose of this example is to show the transformation of the singular vectors of the model space when applying model weight operators. A general singular value decomposition (GSVD) will allow a look at the nature of the model space before and after model weighting. Let us consider the problem of the deconvolution of a known wavelet. Figure(2.1) shows the a) wavelet, b) the reflectivity series, c) their convolution (data) and d) the data plus noise. Because of the band limited nature of the wavelet and noise in the data, the problem is ill conditioned and requires regularization. The reflectivity series is sparse Chapter 2. Inversion for sparse models: algorithms and examples 35 (a) Ricker wavelet (b) Reflectivity 0.02 0.04 0.06 0.08 0.1 time(sec) (c) data 0.2 0.4 0.6 time(sec) (d) data+noise 0.2 0.4 0.C time(sec) 0.2 0.4 0.6 time(sec) Figure 2.1: (a) Ricker wavelet, (b) reflectivity (c) convolution of the wavelet with reflec-tivity (clean data), (d) data plus noise and we will see that it is important to incorporate this knowledge into the solution. I will use a model covariance function derived from a Cauchy pdf. The example will show that a more parsimonious solution can be obtained by updating our knowledge about the model in every iteration. Parsimony in the solution is important because there is a limitation, given mainly by the noise, on the number of singular vectors that we can use to build the model. Therefore, parsimony leads not only to a more resolved model but also allows us to decrease the influence of noise in the solution. 2.4.1 Deconvolution in terms of inversion The easiest way to deconvolve the wavelet from the data is by division in the Fourier domain, (2.29) Chapter 2. Inversion for sparse models: algorithms and examples 36 Here T is Fourier transform (FT) and T~x is the inverse F T , d is the data vector, m is the desired reflectivity and w is the wavelet. This approach works well only for a full band wavelet and data without noise. A more robust approach is the least squares solution, _! T(w)*T(d) m — T (2.30) {T{w)*T{w) + eY where * signifies complex conjugate and e is a small number to prevent division by zero. Let us cast the problem in terms of a system of equations, A m . (2.31) where m is a vector containing the desired reflectivity, and d is the seismic trace. A is a circulant matrix whose columns are shifted versions of the wavelet. If we denote the wavelet by w, the shape of this matrix is ^ iu(0) w(n- 1) w(n-2) ^ w(l) w(0) w(n — 1) iu(0) (2.32) yw(n—l) w(n — 2) . J The multiplication of this matrix with the reflectivity vector is exactly equivalent to the circular convolution between w and m. By padding the wavelet with zeros the product A m is equivalent to the discrete convolution between w and m. Circulant matrices have two very interesting properties (Hansen, 1999): • A l l circulant matrices of order n have the same eigenvectors Wj, with elements W / j = exp(—i2nlj/n) for l,j, = 0, ...n — 1 . (2.33) Chapter 2. Inversion for sparse models: algorithms and examples 37 • The eigenvalues of circulant matrices are given by the F T of the first column, in this case of the wavelet. Therefore, it is possible to build an eigenvector solution with very few operations. The eigenvector decomposition of A = W(JT(w))WI/, where WH denotes the Hermitian of W . The eigenvector solution for m is m = A~xd = W(^r(w))-1Wi/d. (2.34) This equation is exactly equivalent to equation(2.29), which shows that deconvolution of two vectors by Fourier transforms is equivalent to a S V D solution for the equivalent sys-tem of equations. Thus, if a system of equations has a circulant kernel, then the solution can be found very efficiently by fast Fourier transform ( F F T ) techniques. Note also that band limited deconvolution, spectral division of only some frequencies, is equivalent to truncated S V D , a construction of the solution using only some of the singular vectors. In this problem, the benefit of working with the circulant matrix instead of the wavelet comes in terms of the casting of the deconvolution as an inverse problem. A cost function can be formulated as J=||Wd(Am-d)||^ + ||Wm |H, (2.35) where m is the desired model and d is our data vector. The matrices W M and contain the model and data weighting functions respectively. 2.4.2 SVD and GSVD In this section I will show the effect of the model weighting matrix on the singular vectors of the model space. To keep simplicity in the discussion the data weighting function will be set to a simply constant diagonal matrix, W J , = aj1!, where I is the identity matrix and <TJ is the variance common to all the data. Chapter 2. Inversion for sparse models: algorithms and examples 38 (a) Initial model 1 0.5 I ..A L A J 0 0.5 1 , I 0 20 40 60 80 100 120 140 (b) Solution for exact data: SVD and 127 singular vectors 2, . 1 1 1 1 1 1 0 20 40 60 80 100 120 140 (c) Solution for noisy data: SVD and 29 singular vectors '0 20 40 60 80 100 120 140 Figure 2.2: (a) True reflectivity (b) Truncated S V D solution for noise-free data, (c) T S V D for noisy data Given a matrix A, its singular valued decomposition (SVD) is A = USVT, (2.36) where U contains the singular vectors for the data space, V contains the singular vectors for the model space and S is a diagonal matrix containing the singular values of the matrix. A regularized solution for m that minimizes the cost J = A||m||^ + ||d- Am|2 (2.37) can be obtained by the singular vector expansion of m (Hanke and Hansen, 1993) as p u7d. m = E^ 2)^V i + £ tfdV,, (2.38) i=l ° l i=p+l where n is the size of d, p is the rank of A, CTJ are the singular values of A, the (p((Ti) are scalars between 0 and 1 that limit the contribution of the singular vectors in the rank Chapter 2. Inversion for sparse models: algorithms and examples 39 of A to the solution. The second term in the summation gives m*, the nullspace of the kernel A , whose dimension is n — p. The filter functions (p(o~i) can be computed as = ~T~7~T • (2-39) °i + A A simpler approach is the truncated S V D that uses only those singular vectors asso-ciated with the singular values above some numerical threshold. In this case < (^<7j) = 0 or 1, and the expansion is limited to k ufd m E (2.40) i=i ^ where for example k can be chosen as rank(A). Even though these two solutions are different, both minimize the l2 model norm for the model space subject to the data constraints. The difference is that in one case we minimize the model norm in the model space spanned by all singular vectors, and in the other case, the minimum model norm solution is in the model subspace span for only the singular vectors we use. Figure 2.2a shows the true reflectivity (real solution) and Figure 2.2b is the solution obtained with all eigenvectors but one. The solution is almost exact because there is no noise in the data. Figure 2.2c gives the T S V D solution when noise is added to the data and we take only 29 singular vectors (out of 128) are used to avoid the dominance of noise in the solution. The number of singular vectors chosen for both cases were determined by minimizing a Generalized Cross Validation function ( G C V ) . These very different results can be explained in terms of the Picard condition (Hansen, 1998). The Picard condition states that in order that there exists a square integrable solution the singular values must decay less quickly than the projection of the data onto the data singular vectors. For example if we use an infinitely small discretization of the data space, the Picard condition is satisfied if i = N f u r d ) N p I U J O J < ^ ( 2 4 1 ) i=0 ai Chapter 2. Inversion for sparse models: algorithms and examples 40 (a) o°o OCO oo°n ° £ , HT,.?,?. x,v,.** * ***>k. |uTb|/<T. Figure 2.3: (a) Exact data (b) Noisy data This can be understood by noting that the norm of the least squares solution is i=N /„rd.\ 2 l | m w | I I = E ( ^ • (2-42) i=0 This quantity will be very large for small cr; unless d is completely inside the range of A, i.e., u f d j | < o{. (2.43) Figure 2.3a shows how the magnitude of the projection of the data onto the rank of A falls below the magnitude of the singular values for most of the rank of A for exact data, but, for only part of the rank of A for noisy data (Figure 2.3b) . The projection (ufdj) is sometimes called the Fourier coefficients. If the data are contaminated by noise, the Picard condition is only true for the first few singular vectors (Figure 2.3b). Because of this, a truncated S V D has to be built only with the first few singular vectors, and resolution is very poor, as can be seen in Figure 2.2c. Figures 2.4a-d shows four singular vectors for the given kernel. We can see that the singular vectors for a circulant matrix are particularly inappropriate to build a spiky Chapter 2. Inversion for sparse models: algorithms and examples 41 20 40 60 80 100 120 20 40 60 80 100 120 Figure 2.4: Four of the first singular vectors for the model space of the given kernel model, because they are harmonic functions. Therefore what we need is a transformation of the problem that changes the shape of the singular vectors to make them more similar to the desired model. We will see that, by using a l\ or Cauchy pdf as a prior, the eigen-vectors become much more suited for the problem at hand. Such a transformation of the problem can be accomplished in various ways. First I will present a Generalized Singular Value Decomposition ( G S V D ) , that allows a look at the transformed singular vectors. The G S V D (Hansen, 1998) is a way of analyzing the general problem of decomposing a pair of matrices (A, L) as U 0 0 In—p X"1 (2.44) L = V M 0 X"1 (2.45) In this decomposition A is m x n, and L is n x n with rank{L) = p. The columns of U G nmxn and V e Wxp are orthonormal, UTU = I and VTV = I; X 6 Tlnxn is Chapter 2. Inversion for sparse models: algorithms and examples 42 (a) Generalized eigenvalues Y Figure 2.5: Generalized singular values for the convolution kernel. non singular with columns that are ATA orthogonal (they diagonalize the matrix ATA) and £ and M are p x p diagonal matrices with elements Oi and pi respectively. These elements are non-negative and ordered as 0 < <rx < • • • < ov < 1 (2.46) and They are normalized such that 1 > Pi > • • • > PP > 0. (2.47) of + p2 — 1, for i = 1, • • •, p. (2.48) The generalized singular values 7 $ of (A, L) are Oi 7 i = — , for i = 1, Pi (2.49) Figures 2.5a-c show the generalized singular values and the elements cr- and pi for the given kernel A. Chapter 2. Inversion for sparse models: algorithms and examples 43 (b) 0.25 0.2 0.15 0.1 0.05 0 -Figure 2.6: Four generalized singular vectors for the convolution kernel and the sparse constraint. Similar to S V D , the G S V D of a pair (A, L) provides new basis functions to diagonalize both matrices in the pair, r £2 0 XrATAX = 0 In—p and XTLTLX M2 0 0 0 If L is square and nonsingular (as in this deconvolution problem), (2.50) (2.51) AL"1 = U [S M"1] VT . (2.52) This equation shows that the S V D of the matrix quotient AL-1 is closely related to the G S V D of (A, L) as the S V D singular values of the quotient are the generalized singular values of the pair and the singular vectors U and V are the same. Later I will use the S V D of the matrix quotient to obtain a sparse result. Moreover, I will invert the operator Chapter 2. Inversion for sparse models: algorithms and examples 44 100 120 0.5 0: -0.5 - f K f A A I - i X imoiiyiiiini V -20 40 60 80 100 120 Figure 2.7: (a) The solution constructed by using the first 27 singular vectors (no prior information), (b) The solution obtained by adding the last 9 generalized singular vectors (second iteration of IRLS). (c) The same but after the third iteration of IRLS. AL-1 by using a conjugate gradient algorithm, and show that the solution is close to the G S V D solution of the pair (A, L). The truncated G S V D solution is obtained by means of a standard form transformation (Hansen, 1998) A = AL^ , where is the pseudo-inverse of L. Applying T S V D to A, A = ALt = U £ 0 0 I n—p X"XX M-1 0 V1 (2.53) A = U M-1 0 V 7 and r £ 0 0 In—p If Z is the approximation of A obtained with the k largest singular vectors we have v (2.54) Zk= ua^J i=p—k+\ (2.55) Chapter 2. Inversion for sparse models: algorithms and examples 45 (a) (b) 20 40 60 80 100 120 20 40 60 80 100 120 i = 3 i =4 Figure 2.8: First four singular vectors for the modified convolutional kernel A L . The truncated G S V D solution is mk = L^ZJ^d mL,fc = J2 — * + £ UIdx* (2-56) i=zp-k+l i=p+l Let us now examine the generalized singular vectors when the matrix L is incorporated into the cost function to favor sparse models. Figure 2.6 displays the four last singular vectors for the convolution kernel (in G S V D the singular values increase with the order, so that the importance of singular vectors increases with the order as well). Comparing these basis functions with the ones displayed in Figure 2.4, a more parsimonious model is clearly possible. If we were able to use all the singular vectors a similar model will be obtained , but as the Picard condition showed, noise in the data allows the use of only the first few singular vectors. Figure 2.7a shows that the first 27 singular vectors are not enough to build a sparse solution, whereas in Figures 2.7b-c we see that the generalized singular vectors do a much better job (second and third iteration, respectively). The sparseness constraint L is improved at every iteration. In fact it is not necessary to use G S V D for these results. Because of the connection Chapter 2. Inversion for sparse models: algorithms and examples 46 Figure 2.9: Evolution of the solution built by using the modified kernel A L 1 as the iteration proceeds and L is more precisely computed. between the S V D A L " 1 and the G S V D ofthe pair ( A , L ) , all we need to do is to transform the cost function from the general form J = m T W ^ W m m + ( A m - d)T W j W d ( A m - d) to the standard form with J = mTm + (Am — d)T (Am — d) (2.57) (2.58) A = W , A W - \ m = W m m , d = W d d . (2.59) Figure 2.8 displays the first four singular vectors of the modified kernel and Figure 2.9 shows the solution built by using the modified kernel A L - 1 . As L depends on the previous solution, the first iteration is just equivalent to the standard form because there is no prior information to compute L and a constant diagonal matrix is used. As the iterations proceed, the model weight matrix L provides a better prior and the singular vectors have become more similar to the true model. Chapter 2. Inversion for sparse models: algorithms and examples 47 (a) Initial model 1 | 1 1 ' r 1 1 1 1 1 1 1 0 20 40 60 80 100 120 140 Figure 2.10: The solution built by using Preconditioned C G . (a) Initial model, (b) Solu-tion by solving the expanded system of equations. 2.4.3 Preconditioned CG The S V D gives us a comprehension of how the solution is built by the weighted basis functions of the model space, and how prior information changes the characteristics of the basis functions. However, S V D has a high computational cost because it performs a whole diagonalization of the matrix of coefficients in terms of orthonormal vectors, which sometimes is more than required. Hence, for large problems it is advantageous to use iterative algorithms. I will address the deconvolution problem again, using the same general procedure (IRLS). However, this time I will use a more efficient solver, a C G algorithm, for the in-ner linear problem. The data are, as before, the convolved reflectivity series plus noise. In Figure 2.10 we see the solution obtained with the full expanded system. Figure 2.10a dis-plays the original model, and Figure 2.10b shows the solution obtained by expanding the system of equations with the model weighting function, and computing the full solution Chapter 2. Inversion for sparse models: algorithms and examples 48 zero order regularization 1 0.5 X X x X -0 0.5 -1 X X -20 40 60 L = l 80 100 120 1 0.5 0 0.5 -1 i i i i i A ft X X . X 20 40 60 L * 1 80 100 120 1 0.5 - 1 1 A1 i X < j \ X X 0 0.5 -1 X Y 20 40 60 80 100 120 Figure 2.11: The solution built by using left preconditioned C G . (a) First iteration, L = I because there is not prior information yet. (b) fifth iteration, (c) Tenth iteration. at every iteration to this well conditioned problem (equation (2.19)). Figures 2.11a-b-c shows the evolution of the system as the iterations proceeds. By comparing these results with the G S V D results (Figure 2.7) we see that they are very similar. Therefore, I can conclude that weighted C G produces similar results as G S V D . I will use weighted C G in most of this thesis. 2.5 E x a m p l e 2: D e n o i s i n g b y sparse h y b r i d i n v e r s i o n As a second example of a signal processing problem attacked by sparse inversion, I will consider the separation of signal and coherent noise. A similar problem will be considered in Chapter 5, to develop a hybrid R T . The example will serve to explain an alternative to the IRLS algorithm, the log barrier method, that is useful for sparse inversion. Noise is understood in general as the part of the data that is not predicted by a proposed model. Sometimes noise can be modeled as well, and then the difference between Chapter 2. Inversion for sparse models: algorithms and examples 49 signal and noise resides in our understanding or definition of the physical process under consideration. In this section I will discuss a simple method of the separation of signals that possess different properties. The method addresses separation of data and noise by using hybrid inversion plus signal or noise modeling. Hybrid inversion means that the connection between model and data is not limited to one single class of operators or basis functions, but rather a combination of models, each one tailored to generate one particular feature of the data. By using only one of the operators, the signal or the noise can be modeled after the inversion. The different sets of basis functions are not required to be uncorrelated, which makes the process very flexible. Freed of this constraint, we rely on our inversion solver to produce an intelligent choice of bases for the different features of the model. This choice is specified by the following simple set of instructions to the algorithm: if a feature of the data can be explained by different sets of basis functions at the same time, use only the one that is more parsimonious. In other words, find the basis the produces the most sparse model. Sparse inversion is thus essential to this process. I will use a mixed cost function, with an £ 2 norm for the data misfit and an l x for the model norm. A l\ norm for the data misfit can equally be applied, and this is the key idea behind robust methods (insensitive to outliers), because the weight or leverage imposed by large residuals on the solution is smaller for the t\ norm than for a standard l 2 norm. In this example I will compare three common solvers, preconditioned C G with IRLS, non linear C G and a log barrier interior point method. I will present brief explanations of the IRLS method for the hybrid problem and of the log barrier method, but I will omit details of the non linear C G since it is well explained in the literature (see for example Shewchuk, 1994). Chapter 2. Inversion for sparse models: algorithms and examples 50 L i Lis (2.62) 2.5.1 P r o b l e m Given a signal or data set, the problem is to obtain a decomposition into several sets of basis functions at the same time. For example, a simple two component hybrid model is, d(x) = zZWiW + £ & A - ( z ) , (2-60) * j where d(x) are modeled by the basis functions gi(x) and hj(x) with weights given by aj and bj respectively. The weights can be computed by a correlation process between the data and the model basis functions. In matrix formulation the problem posed is, d = Lim.i + L 2 m 2 , (2.61) or equivalently r m i m 2 where the forward operator is [Li L 2 ] , and the data vector d is modeled by the basis functions contained in the columns of operators Li and L 2 with weights given by the model vectors m x and m 2 . The hybrid model can be found by satisfying the following statement of goals, minimize | |mi | | i and | | m 2 | | 1 , subject to | |d — L^rri! + L 2 m 2 | | 2 = (f>, (2.63) where </> is some estimation of the noise level in the data plus failure of the proposed model to explain the data. Once the hybrid model has been obtained, if the signal and noise map to different models, then modeling can be used to separate them, d = L i m i . (2.64) The problem is then, how to force the signal and noise to separate? One of many possible recipes is: use basis functions similar to the signal, and basis functions similar to the noise and apply sparse inversion. Chapter 2. Inversion for sparse models: algorithms and examples 51 2.5.2 T h e i terat ive re-weighted least square m e t h o d for the h y b r i d invers ion The goals in equation (2.63) can be put together in a £2 —1\ cost function J = | | W d ( L i m i + L x m 2 - d ) | | | + | |mi | | i + | | m 2 | | i , (2.65) where is a matrix of diagonal weights containing the inverse of the standard deviation of the data. The previous cost function can be cast as a £ 2 — ^2 cost function by the introduction of model weights W m i , W m 2 as mi m lWrn2]U = -1=, (2.66) such that l|m||i = £ K | = m T W ^ W m m = | | W m m | | 2 . (2.67) i Hence, the cost function becomes J = | | W d ( L i m i + L i m 2 - ci)Hi + | | W m l m i | | 2 + | | W m 2 m 2 | | 2 . (2.68) We first transform equation (2.65) to the standard form, J = IKLimi + L i m 2 - d ) | | 2 + \\nnWl + \\rh2\\l (2.69) where L i = W d L i W " 1 ! , L 2 = W d L 2 W ~ 1 2 , m i = W m i m i , m 2 = W m 2 m 2 and d = Wdd. In this way, the model weight matrices are incorporated into the forward and adjoint operators, that now become nonlinear operators. Chapter 2. Inversion for sparse models: algorithms and examples 52 To simplify notation, let us define L — [Li L2], (2.70) m = (2.71) m 2 and 0 (2.72) 0 Wm2 By using this notation, the minimization of the general cost function in equation (2.68) leads to the same expression in equation (2.16), and in the standard form, equation (2.69) leads to the expression in equation (2.17). The IRLS method can be used to solved both of these equations. 2.5.3 L o g B a r r i e r i n t e r i o r p o i n t m e t h o d Sparse inversion can often be carried out by minimizing the l\ norm of the model subject to the data constraints, a problem commonly addressed by linear programming tech-niques. The goal of linear programming is to determine the values of decision variables that maximize or minimize a linear objective function, where the decision variables are subject to linear constraints. It is a particular case of a general constrained optimization problem. Following Chen et al., (1998), the linear program in standard form is a constrained optimization problem defined in terms of a variable x e R m where cTx is the objective function and Ax = b is a collection of constraints, and x > 0 is a set of bounds. min cTx subject to Ax = b, x > 0 (2.73) Chapter 2. Inversion for sparse models: algorithms and examples 53 The collection of feasible points [x : Ax = b, x > 0] is a convex polyhedron in R m called the simplex. A common method for linear programming is also called the simplex method. The first step in the simplex method is to find one feasible solution (the constraints are fulfilled). Then, iteratively the method jumps from one vertex of the simplex to an adjacent vertex where the cost function is smaller. This is equivalent to taking a subset of the kernel, solving that subsystem of equations, subsequently discarding one of the equations and replacing it by another not in the subsystem. The final solution is the one that minimizes the cost function. The simplex produces feasible solutions in every iteration, which means that it follows the borders of the simplex. The number of operations required grows exponentially with the number of variables. Karmarkar (1984) developed a method that requires a number of operations bounded by a polynomial in the number of variables. As result of this work, a whole family of algorithms arose, called interior point methods (Chong and Zakh, 1996). Interior point methods start from some initial point well inside the simplex and go through the interior of the simplex. Since the solution of a L P problem is always in a vertex of the simplex, the solution approaches the boundary of the vertex as the iterations proceed. These methods are based on "Duality theory". Following Chen et al., (1998), the primal linear program is the linear program in standard form min cTx subject to Ax = b, x > 0. (2.74) The related dual linear program is min bTy subject to ATy + z = c, z > 0. (2.75) x is called the primal variable; y and z are called the dual variables. The quantity |b — Ax||2 is the primal infeasibility, and |c—z —ATy||2 is the dual infeasibility. The difference between the primal and dual objective is call the duality gap. The variables (x, y, z) solve Chapter 2. Inversion for sparse models: algorithms and examples 54 the linear program if and only if the primal infeasibility, the dual infeasibility and the duality gap are all close to zero. A primal-Dual Log Barrier LP method The problem min ||x||i subject to Ax = b, (2.76) can be solved with a linear program like min cTx subject to Ax = b, x > 0. (2.77) where the variable x is the positive version of the variable x, i.e., x = [x+ (—x-)] and A = [A (-A)] The problem min ||x||l + A||Ax-b||| (2.78) can be solved with a perturbed linear program like min cTx + 1 /21[p>|[ 1 + 1 /2|-yx|| subject to Ax + <5p = b with x > 0 and 5 = 1. (2.79) Perturbed linear programming is really quadratic programming, but it preserves the structure of linear programming (Chen et al., 1998). 2.5.4 Results The deconvolution problem addressed in the previous section becomes more interesting with the introduction of coherent noise. A simple regularization approach is not enough to remove or attenuate this noise, simply because its contributions to the solution are not necessarily small for large eigenvalues, because of its coherence. One way to separate signal and noise in this problem is just to map all the data to a different (model) space Chapter 2. Inversion for sparse models: algorithms and examples 55 (a) Ricker wavelet (b) Reflectivity 0 0.1 0.2 0.3 0.4 0 0.2 0.4 0.6 0.8 time(sec) time(sec) Figure 2.12: Basis functions and models, (a) Signal basis, (b) True model, (c) Noise basis, (d) Model for noise where the signal and noise live separately. This is possible if a set of basis functions can be found for both. Once they are separated in the model domain, the separation in the original domain is obtained elegantly by modeling with only one of operators. I will present two examples with different coherent noise but without random noise. Even without random noise the problem is difficult, because it highly under-determined and ill conditioned (condition number of the Hessian matrix is about 1 0 1 9 ) . From now on, all figures will have a similar structure. As a general guide, in all the results a continuous line will represent predicted data or computed model and crosses or dots will represent given data or true model. In every figure, the plot at the top will show the recovered models, [mim 2 ] , as a continuous line, and the (unknown) correct models as a dotted line. The first half of the vector is m i and the second half is m 2 . The middle plot will contain the (unknown) true signal (dotted line) and the predicted signal (continuous line). The plot at the bottom will contain the same information for the coherent noise. Chapter 2. Inversion for sparse models: algorithms and examples 56 E < -1 --2\z i _ 0 0.05 1 — a . 1 " i o a. e < -1 --2h i _ 0 0.05 1— -a i o Q. E < -1 --2h i _ 0 0.05 Figure 2.13: (a) Signal, (b) Noise, (c) Data The forward model is d = L i i r i ! + L 2 m 2 , (2.80) where L i is the same circulant matrix from the previous example, constructed with the vector w, a Ricker wavelet. The second operator, L 2 , will be the generator of coherent noise. In the first example the noise will be generated with blocky functions. This noise could occur, for example, as a consequence of variations in the D C level of an electronic recording equipment. Figure 2.12a shows the wavelet or basis function for the signal, Figure 2.12b is the reflectivity series or model for the signal (unknown), Figure 2.12c shows the basis function for the coherent noise and Figure 2.12d is the model for the coherent noise (unknown). We expect to obtain the clean signal as shown in Figure 2.13a, and the coherent noise as in Figure 2.13b, from the given data in Figure 2.13c. It is important to start by noting that prior information about the model is essential in this problem. The prior information here is used as a weighting function for the model. (a) Signal \ I I \ Chapter 2. Inversion for sparse models: algorithms and examples 57 zero order regularization X X X ] X 1 x 1 -X X X ,x I 1 x A I 1 r i 11* * " * u Figure 2.14: (a) Zero order regularization (b) IRLS with correct model weights (perfect prior), (c) IRLS with perfect model weights for the noise, (d) IRLS with model weights obtained from the data. Figure 2.14 shows four plots, in the top panel we see the results obtained by ignoring any sparseness constraint, and trying to fit the data with a smooth model. The second plot shows that by using the true model for the prior produces a perfect result in one iteration (similar to using a perfect preconditioner). The third plot shows the result when prior information about the coherent noise is available. The bottom plot is a more realistic approach, there is no model information available at the beginning, but it can be bootstrapped from previous iterations. Figure 2.15 shows the results of IRLS. Each external iteration has basically two steps, the model weights are recomputed each time, and a linear system of equations is solved. The linear solver is a left-preconditioned C G algorithm for the extended system of equa-tions. The linear inner system was solved almost completely avoiding over fitting by using an appropriate hyperparameter. The panel at the top of Figure 2.15a shows the inverted model (full line) and the true model (dots). The recovered amplitudes and locations are Chapter 2. Inversion for sparse models: algorithms and examples 58 -0.5 - I -1 -1 -0 -1 -1 -0 --1 -Figure 2.15: Results for IRLS with W T C G L S (blocky noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line). correct, however there are a few spikes not recovered. As can be seen in the middle and bottom plots, the signal and coherent noise can be separated very well. Let us now see how the log barrier algorithm deals with the same problem. This method has six parameters to set, the primal, dual and duality gap tolerances, the hy-perparameters 7 and 8, and the number of iterations in the inner system of equations. As a consequence, there is a large flexibility in the settings allowing different results. When these parameters are correctly tuned, it is possible to obtain very good results very rapidly. The primal, dual and duality gap tolerances are set to small values, be-cause they give a measure of the precision required in the solution. 7 is also set to a small number because it gives a tradeoff between a £2 — h solutions (barrier term). 8 is set to 1, for obtaining a finite misfit (noisy data). The number of iterations in C G have to be set large enough to obtain some precision in the calculated steps. Figure 2.16 shows the results with this algorithm. The model has been inverted Chapter 2. Inversion for sparse models: algorithms and examples 59 Log Barrier Noise= 0.0 % , required time 14.67 s 50 100 150 200 250 signal and predicted signal 0 20 40 60 80 100 120 noise and predicted noise 0 20 40 60 80 100 120 Figure 2.16: Results for Log Barrier method (blocky noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line). successfully, slightly better than before, with less computation time. The amplitudes are correct and the signal and noise can be separated almost perfectly. Some of the small spikes in the reflectivity have not been recovered, but this seems to have little effect in the prediction of signal. The choice of the hyperparameters 5 and 7 changes the solution completely and probably better results can still be obtained. For these results the computing time was almost half that for IRLS, but very good solutions were also computed in about one fourth of the time for IRLS. The nonlinear C G ( N L C G ) algorithm produced good results as well (Figure 2.17), but more time was required, initially about a factor of ten. By using the Hessian matrix as a preconditioner, however, the computing time became similar to IRLS. For larger problems this preconditioner can be computed only approximately, for example by using incomplete L U factorization. The problem with N L C G is the slowness of the line search. Other implementations with faster line search algorithms needed less computing time but Chapter 2. Inversion for sparse models: algorithms and examples 60 Nonlinear C G : Noise= 0.0 % , required time 25.37 s 1 1 1 . . r 1 ; ; ' ' • 50 100 150 200 250 signal and predicted signal \y 0 20 40 60 80 100 120 noise and predicted noise i i i I 1 1 \*T?T» / ^ 0 20 40 60 80 100 120 Figure 2.17: Results for nonlinear C G method (blocky noise test) with preconditioning. Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line). did not give good results. The solution contained some leakage between signal and noise models and some extra spikes in the coherent noise model. A line search with bracketing (Press et. al., 1992) gave also poor results. The line search for N L C G is very sensitive. For this problem this method ranked with very low efficiency, at least compared to the other two approaches. Chapter 2. Inversion for sparse models: algorithms and examples 61 Figure 2.18: Basis functions and models, (a) Signal basis function, (b) True model, (c) One of the basis functions for the coherent noise, (d) Model for noise Chapter 2. Inversion for sparse models: algorithms and examples 62 (a) Signal 2 -1 I _ _ 0 0.05 0 1 0.15 0.2 0.25 0.3 (b) Noise 0.35 0.4 0.45 0.5 \ rA A A /V A rA A / 0 0.05 0 1 0.15 0.2 0.25 0.3 (d) Data+noise 0.35 0.4 0.45 0.5 \ i J -0 0.05 0 1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 time(sec) Figure 2.19: (a) Signal, (b) Noise, (c) Data The second test is the same problem but the coherent noise has harmonic basis func-tions (Figure 2.18). This kind of harmonic noise could be found for example in electro-magnetic data contaminated by power line fields. Many times I have found this problem in the magnetotelluric method because the data are often gathered close to populated areas, and the electromagnetic fields generated in power lines are much stronger than the weak natural electromagnetic fields generated in the earth. Figure 2.19 shows the signal, the noise and the data for this example. The results for this test are excellent using both the IRLS and the log barrier al-gorithms (Figures 2.20 and 2.21 respectively). Computation times with the log barrier algorithm were twice as fast as with the IRLS method, although this difference may be different for larger size problems. Chapter 2. Inversion for sparse models: algorithms and examples 63 IRLS wpcgnr: Noise= 0.0 % , required time 27.19 s I I I 1 I : 1 i . l 1 1 i 1 : - ! 50 100 150 signal and predicted signal 200 250 -_ i • 1 0 20 40 60 80 noise and predicted noise 100 120 0 20 40 60 80 100 120 Figure 2.20: Results for IRLS with W T C G L S (harmonic noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line). Log Barrier Noise= 0.0 % , required time 9.93 s signal and predicted signal 0 20 40 60 80 100 120 noise and predicted noise 0 20 40 60 80 100 120 Figure 2.21: Results for Log Barrier method (harmonic noise test). Top subplot: real model (dots) and inverted model (line). Middle subplot: real signal (dots) and predicted signal (line). Bottom subplot: real noise (dots) and predicted noise (line). Chapter 2. Inversion for sparse models: algorithms and examples 64 2.6 E x a m p l e 3: R e - s a m p l i n g w i t h the F o u r i e r t r a n s f o r m ( F T ) The last example illustrating sparse inversion consists of the regular resampling of an irregularly sampled time series by inversion with the F T . Given a vector d obtained by sampling a continuous function, we can think of its elements as the data produced by the action of the inverse Fourier transform operator F on some unknown model m , where m is actually the F T of these data. Accordingly, the F T of the data can be obtained by applying the inverse operator, F - 1 , to the data. For regular sampling in the time and frequency domains, frequency interval Af — l/(NAt), and fmax = —fmin = 1/(2 A i ) , the inverse operator becomes the adjoint of F . Because of this, we always think of the F T as an orthogonal operator. However, when these conditions are not fulfilled, for example if the sampling is irregular, the F T of the data has to be obtained by a different approach. In the discrete Fourier transform ( D F T ) , the difference between inverse and the adjoint operator is neglected, except through introduction of the sampling interval as a weighting factor. As a result, the F T becomes very poorly resolved and contains numerical artifacts. In the Least Squares method, on the contrary, the inverse operator is either approximated by the pseudoinverse or a regularized inverse. Thus, the F T becomes The hyperparameter e provides a zero order regularization. If the condition number of F ^ F is not very large, e can be reduced to zero, and the inverse operator is exact. This situation happens when the sampling and aperture is carefully chosen to fulfill some stability criteria (Schonewille, 2000). By replacing the constant weighting with model dependent weights we can obtain a more localized F T and eliminate some computing artifacts along the frequency axis. If the true model is sparse, this method is superior d = F m (2.81) m = ( F H F + e I ) _ 1 F H d . (2.82) Chapter 2. Inversion for sparse models: algorithms and examples 65 to the simple LS approach. Otherwise, both methods produce similar results. In all these methods, the data can be resampled to a regular grid by the action of the forward operator (equation 2.81). In Figures 2.22 and 2.23 we can see the differences between all these cases. Figure 2.23a shows a time series (solid line) constructed by the addition of two simple harmonics. The amplitude spectrum of the F T of these data, calculated by simple multiplication of the data with the adjoint operator, is shown on Figure 2.22a. The predicted data (dots in Figure 2.23a) fit the original data exactly. In Figure 2.23b the same data (solid line) have been irregularly sampled, and applying the D F T produces the amplitude spectrum in Figure 2.22b, whose predictions onto a regular grid do not fit the data very well since the model is not correct. A LS solution, Figure 2.22c, does a better job, and the data are predicted fairly well. Nevertheless, as a result of the natural sparseness of the model, there is a large improvement by using the sparse least squares F T (Figures 2.22d and 2.23d). 2.7 Conclusions Many seismic processing techniques can be cast as inversion problems, i.e., to find the model that under the action of a particular operator produces the given data. Because the characteristics of the earth are almost never completely determined by the data, partly because of the physics of the experiment, partly due to the incompleteness of these mea-surements, we face the problem of non-unique models. To deal numerically with this problem we use regularization, avoiding in this way division by zero in our computations. However, regularization is certainly much more than just avoiding division by zero. It represents our knowledge about the missing information. By simple zero order regular-ization we adopt a non compromised position, where all the model parameters have the Chapter 2. Inversion for sparse models: algorithms and examples 66 same probability. This is sometimes the most we can say about the problem, and this approach is, therefore, correct. On the other hand, many times we do know something about the missing information, or we want a solution with certain characteristics. Here is where the model covariance matrices, model weights functions, spectrum of the model, or whatever mathematical entity we use to explain our view of the world, play an im-portant role. This additional information can be incorporated into the solution by using different numerical approaches. I have shown here a few methods, preconditioned C G , transformation to the standard form, generalized singular value decomposition and linear programming. In fact, the particular numerical technique that we use is not nearly as important as the details of the prior information and how it is expressed in mathematical terms. The results, as I have shown in these examples and will show throughout this thesis, often merit the extra work required that incorporating prior information requires. In the first example, that of deconvolution, we have seen how the singular vectors for the preconditioned model space acquire a resemblance to the model. This allows us to span the solution by using a few singular vectors and the noise contribution to the solution is in this way diminished. In the second example I presented a technique to separate signal and coherent noise, based on a combined or hybrid inversion. By using sparse inversion the models for signal and noise are well separated. If the kernels that we use for the models are orthogonal, a more sophisticated combination of basis functions can be applied, because the necessary inverse operators are very simple to compute. This is the approach taken by the technique called "basis pursuit" (Chen et al., 1998). This hybrid method will play an important role later in this thesis in the developement of a R T that allows the separation of linear and hyperbolic events. These simple examples have shown clearly that sparse inversion can be a necessity rather than luxury. Chapter 2. Inversion for sparse models: algorithms and examples 67 F i g u r e 2.22: Four i e r t r ans fo rm of i r r egu la r ly s ampled da ta , a) F T of r egu la r ly s a m p l e d da ta , b) D F T of i r r egu la r ly s ampled da ta , c) Leas t squares F T of the same da ta , d) Sparse F T of the same da ta . Chapter 2. Inversion for sparse models: algorithms and examples 68 (a) Pred ic t ion f rom F F T (reg s a m p l e d data) 0 0 .05 0.1 0 .15 0.2 0 .25 0 .3 0 .35 0.4 0 .45 0 .5 (b) Pred ic t ion by D F T 0 0 .05 0.1 0 .15 0.2 0 .25 0 .3 0 .35 0.4 0 .45 0.5 (c) Pred ic t ion by L S F T 11 i i i ! i i i i i i I 0 0 .05 0.1 0 .15 0.2 0 .25 0 .3 0 .35 0.4 0 .45 0.5 (d) Pred ic t ion s p a r s e L S F T J I L " 0 0 .05 0.1 0 .15 0.2 0 .25 0 .3 0 .35 0.4 0 .45 0 .5 Figure 2.23: Fourier transform and irregular data, a) true signal with regular sampling (line) and predicted by F T (dots), b) irregularly sampled data (line) and predicted data by D F T (dots), c) data (line) and predicted signal by L S F T (dots), d) data (line) and predicted signal by sparse F T (dots) Chapter 3 T i m e invariant R T Time invariant Radon transforms use basis functions with a general expression t = T + pg(h), where g(h) is a function of the offset, h, that defines the shape ofthe curve and p is some curvature or slope parameter. As a consequence, the implementation for all time invariant RTs is similar, with the only change being in the function g(h). Because of the time invariance, the transform can be calculated independently for every frequency. Consequently, a very large problem is decomposed into a series of small problems (Figure 1.6). In general, inverse problems require a number of computations of G(MK), where M is the number of unknowns and k is around 2 or 3. Therefore it is always faster to solve one small problem per frequency than a big problem for all frequencies (or times) simultaneously. The ability to perform computations per frequency is not the only computational advantage of time invariant curves. Operators designed in the frequency domain have the nice property of preserving very well the inherent waveforms. On the contrary, time domain methods are in general difficult to implement in such a way that preserves the waveforms, requiring filtering and interpolation. This is similar to what happens in seismic migration and dip moveout (DMO) methods, where frequency domain algorithms produce in general better images. On the other hand, time domain methods have more flexibility for handling time and space variant problems. Although both the major applications of the RT, as a tool to separate multiples from primaries and as a tool to fill gaps, are enhanced by increasing resolution due to prior 69 Chapter 3. Time invariant RT 70 information, the impact on multiple removal is more evident. In this chapter I address different implementations of time invariant RTs and show applications for synthetic and real data. The emphasis of this discussion will be the possible improvements obtained with the sparse R T as compare to the standard R T . 3.1 S y n t h e t i c d a t a e x a m p l e s This section explains and compares different R T implementations. I use the parabolic R T because of its applicability to multiple removal, but linear and pseudo hyperbolic RTs share the same implementations. For every implementation there is one synthetic example for multiple removal to test resolution, and one example for interpolation. The test data sets are gathers displaying hyperbolic curves convolved with a Ricker wavelet. A normal moveout correction (NMO) partially flattens the hyperbolas before the R T is applied and the inverse N M O moves them back to their original moveout. These moveout corrections are not shown in the examples. The synthetic gather for the resolution tests consists of four events following a given velocity trend, simulating primaries, and three events with the velocity of the first pri-mary, simulating multiples. Similar to what happens in surface reverberations, one of these multiples has its polarity reversed. The destructive interference of this event with a primary located at the same two way travel time produces almost zero amplitude at near offsets and a change of shape of the wavelet as the different moveouts increase with offset. In all the resolution examples I apply an automatic muting in the Radon space to eliminate the primaries, and from this filtered space I recover the multiples. These multiples are subtracted from the data, as is usually done in standard seismic processing. If the transformation has a non negligible nullspace, this procedure leads to less distortion Chapter 3. Time invariant RT 71 ofthe signal than direct removal of multiples. Events with residual moveout smaller than 0.910~ 8 s/m 2 (the units represent residual moveout per offset square) have been removed by the mute process. This procedure is intended to make a fair comparison, but in general the results could be improved by applying a careful mute following a data-dependent polygonal curve. The synthetic gather for the interpolation tests is very similar but I have increased the moveout differences between primaries and multiples. These large moveout differences make the test more difficult because the N M O correction will flatten the primaries only, leaving the multiples with considerable moveout. As a consequence, the multiples after N M O will still depart considerably from a parabolic shape and their mappings will not be naturally sparse. 3.2 Non sparse R T Inversion with zero order regularization can be thought of as a deconvolution of the spreading function, i.e., the impulse response of the forward and inverse transforms. The spreading function matrix, L H L , has a Toeplitz structure (for regular sampling in the Radon space). Therefore a non informative prior leads to a very fast implementation by applying the Levinson algorithm, which requires a number of operations of 0(N2), where N is the number of unknowns. This is the standard R T , available in most seismic pack-ages, widely applied in the industry to remove multiples. The implementation consists of taking the Fourier transform of the data, and for every frequency solving the Toeplitz system of equations given by ( L " L + eI)m = L " d , (3.1) where the constant regularization factor e is usually taken as the percent of the noise in the data. Chapter 3. Time invariant RT 72 Figure 3.1: (a) Synthetic gather, (b) Non-Sparse RT. (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a). Chapter 3. Time invariant RT 73 Figure 3.2: (a) Synthetic data, (b) Synthetic gather with gaps, (c) Non-Sparse RT. (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a). Chapter 3. Time invariant RT 74 Figure 3.1a presents the synthetic gather for the resolution tests, Figure 3.1b shows the R T of the data, and Figure 3.1c its spectrum. The primaries and multiples are not well resolved in the R T space, they overlap. By applying a mute in the Radon panel, an approximate separation is obtained. Figure 3.Id shows the reconstruction of the multiples, and Figure 3.1e shows the residuals, i.e., the desired primaries obtained by subtracting panel (d) from (a). The interpolation test reveals another negative aspect of the smooth R T , perhaps less serious than low resolution. Figure 3.2a presents another synthetic data set used for the interpolation test. I have created three large gaps in this gather, as can be seen in Figure 3.2b. The R T of the incomplete gather exhibits a clear separation of the primaries and multiples (Figure 3.2c) because of the large moveout differences between them, but several truncation artifacts appear as well. The data predicted from the R T space tends to fill the gaps, but the amplitude of the interpolated events decreases within the gaps, as Figure 3.2e demonstrates. The residual panel, Figure 3.2f, obtained by subtracting gather (e) from gather (a), reveals clearly this artificial and incorrect amplitude variation with offset (AVO) . These problems, truncation artifacts and interpolation with low amplitude, arise be-cause, of the several models that fit the data, we have chosen a smooth rather than a sparse fit. These models are equivalent in terms of reproducing the data available, but not in terms of the unavailable information within the gaps. Where there is no information, smooth models correspond to zero amplitude traces (minimum energy), whereas sparse models, because of the nature of the R T , correspond to traces with the same amplitude and character ofthe neighboring traces. This is similar to the difference between standard Fourier and maximum entropy spectra. Zero order regularization enforces smoothness in the model, when actually the smoothness is desired in the transition between data and gaps, in the data domain. Furthermore, smooth solutions obviously lead to poorly Chapter 3. Time invariant RT 75 resolved models. In geophysical applications, smooth models are often a good choice, but not in this case. Smooth models have great flexibility to adapt to variations in the data, which is an advantage, but where the data do not impose any constraint, the predicted data tend to zero providing minimum energy smooth models. 3.3 S p a r s e R T In the statistical approach, the mentioned problems of poor resolution and interpolation are a consequence of assuming that every model is as probable as any other, i.e., a non informative covariance matrix. Even when we do not know the shape of the true covari-ance matrix, we can impose our goal of sparseness by using an informative covariance matrix, or equivalently, weight functions leading to sparse models. Using a sparse con-straint during the inversion produces higher resolution transforms that deal much better with the two mentioned problems. Sparse models obviously are easier to separate, as they occupy a smaller area in the model space. Also sparse models are consistent with smoothness in the data space, so smooth amplitude variations between data and gaps are favored rather than abrupt jumps in gaps. The price that one pays for an informative model covariance matrix is computation time. From a computational point of view, a model dependent covariance or weight function is punished twice. First, the Toeplitz structure in the system of equations is destroyed by the summation of the non constant diagonal matrix C m . Hence, Levinson recursion is at first not possible. The consequence of this can be appreciated by noting that a more general solver, like Cholesky decomposition, requires C ( M 3 / 6 ) operations per frequency, where M is the number of traces in the Radon space (Levinson requires only 0(M2)). Roughly speaking, if 3 minutes are required using Cholesky decomposition, about 3 seconds would suffice using Levinson recursion. Secondly, if the information Chapter 3. Time invariant RT 76 in the model covariance matrix or weighting function is bootstrapped from a previous estimation of the model, then the computation time increases even more. The linear problem becomes nonlinear and iterations are required, so computing times are increased by about 3-5 times. Also, the need for iterations turns out to be a problem when there is insufficient aperture in velocities, as the R T increases unnecessarily at the boundaries of the model space when trying to fit data that map outside of the R T space (Wiggins and Abma, 1990). Fortunately, instead of using a direct solver of order 0(M3), we can use an iterative method such as conjugate gradients or some other subspace method. 3.3.1 W e i g h t e d conjugate gradients least squares ( W T C G L S ) This is a simple approach that has proven to be quite efficient and stable. It is equivalent to applying left and right preconditioning to the forward modeling equation d = L m , (3.2) that becomes W d d = W d L W " 1 W m m . (3.3) A different explanation can be given in terms of a transformation of the corresponding cost function (Hansen, 1998). The cost function in general form that arises by imposing model dependent model weights is transformed to the cost function in standard form [see equations (2.15) and (2.17) ], and the hyperparameter A multiplying the identity matrix is set to zero. This produces the system of equations d = L m , (3.4) with L = W ^ L W " 1 , m = W m m , and d = Wdd. The least squares solution is LHd = L f f L m . (3.5) Chapter 3. Time invariant RT 77 Figure 3.3: (a) Synthetic gather, (b) Sparse RT obtained with the W T C G L S method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a). Chapter 3. Time invariant RT 78 Figure 3.4: (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse RT obtained with the W T C G L S method, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a). Chapter 3. Time invariant RT 79 By putting back the original variables we obtain the system W - T L r W j W d d = W - ^ W j W . L W - 1 ™ , (3.6) that must be solved at every external iteration. This system of equations is equivalent to that resulting from the least squares solution of equation (3.3). As the regularization hyperparameter A is set to zero, the number of iterations plays the role of zero order regu-larization, while the model weights, incorporated into the forward and adjoint operators, forces the sparseness. This method is usually called Weighted Conjugate Gradient Least Squares ( W T C G L S ) . This implementation is very stable and produces excellent results. There are two reasons for this, 1- The hyperparameter for the model space acts by multiplication of the forward and adjoint operators, rather than by addition. As a result, the action of the model weights depends only on relative rather than absolute amplitudes. 2- The number of iterations in the C G algorithm is estimated by using Generalized Cross Validation ( G C V ) for subspace methods (Haber, 1997), that stops the iterations at the minimum of the cross validation function, i.e., when the contributions to the solution start to fit mostly local noise rather than consistent global information from the data. The G C V function is usually expensive to compute, involving several matrix multiplications, however this approach for subspaces is very simple and fast. It exploits the connection between singular values and Krilov spaces. Figure 3.3 shows a) the data, b) the R T , c) the spectrum of the transform domain, d) the reconstructed multiples, and e) the primaries obtained by subtraction. Compare this with Figure 3.1. Clearly, events can be separated very well. Figure 3.4 shows the interpolation results. There is still some decreasing in amplitude in the interpolated traces, but less important than the results shown in Figure 3.2. Chapter 3. Time invariant RT 80 In terms of computing time this approach is fast, but still requires about 5-10 times more operations than using the Levinson algorithm and represents a large improvement compared to the time required by a direct method. 3.3.2 C G - F F T R T This method was first applied to the Radon transform by Sacchi and Porsani (1999). Schonewille and Duijndam (2001) applied a similar idea to the Radon and 2D Fourier transforms. The algorithm exploits the property that the multiplication of a circulant matrix with a vector is equivalent to circulant convolution of a column of the matrix with the vector (Strang, 1986). The technique is certainly very fast, and essentially consists of a C G approach where the slowest part of the algorithm, the matrix vector multiplication, is carried out as a convolution, i.e, as a multiplication in the Fourier domain. Because the R T spreading function matrix, LHL, is not circulant, Sacchi and Porsani (1999) extended this matrix and the model space. The extended part ofthe model can be discarded after every matrix-vector multiplication. The right hand side of the system of equations, the adjoint, has to be padded with zeros as well, to agree in size with the extended kernel. The system to be solved is (3.7) where R e is an expanded form of the spreading function matrix R = LHL, [Re + D] m me 0 Chapter 3. Time invariant RT 81 Figure 3.5: (a) Synthetic gather, (b) Sparse RT obtained with the C G - F F T method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a). Chapter 3. Time invariant RT 82 Figure 3.6: (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse RT obtained with the C G - F F T method, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a). Chapter 3. Time invariant RT 83 R(0) R(n-1) R(n-2) . R*(n-2) R*(n-l) R*(0) R(l) R(0) R(n-1) . R*{n-1) R*(0) R*(l) R{0) . R*(0) R*(0) (3.8) R(n-1) R(n-2) R*{n-2) R*(n-1) In this matrix the R*(i) is the complex conjugate of R(i), added to the original matrix to make it circulant. For every internal iteration in the C G method, a matrix vector multiplication where R e m e is calculated in three steps, 1 - F F T of m e 2 - Multiplication in the Fourier domain of the first column of Hext with m e 3- inverse F F T of the result. The part D e m e is very fast because D e is a diagonal matrix containing the com-ponents of the model covariance matrix. Finally the two vectors are added together, completing the multiplication. The extra elements computed for m are not required and can be discarded. The number of operations required for this method is of order 0(2kextkint2Mlog(2M)), where kext is the number of external iterations, kint is the num-ber of internal iterations and M is the number of traces in the Radon domain. R e + D]m e (3.9) is required. This is carried out as R e m e -I- D m e (3.10) Chapter 3. Time invariant RT 84 Figure 3.5 shows the equivalent of Figures 3.1 and 3.3 for this method. Figure 3.6 shows the equivalent of Figures 3.2 and 3.4. In practical terms, the times required for this method are comparable to Levinson for small and medium sized data. For larger data sets it is even faster. Compared with the previous method ( W T C G L S ) , it has two disadvantages. It requires more memory, as the spreading function matrix is folded and padded with zeros to apply Fast Fourier methods. This disadvantage is negligible in modern computers because of the low price of random access memory ( R A M ) , and the use of F F T algorithms that operate on vectors with a number of elements factorizable by prime numbers rather than powers of 2. The second disadvantage is that the hyperparameter acts on the system by addition, therefore the scale factor is important. This requires a more careful estimation of the hyperparameter than the previous W T C G L S method. However, by applying the transformation to standard form, the matrix vector multiplication in the C G algorithm becomes a multiplication of the extended model vector times the following operator A=(XI + W - T L r W j W (3.11) where A is set to zero for a faster computation and regularization by the number of iterations. Therefore, the action of the operator A can be computed efficiently in three steps • A t a given iteration, multiply the residual vector by the diagonal model weight matrix W" 1 . • Multiply the resulting vector by the matrix L T W j W d L using F F T s , which requires first to extend this matrix to make it circulant (Sacchi and Porsani, 1999). • Multiply the resulting vector from the previous step by the diagonal model weight Chapter 3. Time invariant RT 85 matrix W M . 3.3.3 De-aliased R T Another variant of sparse R T was developed by Herrmann et. al. (1999). This implemen-tation was created with two purposes, to attenuate alias artifacts and to avoid external iterations. The method takes the model information not from a previous iteration but from previous frequencies. To understand why this simple approach eliminates aliased artifacts we need to look at the aliasing condition for the parabolic R T , A h = v j — f — l — — ( 3 - 1 2 ) Jmax \ Qmax Qmin ) where q is the curvature parameter, and H is the total aperture in the offset dimension (aliasing will be further treated in a later chapter). We can see that the frequency appears in the denominator, therefore low frequencies are less demanding in terms of offset sampling and they are in general not aliased. If the model weights are transported from low frequencies to high frequencies it is less likely that they will be affected by aliasing. The assumption behind this method is that the spectrum of the wavelet changes slowly, which is a strong assumption for real data. This method can be implemented by using any of the previously discussed methods, but it works better with the system of equations derived from an underdetermined problem, m = L f f ( L C - 1 L H + A C d ) - 1 d . (3.13) The reason for this is because the model hyperparameter is embedded inside the matrix multiplication and a wrong scale due to variations in the spectrum is less harmful. M y implementation of this method calculates the first 10-15 frequencies of the R T using IRLS. For the following frequencies, the model weights are computed as a function of an average of the model space inside a window that extends from the last computed Chapter 3. Time invariant RT 86 frequency toward low frequencies. Therefore external iterations are required only for the first frequencies. The degree of sparseness can be controlled by the number of external iterations over the first frequencies, by the number of frequencies computed using IRLS and/or by the length of the window used for the model weights. Figure 3.7 shows the multiple attenuation example by applying this method. 3.3.4 H y b r i d t i m e - f r e q u e n c y d o m a i n R T The sparse R T calculated in the frequency domain has a serious problem because the model weights, also computed in the frequency domain, are coupled for all times when in fact they should not be. In the data and Radon spaces, the same weights are imposed for all events in a particular trace. One event at one particular time in the Radon space favors events in that trace for all times. As a result, artifacts are often generated for traces with high energy events. These artifacts can be seen in the sparse R T , in areas where the smooth R T shows no energy. Cary (1998) recognized this weakness and proposed to compute time invariant RTs in the time domain but applying frequency domain operators. As was mentioned be-fore, frequency domain operators have an important advantage over their time domain counterparts; if the computation is done properly inside the bandwidth of the signal, the waveforms are well preserved. Thus, this approach combines the best of both methods, flexibility in the model and data weights and time sparseness and good waveform preser-vation. Because it imposes sparseness in the Radon coordinate and in time, the R T is very clean, with few artifacts. The implementation of this method has some differences with respect to the previously explained methods. The Radon frequency domain operator is applied in cascade with a forward and inverse F F T , Chapter 3. Time invariant RT 87 (d) (e) Figure 3.7: (a) Synthetic gather, (b) Sparse RT obtained with the de-aliased RT method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a). Chapter 3. Time invariant RT 88 (a) q ( s / m A 2 ) 0 2 (b) ( C ) - 1 0 0 0 0 (d) Figure 3.8: (a) Synthetic gather, (b) Sparse RT obtained with the time-frequency domain C G method, (c) Spectrum of (b). (d) Recovered multiples after muting in (b). (e) Primaries obtained by subtracting (d) from (a). Chapter 3. Time invariant RT 89 Figure 3.9: (a) Synthetic data, (b) Synthetic gather with gaps, (c) Sparse RT obtained with the time-frequency domain C G method, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a). Chapter 3. Time invariant RT 90 (3.14) In this equation, d and m are vectors in the time domain, T and T~l are forward and inverse Fourier transforms. The same C G method can be used, but every time the forward or adjoint operators are required in the algorithm, the three operators are applied in sequence. The sparseness in time is sometimes an unrealistic constraint. To attenuate its effect another operator is applied, convolution with a wavelet in the forward operator, and correlation with the same wavelet in the adjoint. Hence, in every iteration of C G L S it is necessary to apply four operators for the forward sequence (convolution, F F T , Radon transform, and inverse F F T ) and for the adjoint ( F F T , Radon transform, inverse F F T and correlation). The price to pay is, of course, more computing time. However, application of all these operators is very fast, and the main burden is that the convergence of C G is slower due to the large number of variables in the time domain. It also requires some extra memory space for storing the three dimensional Radon operator (L(oj,h,q), UJ is frequency, h is offset, and q is the Radon parameter). Figure 3.8 shows the results for the same data set as before, and Figure 3.9 present the interpolation results. The results of this method are perhaps the best of all, but also the slowest to compute. Again it is possible to speed up this approach by using the C G - F F T method explained before. Similar as before, the matrix vector multiplication inside C G is performed by the action of the operator where we have not included W j W d to keep simplicity. We can simplify equation (3.15) A={XI + W - r ^ L T ^ - T ^ - 1 L ^ W - 1 ) , (3.15) Chapter 3. Time invariant RT 91 because for the Fourier transform with sampling Af — (NAt)' (the commonly used F F T ) T~x = TT and T~T = T. (3.16) Hence, A — (AI + W " T J P - 1 L T L j r w - 1 ) . (3.17) Thus, the C G matrix vector multiplication can be applied as • At a given iteration, multiply the residual vector by the diagonal model weight matrix W " 1 . • Take the F F T of the result. • Multiply the resulting vector by the matrix L T W j W d L using F F T s , which requires first to extend this matrix to make it circulant. • Take the inverse F F T of the result. • Multiply the resulting vector from the previous step by the diagonal model weight matrix W " 1 . 3.3.5 Other methods Irregular R a d o n space Schonewille and Duijndam (2001) considered the similarity between the parabolic R T and the 2D F F T , and obtained many useful conclusions about sampling in the transform domain. I found this work highly illuminating, so let me explain their basic results. First of all, they noted that the standard method of using a constant Aq for the R T , based on the maximum frequency of interest, produces a sampling interval which is too small at low frequencies. This is a problem because sampling in the transform Chapter 3. Time invariant RT 92 Figure 3.10: (a) Synthetic data, (b) Synthetic gather with gaps, (c) Irregular RT obtained with a frequency dependent q sampling interval, (d) Spectrum of (c). (e) interpolated data from (b). (f) Residuals obtained by subtracting (e) from (a). Chapter 3. Time invariant RT 93 space implies the existence of some periodicity in the offset domain. If the aperture in the offset dimension is less than this periodicity, then the sampling in the transform domain is equivalent to the introduction of an artificial gap that requires regularization. As the sampling at low frequencies is smaller than the required in terms of the aperture of the data, the introduced gap is larger, and the regularization has to be increased. As a solution, they proposed to use a frequency dependent Radon sampling. Because the aliasing condition for the Radon parameter is inversely proportional to frequency, it is possible to use a Radon sampling Aq proportional to us'1. Wi th this sampling, the parabolic R T is equivalent to a 2D F F T after quadratic stretching of the spatial axis. This allows us to use the same Radon operator for all frequencies. If the R T has to be transformed to time, for example for multiple removal, the transform space must be interpolated to obtain constant sampling. To fill gaps, on the contrary, it is not necessary to obtain the R T in the time domain, and therefore the non constant Aq sampling is not a problem. Two advantages result from using a frequency dependent Aq instead of a constant q interval. First, regularization can be diminished, or even avoided, which results in an improvement in resolution and interpolation inside gaps. Secondly, the operator does not need to be recomputed for every frequency, allowing a series of improvements in efficiency. For example, the actual inverse of the operator can be computed and the solution for every frequency becomes a simple matrix vector multiplication. Figure 3.10 shows the interpolation example using this method. The R T does not show localization in the t — q (time-curvature) domain (Figure 3.10c), because q is frequency dependent along time. Sparseness can still be enforced in the f — q (frequency-curvature) domain (Figure 3.10d). Chapter 3. Time invariant RT 94 Interior point methods This method was explained previously in the chapter about sparse inversion. Even though for small size problems it produces excellent results, for the Radon transform problem it demands more computation time. This is because the linear problem inside the external iterations needs to be solved quite precisely, and therefore an early stopping of the C G algorithm leads to poor convergence. This demands more computing time than for the IRLS method, where the inner linear problem is truncated at the very beginning of the iterations, when only a few elements of the Krilov space contribute to the solution. N L C G w i t h linear search The Nonlinear C G method allows for the computation of the sparse R T by minimizing the £ 2 — 1^ cost function using only the cost function and its gradient. In my experience, however, this method does not produce very good results for this problem, but I mention it here because it might be useful for adding more constraints and experimenting with different norms. Because the Hessian matrix is not required, the complexity of the cost function can be increased, as long as the gradient can be computed. 3 . 4 R e a l data examples for marine and land data The previous examples on synthetic gathers have shown notable improvements in multiple attenuation and interpolation using a sparse R T . However, real data tests are much more demanding because when complicated events map to the transform domain, they produce complex signatures that spread over the R T space. Often, these signatures do not map completely inside the transform space and serious artifacts arise at the boundaries of the R T . Because the increase in computation time for the sparse R T is not negligible compared Chapter 3. Time invariant RT 95 Figure 3.11: (a) Marine CDP. (b) Primaries obtained by subtraction, (c) Multiples obtained from the inverse RT of (e). (d) Non-Sparse RT. (e) RT after mute of primaries. Bottom - Semblance analysis for: (f) data, (g) primaries, and (h) multiples Chapter 3. Time invariant ET 96 Figure 3.12: (a) Marine CDP. (b) Primaries obtained by subtraction, (c) Multiples obtained from inverse RT of (e). (d) Sparse RT - W T C G L S . (e) RT after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples Chapter 3. Time invariant RT 9 7 o f f s e t ( l e e t ) - 1 . 5 - 1 . 0 - 0 . 5 o t f s e t ( f e e t ) -1.0 -0.5 X10* 0 o f f s e t ( l e e t ) -1.0 -0.5 X 1 0 - 4 0 (b) Primaries (d) RT sparse: Time-Frequency domain V e l o c i t y { m / s e c ) 4000 5000 6000 7000 8000 0 i MUTE (c) Multiples q ( s / f e e t * 2 ) V e l o c i t y ( m / s e c ) 5000 6000 7000 (e) RT after mute V e l o c i t y ( m / s e c ) 4000 5000 6000 7000 8000 (f) Semblance: Data (g) Semblace: Primaries (h) Semblance: Multiples Figure 3.13: (a) Marine CDP. (b) Primaries obtained by subtraction, (c) Multiples obtained from inverse RT of (e). (d) Sparse RT-Time frequency domain, (e) RT after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples Chapter 3. Time invariant RT 98 to non-sparse RTs, it is important to be aware under what circunstances the method will improve the results considerably. Even though the method is always better for synthetic and well behaved data, it often does not produce considerable improvement for noisy data. The improvements of the results depend on, among other things, how noise and aliasing affect the algorithms and the validity of the assumptions about sparseness. For noisy real data with gaps and non hyperbolic moveouts, the ideal R T is not necessarily sparse and attempts to obtain sparseness in the transform could result in increasing noise. Another problem we can expect is that often, during the inversion, we have to stop the algorithms far from the optimal solution for the proposed cost function. In these cases the results resemble the least squares solutions and the improvements, if any, are difficult to see. Figure 3.11 represents the cycle for multiple attenuation in a C M P using the R T . Figure 3.11a shows a marine C M P gather from the Gulf of Mexico (Mississippi Canyon data, courtesy of Western Digital). The gather is transformed by using the R T after N M O correction, such that the primaries will have been flattened whereas the lower velocity multiples will show a larger residual moveout. Therefore they will map towards the right side of the zero moveout central trace. The N M O / R T procedure puts primaries and multiples in two different areas of the transform domain. By applying a mute to the primaries (see Figure 3.l ie), and applying an inverse R T to this space, the multiples are recovered (shown in Figure 3.11c) and subtracted from the data to obtain the primaries (Figure 3.11b). A final semblance analysis or velocity spectrum is shown at the bottom of Figure 3.11 to judge the quality of the procedure. In general the efficiency of the method depends on: 1 - The differences between the stacking velocities of primaries and multiples. 2 - The resolution of the Radon transform used. If the N M O correction is not accurate, multiples will still map towards the right hand Chapter 3. Time invariant RT 9 9 side of the primaries and separation is possible. However, a N M O correction based on high velocities may give problems, because multiples could map outside the Radon space, to the right hand side of it, and primaries might be confused with multiples. O n the other hand, choosing low velocity may flatten multiples, and map primaries outside the Radon space, to the left hand side of it, causing multiples to be confused with primaries. The necessity of some velocity trend estimate is common to many seismic processing tools. This example with the non-sparse R T shows that even when primaries and multiples may have very different moveout, the lack of localization of the R T can make it difficult to mute the signal without touching the multiples. In this situation, part of the multiples will remain in the final filtered C M P gather. Figure 3.12 shows a similar example as before but using the W T C G L S method for the sparse P R T . Even though the result is not much better, we can see that there is a better localization of events between 4.5 and 7 s. Figure 3.13 shows the result using the time-frequency domain. Now it is possible to see a differentiation between the primary and multiple trend between 4.5 and 7 s. Figure 3.14 presents the Mississippi Canyon data after N M O correction and stack. Between C M P s 1500 and 2400, the apparent horizons at 3.5, 4 and 4.5 s are multiples. Figure 3.15 presents the same data after multiple removal using the P R T with W T C G L S , N M O correction and stack. Most of the multiples have been attenuated by subtraction. The processing of the total line took about 20s per C M P (in a Pentium III, 700 MHz) but it can be reduced by a factor of three by using model weights from contiguous C M P s . Figure 3.16a presents the same C M P as before, but after some artificial gaps have been created. Figure 3.16b shows its non-sparse R T , Figure 3.16c displays the interpolated data, and Figure 3.16d reveals the residuals between original data without gaps and data generated from the Radon panel. Figure 3.17 shows the same process for the sparse R T . By looking at the residuals it seems that the interpolation is very poor in both cases. The Figure 3.14: (a) The Mississipi Canyon data after N M O correction and stack. Figure 3.15: (a) The Mississipi Canyon data after multiple removal with P R T ( W T C G L S ) , N M O correction and stack. Figure 3.16: (a) Marine C M P . (b) Non sparse R T obtained with Levinson recursion, (c) Interpolated data obtained from (b). (d) Residuals Figure 3.17: (a) Marine CMP. (b) Sparse RT obtained with the W T C G L S method, (c) Interpolated data obtained from (b). (d) Residuals. Chapter 3. Time invariant RT 104 residuals show high energy mainly because there are many local features on the original data that cannot be reproduced in detail by the R T interpolator. The cause of this problem is due to the fact that either these details do not map to the R T , for example if they are not recorded outside the gaps, or these fine data features can be reconstructed only by using very small singular values of the transformation. Still, by looking at the predicted data in the gaps we can see that the interpolation is partially satisfactory for all offsets in the sparse R T example, and for near offsets in the non sparse R T example. Better results can be obtained by working with hyperbolic RTs in the time domain, as will be demonstrated later. Another example with lower fold data is shown in Figure 3.18. The data gather in Figure 3.18a corresponds to a C M P gather from a 3D land data set. As usually happens with C M P s from 3D data, the sampling along the offset axis is very irregular. Figure 3.18b displays the standard non-sparse Radon transform of this gather, Figure 3.18c is the sparse R T obtained with the W T C G L S method. Even though it is difficult to speak of sparseness in such a complicated R T , we can see that the lateral resolution of the events has increased considerably. The problem with this example is the introduction of artifacts in the transform domain. The increase in lateral resolution brings some distortion in the time direction and a loss of time resolution. A time domain method eliminates this problem by enforcing sparseness in the time axis as well. Figure 3.18d shows the sparse R T computed using the mixed time-frequency domain method ( T F D ). Interpolation of this same C M P can be achieved with any of the previously discussed methods with good results. The symmetry of the P R T allows the algorithm to fill gaps at one side using information from the other side of the symmetric axis (Figure 3.19). However, by eliminating some of the traces of the original C M P , we can see how the non sparse R T loses its capability to fill the gaps, whereas the sparse R T still fills them. Figure 3.20 displays the same C D P as before with some traces removed. The interpolation by Chapter 3. Time invariant RT 105 (a) Data (b) Non sparse (c) WTCGLS (d) TFD Figure 3.18: (a) C M P from land data, (b) Non sparse RT. (c) Sparse RT obtained with the W T C G L S method, (d) Sparse RT obtained with the time-frequency domain method. Chapter 3. Time invariant RT 106 offset(m) offset(m) (c) (d) Figure 3.19: (a) C M P from land data, (b-c-d) Interpolation by using b) non sparse RT. c) Sparse RT obtained by the W T C G L S method, d) Sparse RT obtained by the time-frequency domain method. Chapter 3. Time invariant RT 107 Figure 3.20: (a) C M P from land data after eliminating some traces, (b-c-d) Interpolation by using b) non sparse RT. c) Sparse RT obtained by the W T C G L S method, d) Sparse RT obtained by the time-frequency domain method. Chapter 3. Time invariant RT 108 Algorithm synthetic ex 1 synthetic ex 2 marine data ex. 1 L E V I N S O N 1 1 3 W T C G L S 3 2 16 C G F F T 3 2 5 T I M E - F R E Q U E N C Y D O M A I N 27 16 47 Table 3.1: Computing times in seconds for the explained methods in the different exam-ples. using the non sparse R T (Figure 3.20b) has failed in the large gaps, while in contrast, the interpolation by sparse R T (Figures 3.20c and d) has been quite successful. Hence, for this example interpolation has been achieved for different reasons in each method. The sparse R T interpolates the gaps by minimizing the amplitude variations in the transition from traces to gaps and also by preserving symmetry in the events. The smooth R T does not penalize very strongly amplitude variations and therefore depends more on the symmetry of the gathers. Therefore, the smooth R T fails to interpolate when the information is missing at both sides of the gather. In my experience, when the gathers have gaps in one side but information in the other, interpolation can be achieved almost equally by the sparse or the smooth R T ; the contrary applies for missing information at both sides (or one side gathers). The computation times for the different explained methods appear in table (3.4). These tests were performed on a Pentium III, 700 M H z computer. 3.5 Conclusions From the examples considered in this chapter, we can conclude that the performance of the sparse R T in separating primaries from multiples and interpolating gaps, is superior to the performance of the standard R T . Nonetheless, its application for real data is non trivial, and well designed algorithms are required. Not only must the algorithms be fast Chapter 3. Time invariant RT 109 and robust, but they must be easy to implement. The choice of the hyperparameters could be a serious problem as different results can be obtained with different values. In this sense, the combined W T C G L S - G C V algorithm was demonstrated to be very practical. The G C V function stops the inner iterations as soon the algorithm starts to fit noise. The number of external iterations is less of a problem, because the IRLS algorithm tends to converge to a stable solution. The C G - F F T method is the fastest, and it should be applied when the size ofthe data is very large. It is, however, more sensitive to the choice of the hyperparameters. Finally, the hybrid time-frequency domain method overcomes the most difficult problem faced by previous methods, the lack of time localization of the model weights in the frequency domain. Its drawback is a large computation time. Other problems that are important to consider are 1 - Real data can often produce non sparse RTs because of non hyperbolic paths that cannot be N M O corrected appropriately for the P R T . This is very common for data collected in areas of large velocity contrast, for example, data associated with salt domes or areas with strong anisotropic properties. A V O , noise and gaps are other reasons for non sparseness in the transform domain. 2 - Aliasing artifacts increase with iterations and often produce accumulation of energy at the borders of the Radon space. The first problem could be mitigated by the addition of an extra term in the basis functions, similar to the one used for N M O . The second problem often can be diminished by a careful choice of the limits for the Radon space, and a good stopping criteria for both external and internal iterations. These topics are the subject of a future research. Chapter 4 Time variant Radon transforms Seismic reflections in C D P gathers are well approximated in general by hyperbolas. Therefore, a sparse hyperbolic Radon transform (HRT) is ideal to map reflections onto a sparse domain. This allows one, as explained before, to extend the aperture of limited data, filter noise and fill gaps. In the same manner, reflections in slant stack sections (gathers after a linear Radon transformation) can be well approximated by ellipses, there-fore, an elliptical R T (ERT) will give similar results to a H R T when applied to slant stack sections. A problem with these transformations is that they have a time variant kernel that results in slow implementations. However, the computation time can be kept low enough for practical applications by defining a model space in terms of an irregularly sampled velocity space to minimize the number of unknowns during the inversion and using sparse matrices . I implement hyperbolic and elliptical time domain RTs , by using inversion via weighted conjugate gradient methods with a sparseness constraint. The hyperbolic R T performs accurate interpolation in C M P gathers, while the elliptical R T attenuates sampling artifacts in slant stack sections obtained from C M P gathers with poor sampling and gaps. 4.1 PRT vs HRT The Radon transform has been shown to be an efficient interpolator (Kabir and Ver-schuur, 1989). A hyperbolic Radon transform applied to a relatively noise free C M P gather with large aperture produces a sparse model. If the C M P gather contains gaps 110 Chapter 4. Time variant Radon transforms 111 and/or noise, improvement can be achieved by forcing sparseness on the Radon model (Thorson and Claerbout, 1985). Faster parabolic sparse Radon transforms can be ap-plied to N M O corrected C M P gathers to extend aperture and fill gaps (Sacchi and Ulrych, 1995). The quality of the interpolation and extension depends, however, on the similar-ity of the seismic events to the basis functions. Thus, a true hyperbolic R T in general produces better results than the parabolic version for data with large moveout. A prob-lem with the H R T is that it is time variant so it is not possible to use fast solvers like Levinson recursion (Rostov, 1990). Similar to other large ill conditioned problems, the H R T is better solved by iterative methods like, for example conjugate gradient and L S Q R (Saad, 1996), because they produce an approximate solution at each iteration that can be stopped at any moment to limit computation time. However, the time at each iteration must still be reduced as much as possible and a very simple way of achieving this is by using traces in the Radon domain with velocity dependence on the two way travel time r. M y implementation consists of a trace located at a selected position in the model space that contains a velocity trend estimated by semblance analysis, and neighboring traces that contain perturbations of the central trace. These neighboring traces can be either at constant intervals or not. Sometimes, variable intervals produce better results. There are advantages and disadvantages in a hyperbolic R T versus a parabolic R T . Looking at the location of the events in the model space there is a similarity between applying parabolic R T to N M O corrected traces, and applying time dependent velocity hyperbolic R T to the original traces, in the sense that primaries fall close to a particular trace, and multiples and coherent noise fall further away from it. However, a better approximation of reflections can be expected with hyperbolic basis functions, in partic-ular for large offsets, common in marine data. On the other hand, the NMO-parabol ic approximation is computationally more efficient. Another interesting aspect of irregu-lar sampling is how it affects aliasing. Aliasing is normally considered to occur as a Chapter 4. Time variant Radon transforms 112 consequence of regular sampling, but an irregularly sampled Radon space allows the re-laxation of conditions for the minimum sampling required to avoid aliasing (Trad and Ulrych, 1999). 4.2 T i m e variant R a d o n transforms 4.2.1 The Hyperbo l i c R a d o n Transform The hyperbolic Radon Transform (HRT) is defined in terms of summation along hyper-bolic paths. If we denote the data in the time-offset domain by d(t, h) and the Radon transformed data (a Radon panel) by m(r, v), the hyperbolic Radon transform is defined as Climax I m{r,v)= I d(t = ^/(r 2 + h2/v2), h) dh, (4.1) where h indicates offset (the source-receiver distance in a common midpoint ( C M P ) gather). The variable r is the two-way zero offset time and v denotes the root-mean-square (RMS) velocity. Equation (4.1) represents a decomposition of a C M P gather in terms of hyperbolas (Thorson and Claerbout, 1985). I use a variant of this representation with v = J(T). Equation (4.1) can be discretized, rearranged and represented in matrix form as follows: d = L m , (4.2) where d and m represent the C M P gather and Radon panel after lexicographic arrange-ment, respectively. The linear forward modeling operator L can be efficiently manipulated using sparse matrix techniques. The latter is examined in detail later. A low resolution Radon panel can be estimated by using the transpose or adjoint operator L T . In this case, the low resolution Radon panel is given by Chapter 4. Time variant Radon transforms 113 m ^ - = L r d . (4.3) It is clear that since L is not an orthogonal operator, L T does not define the inverse operator. Thorson and Claerbout (1985) proposed the computation of the Radon trans-form by inverting equation (4.2) rather than using the low resolution Radon transform L T . For that purpose they proposed a stochastic inversion technique that is able to retrieve a sparse Radon panel. 4.2.2 E l l i p t i c a l R a d o n T r a n s f o r m Inspired by the H R T , I developed a method to improve the continuity of the elliptical events of slant stack sections and filter out some artifacts introduced by sampling prob-lems in the original time - offset gather. The method I implemented is inspired by the H R T in the sense that a summation is performed along elliptical paths in the R T domain, rather than hyperbolic paths in time. A motivation for such a procedure is the character-istic of slant stack sections of having infinite offset aperture in a limited transform space. This is an advantage for aperture extension. What is required is a transformation from the slant stacks to a new domain where the continuity of the slant stacks is improved by forcing sparseness. The first of such new domains we can imagine is the (r — V) domain, similar to the velocity stacks. To apply the conjugate gradient methodology we need to define the operators to map between the slant stack domain d(r,p) and the model space m(r, V) (adjoint operator) and vice versa (forward operator). In analogy with the previous definitions, the elliptical R T is given by a decomposition of the slant stack in ellipses, fPmax m(r0,V)= d(T(p,V,r0),p)dp, (4.4) ^ Pmin Chapter 4. Time variant Radon transforms 114 where r(p, V, To) is obtained as r(p,V,r0) = ^ c o s ( r 9 ) = T 0 / L - p 2 V 2 . (4.5) In these equations V is velocity, r the two way travel time and p = sm(8)/V is the ray parameter. The elliptical R T proceeds exactly as the hyperbolic R T . During the forward operation, again represented by, d = L e m , (4.6) every point of the model m produces energy in the data d along ellipses given by the operator L e for every particular velocity. A low resolution elliptical transform madj could be computed by using the adjoint operator, L T , as madj = L j d , (4.7) i.e., a summation of the slant stack data following elliptical moveouts. However, the same weighted conjugate gradient least squares algorithm used for the hyperbolic R T can be applied, to invert the forward problem. 4.3 I m p l e m e n t a t i o n 4.3.1 C G so lut ion for the l inear p r o b l e m Performing the Radon transformation by least squares inversion in the time domain requires the solution of a large ill conditioned problem. The size of the (sparse1) matrix of coefficients is (nt x nh) x (nr x nV) where the transform maps from a model space m(r, V) in time-velocity (or some function of velocity) to the data space d(t, h) in the time offset or d(r,p) in the slant stack domain. I use a weighted preconditioned conjugate gradient (CG) algorithm for normal equations to solve the system of equations with weights for lThe word sparse refers here to the property of a matrix of having many elements equal to zero. Chapter 4. Time variant Radon transforms 115 model and data. The preconditioning allows a sparse solution and the weighting allows a softening of the data constraint for bad data. Some of the advantages of this type of method have been already discussed in Chapter 2, but particularly important for the time domain implementations is that the only required operations are forward and adjoint operators. In computing the Radon transform, these operators can be applied as sparse matrix vector multiplications and can be performed very efficiently using parallel architecture, (Hansen, 1998). The algorithm I use is a left preconditioned version of conjugate Gradient for the Normal Equations denoted C G N R (N for "Normal" and R for "Residual") (Saad, 1996) with data weights. The symmetry for the matrix of coefficients M_1ArA is obtained by replacing the Euclidean inner product for a preconditioned inner product. The precon-ditioner is the inverse of the diagonal model covariance matrix, equations (2.8) or (2.10). The pseudo-code is shown in Figure (4.1), where x is the solution vector, b are the data, r the residuals, g is the gradient vector, z, is the preconditioned gradient and p is the search direction vector. 4.3.2 Fast sparse implementation for the operators As explained in Claerbout (1992) many inversion and seismic processing schemes can be implemented by C G techniques, just by finding the forward and adjoint operators. The operator that produces the forward and adjoint transforms can be implemented either as a operator, i.e. as a procedure or algorithm, or as a sparse matrix. The first approach saves memory, the second one saves time because it allows fast matrix vector multiplications. The operator implementation gives more freedom in the design of the algorithm, for example applying sine interpolation for estimating the actual value of a variable at a time between samples. Also, memory for the matrix is not required as every operation is carried out directly on the data and model vectors. The drawback of this approach is Chapter 4. Time variant Radon transforms 116 Algor i thm 1 Left Preconditioned C G N R with data weights r 0 = W d ( b - Lx 0 ) , go = L r W j r 0 , z 0 = M _ 1 g 0 , p 0 = z 0 while (r(k) > e) A N D (k < k m a x ) do w fc = WdLp*. a k ~~ llwjfcHi x f c + i = x fc + akpk r k + i = rk- akwk g f c + i = LTWdrk+1 zk+1 = M~lgk+1 a _ ( z f c + i i g f c + i ) - (z*,gfc) Pfc+1 = zk+l + PkPk k = k + 1 end while Figure 4.1: Algorithm P C G N R that it requires much more computing time because of repeated calculations from input to output. The sparse matrix implementation, on the other hand, is much faster than the operator implementation but requires extra random access memory ( R A M ) to store the sparse matrix. A few years ago, there was no doubt that saving memory was very important. Nowa-days memory is cheaper than time, and the actual computation of the matrix is more advantageous. Most sparse applications require the saving ofthe indices for every nonzero element and its value. In the R T kernel the sparse matrix has all nonzero values equal to one and only the indices have to be saved in the R A M , not the matrix coefficients themselves. If linear interpolation is applied, the indexes are multiplied by ten before rounding to integers to preserve the first decimal digit. For a typical data set (for exam-ple, the marine C D P presented later with 1100 time samples and 132 traces), about 100 M b of R A M are enough. In the pseudo-code in Figure 4.3, index is a matrix with two columns of nsparse length Chapter 4. Time variant Radon transforms 117 fdr(ti=(y,it<nspars§it++) d[inde40] [it] ]+=m[inde^\] [it] S i m p l e c o d e : A d j o i n t •fpr(ii=p.if<nspars$it++) m[inde^l] [/(]]+=d[inde>{0] [it];; Figure 4.2: Sparse matrix operations. A two column matrix connects the data and model spaces. each. The first column contains the index for the output, the second one the index for the input. Thus, every element of index corresponds to the indices of the operator, and the value of the element located there is one. Matrix vector multiplication is performed as shown in the pseudo-code in Figure (4.3) (a few extra steps are required for linear interpolation). For this implementation I have applied linear interpolation rather than sine interpolation in order to decrease computation time. Sometimes I incorporate a convolution with a band pass wavelet into the forward operator and a correlation into the adjoint operator. Both the hyperbolic and elliptical operators are designed in the same manner, as explained above. However, a problem in the elliptical interpolator is that there are elements of the data space that cannot be reached by the range limited operator. A vector with the size of the data space is simultaneously set with zeros and ones during the computation of the sparse operator, assigning zeros for the unreachable part of the data space. The original data values are kept for those locations when the inverse transform Chapter 4. Time variant Radon transforms 118 Algorithm 1 Sparse operator { d vector length ny with all data traces} { m vector length nx with all model traces} if Adjoint = = F A L S E then d = 0 for j = 0 to nsparse do d[mdea;[0][j]] = <i[mc/ex[0][j]] + m[mdea;[l][j]] end for else m = 0 for j = 0 to nsparse do m[mdex[l][j]] = m[mdea;[l][j]] + d[mdex[0][j]] end for end if Figure 4.3: Sparse operator is performed. 4.3.3 Velocity grid To create the non-uniform grid of velocities I define a vector with velocity increments, A V , whose elements increase by a small user defined quantity, e, as we move far from the central velocity trend. Also a minimum velocity increment, AVmin, is used to define the minimum separation between velocity traces. Starting from the velocity law estimated by semblance analysis at the trace n/2, where n is the number of traces in the Radon space, I compute AV as follows AViq+x = AViq(l + e) + AVmin for iq = n/2 to n (4.8) A l V i = AViq(l + e) - AVmin for iq = n/2 to 1 The velocity at every grid point is given by the original velocity law plus the increment. The square of the velocity is required so the grid is Gi,,* = Vllt + AVl + 2ViqMAViq (4.9) Chapter 4. Time variant Radon transforms 119 where Giq^t is the grid stored in the R A M . The velocity increment vector does not need to be constant with time. I actually have obtained better results by multiplying the increments by (l+a^ 4 ) , where tu is time in seconds, and a is a constant. This modification allows the velocity span to increase with time, which is sometimes convenient as the aperture decreases with depth and so does the sparseness. This is a completely empirical approach and is designed to sample more densely the energy of events with the primary velocity trend. 4.4 Examples I make use of a set of hyperbolic events with a central gap and far offset truncation to show the capability of both methods, hyperbolic and elliptical. I used either a Cauchy norm or the li, depending on the result. In general both norms lead to sparse models, and I do not have any special reason to choose one over the other. A l l the computing times are evaluated on a Pentium III processor, 700 M H z , with 256 M b of R A M . The number of traces and iterations have been reduced as much as possible, while still producing acceptable results, to reduce the computation time. Figure 4.4a shows a simulated incomplete C D P gather and Figure 4.4b the irregular r dependent hyperbolic Radon space. Notice that all events fall close to the central trace that contains the background velocity as a function of r as estimated from semblance analysis. In this way the size of the model space is considerably reduced. Figure 4.4c shows the interpolated section obtained by applying the forward Radon operator to the panel in Figure 4.4b. Figure 4.4d shows the residual between the simulated gather without truncation or gaps and the interpolated panel. For this example (nt = 512, nx = 62, J I T = 512, nV = 40) 10 internal iterations and 3 external iterations were carried out for the Radon transform in 2.5 seconds. Figure 4.5 shows the same test using the standard Chapter 4. Time variant Radon transforms 120 parabolic Radon transform with Levinson recursion, which is the fastest method available for this size of data. Note that the horizontal axes for the Radon spaces are different. The R T using 65 traces for the Radon space required 1 second in my implementation. For a small data set like this the interpolation has been considerably improved at the expense of a small increase in computation time. This low cost of the time domain R T unfortunately increases as we work with larger data sets. However, this simple example shows the most important difference between the two methods. P R T with least squares interpolates well for near zero offsets, but not very well for large offsets. A high resolution P R T (see Chapter 3) interpolates all offsets, but we will see later that even this method does not interpolate very well when events with large and different moveout are present for the same two way travel time. Figure 4.6a shows the same synthetic hyperbolic data used in the previous example, with the same near offset gap and far offset truncation (Figure 4.6b). After mapping to the linear Radon domain the complete data set, the ellipses in the slant stack domain appear very regular and continuous (Figure 4.6d). Note, however, some artifacts because of aliasing in the original complete data set, visible in its f-k spectrum (Figure 4.6c). These artifacts are going to be removed as well by the elliptical R T . The Linear Radon transform from the incomplete data set shows several more artifacts, due to truncation and gaps (Figure 4.6e). Figure 4.6f shows the difference between the slant stack computed with the complete gather and the computed with the partial data set. The slant stack obtained from partial data is transformed to the model space (Figure 4.6g) using an elliptical R T , and sparseness in the velocity space is enforced to improve the continuity of the events in the slant stack domain. Figure 4.6h shows the slant stack obtained from Figure 4.6g using the elliptical R T operator. Figure 4.6i is again the difference between slant stack from complete data and the slant stack from the incomplete data + elliptical Radon (inverse and forward). We see that the ellipses are more continuous and, more Chapter 4. Time variant Radon transforms 121 importantly, several artifacts have almost disappeared. The reason for this improvement resides in the regularization properties of C G that first fit the data that resemble the basis functions, and only in later iterations fit the small details like artifacts. Furthermore, many of these artifacts belong to the nullspace of the elliptical R T . A n unfavorable result is that the extrapolation at far offsets is limited. To try to avoid this problem, an elliptical taper is applied to prevent a undesired fit of truncation artifacts at large p values. I have not been completely successful with this approach, and in general such a taper is difficult to design because it is not clear where the truncation artifacts start along the elliptical events. Also a taper has been applied for small p values where the gap produces low amplitude values for the ellipses. Because the elliptical events on the L R T for real data become complicated after the point were evanescent waves appear, a data weighting function that limits the extension of the elliptical basis functions has been applied as well. The next example is a marine C M P gather (Figure 4.7a) from the Gulf of Mexico (the same gather used in Chapter 3). The original gather has full coverage from -4900 to -18 m and, because of this, is an excellent data set to test interpolation methods. I remove traces in three ranges of offsets: near offsets, gap between -460 to 0 m; middle offsets, gap between -2600 to -2140 m; and large offset, gap between -4740 to -4276 m. Comparison between the original and interpolated traces will give us information about the quality of the waveforms and times of the seismic events in the created traces. Figure 4.7b shows the modified hyperbolic Radon transform. It is important to note that the location of the different events in the Radon space depends on the chosen parameters. For example, the central velocity trend can correspond to the velocity of multiples, primaries, or a constant value, depending on which events we are interested in. Hence, it could sometimes be better to sample multiples more densely because of their large moveout. Also, the velocity sampling changes at a desired ratio when moving away Chapter 4. Time variant Radon transforms 122 from the central velocity trend. This velocity trend can itself be located at different positions on the model space. Tests suggest that better interpolation results can be obtained for some reflections by changing these parameters. The extra degrees of freedom in the selection of the model space will probably prove to be useful for other applications than interpolation. Figure 4.7c shows the original gather with the gaps interpolated with the synthetic traces. Figure 4.7d shows residuals between the complete original data set, and the complete synthetic gather. A few comments are in order. First, we see that some of the high energy reflections on the shallow part of the existing traces have not been preserved completely for large offsets and moveouts. This indicates poor reconstruction of the waveforms or times. However, it is not itself a problem because these traces are not going to be used for interpolation, as they are already in the data set. The large moveouts of these events make them difficult to preserve. Their high energy makes it easy to note any distortion relative to the others (for the parabolic approach these events would suffer from N M O stretch). Let us now focus on the interpolated traces in the gaps. The residuals for these traces contain coherent energy together with noise. Nevertheless, we are creating information in certain areas of the data space using information from other areas of it. The data contain irregularities and amplitude variation with offset (AVO) , features that are local in nature, and not always possible to create from the model. Therefore it is important to judge the results relative to what can be obtained with other methods, like the P R T . Comparison with other interpolation techniques that allows irregular original sampling is beyond the scope of this thesis, but also important. For example it would be worthwhile to compare with interpolation in the f-k domain, or time domain dip moveout. I will mention only that in some tests of interpolation in the f-k domain, I found that the interpolation of large gaps in C M P s was worse than with the H R T . Good results can be obtained, however, by interpolating in N M O corrected data or constant offset sections. Chapter 4. Time variant Radon transforms 123 Figure 4.8 shows a series of details into the original and interpolated data for the different offset ranges. To make a fair comparison between P R T and H R T I have used a sparse R T for both cases. A comparison with standard non sparse P R T is unnecessary because it has been proven already that sparseness improves aperture extension and interpolation (Chapter 3). The first column shows the original data, the second column the original data interpolated with sparse P R T , and the third one interpolated with H R T . Every row represents a different offset range, i.e., first row shows near offsets, second row displays middle offsets, and third one contains large offsets. As expected, sparse P R T performs better at short offsets, simply because the frequency domain operators are able to preserve more accurately the waveforms. However, middle and far offsets are better interpolated by the H R T , because the P R T approximation cannot fulfill at the same time the large moveout of the multiples and primaries. The H R T required 22 seconds to obtain the result in Figure 4.7c with 15 C G iterations and 3 external iterations to update the model weights (about one second per iteration). The same calculation took 20 seconds for the result in Figure 4.7b (sparse P R T ) . A better alternative to external iterations is to use the model weights from previous C M P s when processing along the line. The model weights do not have to be precise, therefore smooth changes from C M P to C M P are not important. I use the same marine gather with three gaps to show how an elliptical R T can eliminate several artifacts in the L R T that originate from the missing traces. Figure 4.9a shows the complete C M P gather at location 1000. Figure 4.9b is its L R T . Figure 4.9c is the same C M P gather with gaps and is resampled with one every two traces (decimated in offset). Figure 4.9d is the L R T of the incomplete decimated gather. The Radon model has plenty of artifacts as a consequence of gaps and poor sampling in the offset domain. Figure 4.9e is the velocity model obtained from the L R T in Figure 4.9d by using the elliptical R T . Figure 4.9f is the linear R T of the incomplete gather obtained from the Chapter 4. Time variant Radon transforms 124 velocity model in Figure 4.9e by using the inverse elliptical R T . The H R T presents some advantages over the P R T for multiple attenuation in data with large moveout. In these cases, multiples can be better separated in hyperbolic than in parabolic Radon domains. Figure 4.10 shows the same marine C M P gather. The gather is transformed by using the H R T . The primaries will map close to the central trace of the H R T , which contains the primary velocity trend. Lower velocity multiples will map to the right hand side of this trend. The H R T procedure puts primaries and multiples in two different areas of the transform domain. By applying a mute to the primaries (see Figure 4.10e), and applying an inverse R T to this space, the multiples are recovered (shown in Figure 4.10c) and subtracted from the data to obtain the primaries (Figure 4.10b). A final semblance analysis or velocity spectrum is shown at the bottom of Figure 4.10 to judge the quality of the procedure. Figure 4.11 shows the result of H R T (multiple attenuation) + N M O + stack on the whole line. This Figure has to be compared to Figures 3.14 and 3.15 (NMO+stack and P R T + N M O + S t a c k respectively) to judge the quality of the results. The multiples have been attenuated slightly better than with the P R T . The processing time for P R T and H R T are very similar in this data (about 20s per C M P in a Pentium III, 700 M H z ) , mainly because the multiple prediction by R T does not require to compute the complete Radon space for the time domain algorithms. 4.5 Discussion The examples illustrate that a sparse hyperbolic or elliptical Radon transform in the time domain can be effectively applied to perform interpolation and extrapolation by using • A weighted preconditioned conjugate gradient least squares method. • Sparse matrix upload in R A M . Chapter 4. Time variant Radon transforms 125 • A n irregularly sampled Radon space. The same procedure can be applied for any curve described by an analytical expression and one unknown parameter for every element of the model space (i.e. every point of the model space represents a curve that can be expressed as a function of time, offset and one extra parameter from the model space). Curves with more parameters would complicate the transform but could be useful, for example, to model A V O effects. A n important point to note for the sparse R T is that the events do not, in reality, have to be perfectly hyperbolic or elliptical. Naturally, the closer to hyperbolas or ellipses, the sparser is the representation. There is, however, a problem in the slant stack interpolation that is worthy of con-sideration. When the goal is to interpolate gaps and extend aperture, predicted traces can replace zero or empty traces. The elliptical slant stack domain does not have gaps but only lack of continuity due to gaps in the time offset data and artifacts due to poor sampling in the original data. These artifacts are difficult to separate from the data and a conservative subtraction process, similar to the prediction and subtraction of multiples, is not possible. So far I have replaced original traces with synthetics from the elliptical operator. There is no easy criterion to decide whether a particular filtered trace is better than the original trace. A worthwhile approach, even though more computations are required, might be to go back to time offset, interpolate with the hyperbolic R T and then transform back to the slant stack domain. However I still believe that results with a direct elliptical R T could be useful since the elliptical operator can gather or predict energy from, or to, an infinite theoretical aperture, which is not possible in the time-offset domain. This capability is useful if the aperture can be extended and artifacts attenuated while preserving the original data completely. I do not, as yet, know how to implement this. Chapter 4. Time variant Radon transforms 126 4.6 Conclusions For data with large moveout, as is often the case for marine data, the hyperbolic transform is more accurate than its parabolic counterpart but much slower to compute. I have presented here a method of decreasing the computational workload in order to implement hyperbolic interpolation. I make use of a r dependent model space and a weighted conjugate gradient algorithm implemented with sparse matrix vector multiplication. I have also extended the idea of interpolation by means of a sparse R T to the slant stack domain. A n elliptical R T in the time domain obtained by conjugate gradient techniques and a r dependent non-constant sampling model space is used, which is made possible due to the flexibility of the summation method in the operators. Memory requirements can be minimized by working with operators or computing time can be reduced by one order of magnitude by working with sparse matrices. Results with synthetic and real data suggest that interpolation of central gaps, aperture extension and hyperbolic filtering can be simultaneously achieved. As a final comment, notice that the computation times are about 20s (Pentium I I I , 700 Mhz) for C M P gathers of considerable size. These computing times can actually be reduced about three times if the model weights are estimated from the previous C M P and, of course, even more by using more advanced hardware. Based on these facts I believe that a sparse time domain H R T is not impractical and may be used in real data processing providing that it presents some advantages over a faster R T . M y conclusion is that for data sets with gaps at large and middle offsets, the H R T is an interesting possibility in spite of the extra (but reasonably so) time required. Chapter 4. Time variant Radon transforms 127 offset(m) q(s^m2) x10" (c) (d) Figure 4.4: (a) Synthetic C M P gather with a near offset gap and very limited aperture, (b) Sparse hyperbolic RT. Horizontal axis represents slowness at zero time, (c) Gather produced by the Radon model, (d) Residuals between complete synthetic gather (not shown) and recovered after RT. Chapter 4. Time variant Radon transforms 128 Figure 4.5: (a) The same synthetic C M P gather with a near offset gap and very limited aperture, (b) Nonsparse parabolic RT computed in the frequency domain (after NMO) . Horizontal axis represents residual moveout per offset square, (c) Gather produced by the Radon model, (d) Residuals between complete synthetic gather (not shown) and recovered after RT. Chapter 4. Time variant Radon transforms 129 Figure 4.6: (a) Complete simulated gather, (b) Gather a) with missing traces and trun-cation, (c) f-k spectrum from a), (d) Slant stack from the complete gather, (e) Slant stack from incomplete gather b). (f) Residuals (d-e) (g) Velocity model obtained from e) by inverse elliptical RT. (h) Forward ERT from g). (i) Residuals (d-h). Chapter 4. Time variant Radon transforms offset(m) (c) Figure 4.7: (a) C D P from marine data with Transform with irregular Radon space, (c) model in b). (d) Difference between original generated from b). 130 q (ST 'm 2 ) x1 (d) three artificial gaps, (b) Hyperbolic Radon Original and interpolated traces using the (complete set) and synthetic traces (Gather Chapter 4. Time variant Radon transforms 131 offset(m) -1500 -1000 -500 offset(m) . -1500 -1000 -500 4.0-r (a)data (b)PRT-HRF (c)HRT-TD offset(m) -3000 -2500 -2000 -1500 4.0T offset(m) -3000 -2500 -2000 -1500 offset(m) -3000 -2500 -2000 -1500 (d)data (e)PRT-HRF ( f )HRT-TD offset(m) -4500 -4000 offset(m) -4500 -4000 -3500 offset(m) -4500 -4000 -3500 (g)data (h)PRT-HRF (i)HRT-TD Figure 4.8: (a-b-c) Near offsets: Gap within -460 to 0 m (a) Original data, (b) Inter-polated with sparse PRT. (c) The same with HRT. (d-e-f) Middle offsets: Gap between -2600 to -2140 m (d) Original data, (e) Interpolated with sparse PRT. (f) The same with HRT. (g-h-i) Large offsets: Gap between -4740 to -4276 m (g) Original data, (h) Interpolated with sparse PRT. (i) The same with HRT. Figure 4.9: (a) C M P gather at location 1000. (b) Linear RT of (a), (c) C M P gather at location 1000 with gaps and one every two traces, d) Linear RT of (c). (e) Model obtained from (d) using elliptical RT. (f) Linear RT of (c) obtained from (e) Chapter 4. Time variant Radon transforms 133 Figure 4.10: (a) Marine C D P . (b) Primaries obtained by substraction. (c) Multiples obtained from inverse R T of (e). (d) Sparse H R T . (e) R T after mute of primaries. Bottom - semblance analysis for: (f) data, (g) primaries, and (h) multiples Chapter 4. Time variant Radon transforms 134 cdp# 1000 1500 2000 2500 H RT- N M O-stac k Figure 4.11: (a) The Mississipi Canyon data after multiple removal with HRT, N M O correction and stack. Chapter 5 A hybr id Linear -Hyperbol ic R a d o n transform Where Horse and Rider each can trust the other everywhere, It takes a fence and more than a fence to pound that happy pair; For the one will do what the other demands, although he is beaten and blown, And when it is done, they can live through a run that neither could face alone. Rudyard Kipling As we have seen in previous chapters, one application of the Radon transform (RT) is to map events with different curvature in the data space to different areas of a model space, where these events can be separated or filtered before returning to the data space. This ability of the R T is in practice diminished by the fact that data contain events with different shapes. Linear events due to ground roll (GR), for example, are superimposed on hyperbolic events from reflections. Both events cannot be simultaneously focused with the standard R T . A linear R T localizes the G R and not the reflections, and a parabolic R T for N M O corrected data does not localize the G R . In this chapter I present a hybrid Radon transform to perform a sparse representation in the model space of data containing events with different shape. The R T is performed by two operators, each one with a different set of basis functions. The same idea can be extended to more than two operators for data with more complicated combinations of events. 5.1 The hybr id operator The Radon transform (RT) is usually implemented to focus events with one of three different shapes: linear, parabolic and hyperbolic. The linear R T (LRT) has been used 135 Chapter 5. A hybrid Linear-Hyperbolic Radon transform 136 mainly for processing directly in the slant stack domain (for example multiple attenuation, stack and migration). The parabolic and hyperbolic R T ( P R T and H R T ) are used for multiple attenuation and sometimes for interpolation and aperture extension. In each case there is an operator that maps the model space to the data space, and its adjoint (or conjugate transpose) that performs the reverse mapping, from data to model. Whatever the shape of the basis functions, the Radon Transform (RT) can be defined in terms of summation along paths, and represented in matrix form as where d and m represent the data gather and the Radon transform (model), respectively. The operator L can be applied frequency by frequency when the curve is time invariant (linear, parabolic, and pseudo-hyperbolic R T ) . Conversely, the time variant H R T operator has to be applied in the time domain, usually using sparse matrix techniques. Strong linear events are often present in the data superimposed on the hyperbolic reflections. A n appropriate forward modeling of these data is given by d = L m (5.1) d = L i i r i i + L 2 n i 2 (5.2) or equivalently d = [Lj L 2 m i (5.3) m 2 where the forward operator is [ L i L 2 ] , and the adjoint operator is . These two operators, forward and adjoint, are all we need to solve for the model m i using, for m 2 Chapter 5. A hybrid Linear-Hyperbolic Radon transform 137 example, a least square conjugate gradient algorithm (Claerbout, 1992). The two com-ponents of the model, m l and m 2 , are the linear and pseudohyperbolic Radon transform spaces. In the examples in this thesis I use L i as the Linear R T operator that maps lines to points, and L 2 as the pseudo hyperbolic Radon transform ( P H R T ) operator of Fos-ter and Mosher (1992). This second operator has basis functions following the equation t = T + pg(h), where g(h) is a function of offset and depth. Once a particular depth is chosen, the P H R T maps hyperbolas with apex at this particular depth to points. Deeper hyperbolas in the data transform to short ellipses in the model and shallower hyperbolas in the data to short hyperbolas in the model (Foster and Mosher, 1992; Bickel, 2000). Looking at the general system in equation (5.2) we could think of any combination of Radon operators. This kind of hybrid R T could be applied whenever the data have events with different shapes. Similar approaches are used to decompose images in a dictionary of basis or atoms, called basis pursuit, where the decomposition is done using an opti-mization method with constraints given by a lx norm for the model, leading to sparse representations (Chen et al., 1998) In this chapter I apply a hybrid operator for ground roll removal to understand whether it is possible to map different events to different spaces (Trad et al., 2001a). We will see that this is indeed possible, as long as the two spaces do not overlay, in other words, the basis functions are quite different. In this case, a very simple filtering procedure can be applied by inverting the data for the combined model (with all the operators), and reconstructing the data with only the desired operator. When the op-erators overlap, the same event can be mapped to both spaces and the filtering has to be applied more carefully in each model space to eliminate or separate the events. A multiple representation of the same data event in different model areas is undesired. A further improvement is to focus the solution using model weights. These weights allow the mapping of a particular event with preference for one of the operators. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 138 5.2 In quest of a zero nulspace The process of transforming data backwards and forwards often implies a loss of informa-tion, defined by the nullspace of the transformation. Sometimes this is a useful property to eliminate certain kind of events, which is automatically done when undesired com-ponents live in the nullspace of the transformation. However, this kind of uncontrolled filtering can also lead to loss of useful signal. Hence, in general, it is better to have a process that allows the recovery of data exactly, implying zero nullspace, and to apply a controlled filtering in the model space. Even when this is possible in theory, it is often not possible in practice, as the limited aperture and sampling of the model prevent a perfect mapping between data and model spaces. Sometimes it is possible to obtain nearly a perfect mapping using a large and densely sampled model space. A n important issue when a transformation is implemented is to use a complete set of basis functions, able to map any possible events in the data. Sometimes this goal can be achieved more economically by using different kinds of basis functions at the same time. This is the main idea behind new techniques like decomposition of frames, matching pursuit, best orthogonal basis, and basis pursuit (Chen et al., 1998). Another key property of a transform is the separability of the different events. We want to map different events to different areas of the model space. The separability is directly related to the sampling and aperture of the data, because imperfect information yields a lack of resolution. We have already seen that the lack of resolution can often be diminished by appropriate constraints on the model space, or regularization. From the previous discussion, it is clear that the path to follow is to design a transform, with a wide availability of basis functions suitable for a sparse representation of different events, and apply some sparseness constraint to overcome the lack of information in the model space. M y approach to this problem consists of the implementation of a hybrid Chapter 5. A hybrid Linear-Hyperbolic Radon transform 139 Figure 5.1: Hybrid inversion with sparseness constraints Radon transform to map with maximum sparseness two different kinds of events: linear and hyperbolic. When applying the Radon transform to data that contain events with these or other different shapes we want to map them to different spaces and work in each space separately. Using a combined R T is useful because, given a limited number of traces for the model space, the nullspace of the transform is reduced considerably. In other words, the combined R T requires a smaller number of traces to map the data, than a simple R T does. A n example of this situation is a linear R T that can map any kind of event if the range of slopes scanned is large enough, which implies that the model space has enough aperture and is sampled appropriately. However, by using a H R T and L R T simultaneously, the model space can be reduced, because rather than mapping hyperbolas to ellipses and lines to points, the combined L R T - H R T will map all events to points. A n example of this can be seen in Figures 5.2a-d and 5.3a-d. In Figure 5.3 we can compare the residuals between the original data (Figure 5.2a), and the recovered gather after mapping to the parabolic, linear and the hybrid (linear-hyperbolic) Radon Chapter 5. A hybrid Linear-Hyperbolic Radon transform 140 domains shown in Figures 5.2b, 5.2c and 5.2d, respectively. In all the examples I use the same number of traces in the Radon space for a fair comparison. The data, a shot gather from a vibroseis source [the number 10 in the data examples from Yilmaz (1987)], contain a superposition of different events, mainly linear and hyperbolic: two strong shallow reflectors, ground roll and side-scattered energy. A parabolic R T cannot map the data completely, as can be seen after recovering the data in Figure 5.2b. Many linear events disappear and several artifacts are introduced as a consequence of an incomplete representation, and hence, incomplete destructive interference. A linear R T does not produce the complete representation either, as can be seen by the hyperbolic events at near zero offset in the residuals (Figure 5.3c). Extending the traces to lower slowness values would completely map the data, with the hyperbolas mapped to ellipses, but it would not be possible to separate the linear events without touching the elliptical curves because of aliasing problems. A hybrid Radon is able to map most of the events (Figures 5.2d and 5.3d) and, as we will see later, it is very easy to separate linear and hyperbolic events automatically. 5.3 Inversion for sparse hybrid models Similar to the example presented in Chapter 2, I solve the constrained inversion by first defining a cost or objective function. The inverse Radon operator is found by minimizing the cost function J = | |W d (Limi + L 2 m 2 - d)||2 + | | W m l nn|^ + | | W m 2 m 2 | | 2 , (5.4) where Wd is a matrix of diagonal weights containing the inverse of the standard deviation ofthe data and W m i , W m 2 are matrices of model space weights that we can set to enhance the resolution of the R T . The notation becomes simpler by defining, L = [L 1 L 2 ] , (5.5) Chapter 5. A hybrid Linear-Hyperbolic Radon transform 141 m = mi m 2 and W = v v TO — w m l 0 o w m 2 The minimization of the cost function J leads to the Radon model, m = ( W £ W m + L T W j W d L ) " 1 L T W j W r f d . W m i , W m 2 are diagonal matrices, whose elements are set to (5.6) (5.7) (5.8) /rrii (5.9) such that l|m||i = £ K | = m T W^W m m = ||Wmm||2. (5.10) i The model that minimizes the cost function (5.4) is solved iteratively by using the IRLS method explained in chapter 2. 5.4 Examples There are several methods for ground roll removal, and it is not my intention to design a new one. I believe that the results presented in this chapter can be improved by some other filtering techniques as, for example, a wavelet transform with hyperbolic two dimensional wavelets, or some clever f-k filtering (although my attempts at f-k filtering were not completely successfully because of the aliasing). Nevertheless, the ground roll problem is appropriate to test a hybrid R T . M y intention is to investigate whether the Radon transform can be enhanced by combining different operators. I designed a simple synthetic example, where five hyperbolic reflections are crossed by six linear low frequency, high amplitude events (Figure 5.4a). A low pass filter would Chapter 5. A hybrid Linear-Hyperbolic Radon transform 142 remove the noise completely, together with the low frequency component of the signal. A f-k filter is able actually to partly separate the two kind of events (Figure 5.4c). However, aliasing is conspiring against this procedure, because the alias of the linear events leaks over to the signal region in the f-k domain (Figure 5.4b), and applying a filter removes only the non-aliased part of the signal (Figure 5.4d), leaving in the data the aliased part of the linear events. In applying the hybrid R T (Figure 5.5), the first operator maps the linear events to the first area of the model, the linear Radon domain ( L R D ) , while the second operator maps the hyperbolic events to the pseudo hyperbolic space, the hyperbolic Radon domain (HRD) (Figure 5.5b). To avoid distortion of the linear events, I do not apply any N M O correction, therefore a P H R T approximates better the reflections than a parabolic R T . A mute applied in the hyperbolic space removes the signal, and after an inverse R T is applied, only the linear events are predicted (Figure 5.5c). The subtraction of the linear events from the data leaves the hyperbolic events intact. Only a small amount of noise remains, in particular at the top of the linear events (Figure 5.5d), probably because of non perfect destructive interference of the noise. In the second example I remove ground roll noise in a real data set, the shot gather number 25 from Yilmaz (1987) (Figure 5.6a). Because a shot gather is in general asym-metric, any R T other than linear requires a different model for negative and positive offsets (this is another problem where a R T with two different operators can be used). Figure 5.6b shows 320 traces in the Radon domain. The first 160 traces correspond to the negative offsets, the 161-320 traces correspond to the positive offsets. Each half shows the two spaces side by side, left hand side is the linear R T , right hand side is the pseudo hyperbolic. Applying again a mute in the pseudo-hyperbolic space, or equivalently, using only the linear operator of the hybrid Radon transform, the linear ground roll can be modeled (Figure 5.6c) and separated by subtraction (Figure 5.6d). A very small part of the ground roll remains in the small offset traces, and can perhaps, be removed using Chapter 5. A hybrid Linear-Hyperbolic Radon transform 143 additional filtering. This part of the data has not been predicted for two reasons: the slopes become very large (approaching infinity) and the linear events are heavily aliased. In fact f-k filtering has similar difficulty in removing the ground roll in the near zero offset traces. The last example (Figure 5.7), is intended to show the importance of the sparseness constraint in the localization of events. The synthetic gather has some hyperbolic events together with two linear ones (Figure 5.7a). The slopes of the linear events are inside the range of the pseudo hyperbolic domain, hence the R T will rely entirely on the sparseness constraint to focus the model of the linear events to the linear space. The non sparse R T fits the data with no preference for the linear or the hyperbolic space, and the event is mapped in part to the Linear Radon domain (LRD) , and in part to the Hyperbolic Radon domain (HRD) (Figure 5.7b). Figure 5.7c shows the spectrum (f-p-q, where p is slowness and q is a curvature parameter) of the R T . Figure 5.7d is the recovered data from the hyperbolic part of the model space only. We see that part of the linear events leaks into this space. Figure 5.7e shows the recovered data from the linear space only. No signs of the reflections can be seen here, but the linear events are not complete. Figure 5.8 contains the same plots but this time using the sparse R T . The linear events are mapped into the L R D , where they have smaller size in the t\ norm (Figures 5.8b and c). Note that the linear and hyperbolic events have been separated almost completely (Figures 5.8d and e). 5.5 Discussion and conclusions In this chapter a combined linear - pseudo hyperbolic Radon transform has been intro-duced, the purpose of which is to decrease the size of the model space without increasing the nullspace of the transform. Applications of this hybrid approach could be coherent Chapter 5. A hybrid Linear-Hyperbolic Radon transform 144 and/or incoherent noise removal. A separation can be obtained in reasonable time using an iterative re-weighted least squares method, and, if possible, by choosing a range of parameters for the operators that produce different basis functions. Real and synthetic examples show a successful separation of events with different shape into the two Radon domains. A n advantage of using a hybrid approach is the freedom to choose different parameters for the two parts of the model space, for example, sampling interval in the Radon domain, maximum frequency to use, regular or irregular Radon sampling, etc. If the data contain linear and hyperbolic events, the sampling for the hyperbolic space cannot fulfill the aliasing condition for linear and hyperbolic events simultaneously. The hybrid R T does not have this problem; each part of the model space fulfills the alias condition for the events that it is expected to map. Although a two operator approach is quite simple, the sparseness constraint is es-sential to prevent the linear event being mapped as an incomplete hyperbolic event and vice versa. If a simple £ 2 constraint is used, the method resembles the method of frames, where the solution is an average of all possible solutions. If an t\ constraint is imposed, a theoretically perfect separation is possible. In this case the method resembles a basis pursuit approach (Chen et al., 1998). As the data are scanned for all the basis functions, the separation is more difficult when the same event can map parsimoniously to different areas of the model space. For example I tried, unsuccessfully, to apply the idea to two pseudo hyperbolic operators with different focusing depths. Model weights may be able to overcome this problem, but I have not been completely successful in this regard. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 145 offset(m) offset(m) 1000 2000 3000 1000 2000 3000 (a) (b) offset(m) offset(m) 1000 2000 3000 1000 2000 3000 (c) (d) Figure 5.2: (a) Shot gather, (b) Recovered data after Parabolic RT. (c) Recovered data after LRT. (d) Recovered gather after LRT + PRT. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 146 Figure 5.3: (a) Shot gather, (b) Residuals after Parabolic RT. (d) Residuals after L R T + PRT. (c) Residuals after LRT. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 147 Figure 5.4: (a) Synthetic gather with hyperbolic and linear events, (b) Data in the / — kx domain, (c) / - k Filtered data, (d) Filtered f — kx domain. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 148 Figure 5.5: (a) Synthetic gather with hyperbolic and linear events, (b) Sparse lin-ear-hyperbolic RT. Horizontal axis shows trace numbers. Slopes for the left hand side (linear space) and curvatures for the right hand side increase to the right side, (c) Pre-dicted noise obtained by muting the parabolic space, (d) Signal obtained by subtracting c) from a). Chapter 5. A hybrid Linear-Hyperbolic Radon transform 149 Figure 5.6: (a) Shot gather, (b) Combined Linear + Pseudo Hyperbolic Radon domain, for negative offsets (traces 1-160), and for positive offsets (traces 161-320). (c) Recovered gather from linear space only, (d) Signal obtained by difference between (a) and (b). Chapter 5. A hybrid Linear-Hyperbolic Radon transform 150 Figure 5.7: (a) Synthetic gather with hyperbolic and linear events, (b) Non-Sparse Linear-hyperbolic RT. (c) Spectrum of (b). (d) Recovered data from hyperbolic model, (e) Recovered data from linear model. Chapter 5. A hybrid Linear-Hyperbolic Radon transform 151 offset(m) #trace #trace (d) (e) Figure 5.8: (a) Synthetic gather with hyperbolic and linear events, (b) Sparse Lin-ear-hyperbolic RT. (c) Spectrum of (b). (d) Recovered data from hyperbolic model, (e) Recovered data from linear model. Chapter 6 Aliasing It is commonly assumed that given a sampling interval, the corresponding Nyquist fre-quency imposes a strict limitation on the maximum frequency allowed in a sampled signal. There are several situations where this is not true. It is often possible to recover information above the Nyquist limit when facing irregular sampling, high dimensionality, availability of prior information and sparse inversion. Aliasing effects for the R T can be explained in terms of summations in phase per-formed by the discrete R T operator. For some operators, like the L R T , irregular sam-pling leads to some degree of alias attenuation. In seismic processing, irregular sampling is sometimes imposed by practical considerations. Therefore, irregular sampling can be used to attenuate aliasing for the L R T . Even when further processing requires regular sampling in offset, modeling in the Radon domain allows for re-sampling to a regular grid. Another situation where aliasing is attenuated is when the solution is constrained, such that aliased events are penalized. This is the case for the time domain sparse R T . 6.1 Aliasing in diferent flavors The R T shares the same aliasing problems as the Fourier transform ( F T ) . Usually low velocity linear events appear as a consequence of air waves, or ground roll. In such cases, where the linear R T is a useful tool in separating low velocity linear events from the data, we have a serious constraint given by aliasing. Any attempt to remove these low velocity events must also mute the aliased events that for these situations, could be located in 152 Chapter 6. Aliasing 153 the space of useful information. In the P R T aliasing prevents the complete elimination of multiples. In seismic processing a distinction is often made between data and operator aliasing. Data aliasing is understood as the effect of under sampling a signal with respect to the Nyquist frequency, either along the time or space coordinates. The term operator aliasing defines a situation where the action of the operator on the data introduces aliasing artifacts, independently of the existence of aliasing of the data. This type of aliasing can often be attenuated by careful implementations of the operators, but first it is necessary to understand the origin of this aliasing. In this sense it is very useful to analyze the effect of the impulse responses of the forward and inverse operators. Aliasing arises in the Radon transform space due to poor sampling in offset, which is a common problem. Aliasing in the data space due to poor sampling in the transform space is also possible, but this problem is easy to correct because it simply requires a decrease in the sampling interval for the model space. However, when the operator is inverted to obtain the transform, decreasing the model space sampling considerably increases the computational cost of the transformation and requires careful regularization. When the offset is irregularly sampled, the conditions under which aliasing appears are not very clear. Some interesting comments about the effect of irregular sampling on aliasing, in the context of N M O and stacking, can be found in Wisecup (1998). This paper emphasises the difference between input driven processes that produce irregular output sampling, and output driven processes that produce regular output sampling. We can certainly look at this phenomenon in the context of the Radon transform. Chapter 6. Aliasing 154 ( a ) ( b ) 0.8 0.6 0.4 0.2 0 I 1 0.8 0.6 0.4 0.2 0 200 400 ( c ) 600 u JL u 0.5 ( d ) 500 400 300 200 100 200 40C 600 .5 0 0.5 Figure 6.1: (a) Regular sampling function, (b) Fourier transform of sampling function, (c) Irregular sampling function, (d) Fourier Transform of Sampling Function 6.2 Diferent views of aliasing For a given continuous signal we avoid aliasing by discretizing the signal in such a way that all of the frequencies contained in the signal fall below the Nyquist limit. In fact it seems that embedded in our heads is Mr Nyquist saying, " If you cross my limit you will fall into a trap of chaotic aliased data". Indeed, when Christopher Columbus decided to sail around the earth, people said he would fall off its edge. However, Columbus found that the earth was not as one-dimensional as people thought it was, just as the aliasing problem is not as one-dimensional as people are led to believe. This is literally true when we move from one-dimensional to two-dimensional transforms. Data aliasing refers to aliasing with respect to the Fourier transform operator. It Chapter 6. Aliasing 155 can be easily explained by noting that sampling a continuous function is equivalent to multiply the function by a Dirac comb function (a bed of nails). Noting the equivalence of element-wise multiplication in the time domain and convolution in the frequency domain, we see that the F T of the digitized signal will consist of repeated copies of the F T for the continuous signal. The reason is that the F T of the sampling function is also a spiky function (Figure 6.1a-b). Therefore aliasing appears because of the overlapping of copies of the same spectrum due to the convolution. In mathematical terms, a sampling function with a period of A s , has a sampling period in the Fourier domain of Ks = 1 /As . Hence, if the bandwidth of the signal Kx = kXmax — kXmin (where kx is the F T variable of x) is greater than Ks, then the signal is aliased. Information theory provides us with a tradeoff relationship between the bandwidth of a signal and its sampling rate. In order to avoid aliasing, the bandwidth Kx of a function f(x) sampled at a unit sampling interval, must be less than 2TT, which can be easily visualized by taking the z transform of the signal. If the sampling interval is not one, then we have to scale the bandwidth by the actual sampling rate Ax. Therefore, the bandwidth must satisfy AxKx < 2vr, (6.1) and the sampling rate must be A x < 2_" • (6.2) This condition can also be expressed in terms of the slope, p, of a linear event and its frequency. Since kXmax - upmax and kXmin = copmin, then Ax < -j—^L - , (6.3) ^ \Pmax Pmin ) which is the maximum offset interval for the L R T . The maximum offset interval for the P R T can be obtained using the stretching transformation x —> x2 and the approximation Chapter 6. Aliasing 156 Ax = 2xAx. By replacing these variables into equation (6.3), we have A x < 9 ~ m — T ~ — z — v ( 6 - 4 ) Let us examine the maximum sampling intervals required in the Radon space. Given the projection operator L and its adjoint LH, the adjoint model madj is given by m a d j = LHd = L H L m t r u e . (6.5) This adjoint model is the result of convolving the actual model m t r u e with the smearing filter a = LHL whose expression is for the L R T , and ,(u,p)= [ e-^xdx, (6.6) Jx .(u,p) = I e-l"px2dx, (6.7) Jx for the P R T . Hugonnet and Canadas (1995) derived the maximum sampling rate intervals for the Radon parameters for both the L R T and P R T . They used the convolution/deconvolution formalism in the frequency domain, calculating the bandwidth of the operator itself and using equation (6.1). For example, the bandwidth of the smearing filter aun can be estimated as follows. The F T of a function f(p) with respect to p is Fs{kp) = [ s{p)e-lkvpdp. (6.8) Jp The spectrum for the smearing filter aiin is Falin(co, kp)= f [ e ^ - ^ d p d x = [ 2TT5{UJX - kp)dx (6.9) Jx Jp Jx Chapter 6. Aliasing 157 which is f 2ir/iU if COXmin < kp < cox Faiin(u,kp) = < (6.10) 0 else Hence, the bandwidth for the operator is Kp = co(xmax — x m i n ) . From equation (6.1) the critical sampling interval is A p < n x '-,•)• ( 6 1 1 ) Applying a similar derivation to the smearing filter apar Hugonnet and Canadas (1995) show that J \ max ^minJ Equation (6.3) can be derivated in a similar way, using the filter LLF, but equation (6.4) is more difficult to obtain because it involves the solution of an integral without an analytical solution. Now let us look at what happens when these conditions are not fulfilled. The situation is more complicated than for the F T because the spacing of the aliased events in the Radon domain is frequency dependent, and therefore, the aliased events pa for an event p t r u e are located at a frequency dependent slope (Turner, 1990) P ( / ) = J W + (6-13) where n is a positive or negative integer. For example, to avoid aliasing we need to keep the maximum slope in our data below the value given by Pmax ~ "7 T I" Pmin- (6-14) Jmax^--^ If we have an event with a slope ptrue greater than pmax = l/(Axf) + pmin, the aliased event will appear first. If the slope p t r u e is smaller than pmax but the Radon space goes beyond p(f) = ptTae + n/(Axf) again the aliased event will appear. To avoid interference Chapter 6. Aliasing 158 between true events and aliased events we must keep the Radon space between limits imposed by the aliased events (Marfurt et al., 1996). For the parabolic R T , the aliasing conditions are only approximated because the sampling is now irregular. There are local summations in phase for 77 ltf)=^e+2A _ y (6.15) In other words, the "pNyquist" o r uQNyquist" are different at every frequency. Hence, an aliased event will not be localized in the time domain, but will be dispersed throughout the whole Radon domain. 6.3 O p e r a t o r a l i a s i n g a n d i n v e r s i o n Operator aliasing can also be explained in the context of inversion as an effect of the non-uniqueness of the model. Aliasing in this context arises because the operator allows a single event in the true model to produce several events in the estimated model. Note in equation (6.5) that the true model is free of aliasing, and consequently, aliasing is introduced into the adjoint model by the smearing filter. Therefore, the deconvolution of the smearing filter should, in theory, remove the aliasing. However, aliasing remains because the deconvolution of the smearing filter is not perfect due to the finite aperture and discretization. More importantly, perfect deconvolution of the smearing filter is not possible because aliasing introduces ill conditioning in the kernel of the transformation. The effects of the undesired spreading function filter can be eliminated from the adjoint madj by (Tarantola, 1987) m = (LHC^L + C " 1 ) - 1 ! ! ! ^ , (6.16) where, as before, C m and are the model and data covariance matrices. Any difference Chapter 6. Aliasing 159 between the true and recovered models (for example aliasing) are given by mdiff = (I - R)mtrue (6.17) where I is the identity matrix and R is the resolution matrix given by R=(L^ C,-1L + C-1)-1L^C,-1L, (6.18) such that m = Rm, ltrue- (6.19) Hence, the remaining aliasing depends on the invertibility of the smearing filter given by the matrix LHL. If this matrix is well conditioned, then = I, C " 1 = 0 (matrix of zeros) and the true model is recovered exactly. When zero order regularization is applied, the larger the regularization the larger the remaining artifacts. For the L R T with regular offset sampling, the spreading function matrix has lower condition number as the frequency increases, but after aliasing appears becomes ill con-ditioned (Sacchi, 1995). Regularization is now required, and the aliasing effects remain. Schonewille and Djuindam (2001) show that the increasing ofthe condition number when the frequency decreases is a consequence of using constant Radon sampling for all frequen-cies. At low frequencies, the condition number is large because there is an oversampling of the R T . Because this sampling interval is calculated for the maximum frequency of in-terest, the applied sampling interval is optimal only for this frequency, where as a result, the condition number is the lowest. If the Radon sampling is chosen independently for every frequency, then the condition number of the matrix is constant for all frequencies and can be set to the minimum subject to equations (6.11) or (6.12). However, after aliasing appears, the condition number increases and careful regularization is required. Figure 6.2a shows the decreasing of the condition number for the (regularized) smear-ing filter, when dh = 25m and pmax = l x 10"3, and fnyquist = l/{Ahpmax) = 40Hz. Chapter 6. Aliasing 160 0 10 20 30 40 50 60 70 80 90 100 frequency Frequency Figure 6.2: Condition number as a function of frequency, a) L R T . Note the jump at / = 40Hz where aliasing arises - b) P R T When aliasing appears, the condition number increases and only regularization keeps it at reasonable values. The P R T instead shows a smoother behavior in the condition number (Figure 6.2b). 6.4 Irregular sampling A first inspection of the effect of irregular sampling for the F T operator is a good be-ginning to understand the effect of irregular sampling in general. Because the F T of an irregular sampling function is a spike with a background level (Figure 6.1c-d), when a continuous signal is digitized with an irregular sampling function, the resulting spectrum does not contain repeated copies of the same spectrum, but a distortion due the non null background. Hence, we see that irregular sampling can eliminate the aliasing for the F T operator. This concept is very well known, but not very often considered as particularly useful, as most of the time we use regular sampling in the time domain. Also, for the Chapter 6. Aliasing 161 (a) 300 200 r i 100 o -100 300 200 r | 100 0 -100 I I 1 I i i i i -4 -3 0 (b) 4 5 x10~3 I I I -I I I Figure 6.3: Smearing filter L^L (linear r — p ) a) / = 10Hz, regular offset no aliasing, b) / = 2bHz, regular offset, aliasing, (c) / = 25Hz, irregular offset, no aliasing. F T , irregular sampling brings more disadvantages than benefits. The new F T operator is not orthogonal anymore and a full inversion is required to recover the data. Further, the efficient F F T algorithm cannot be used. Because the effect of irregular sampling is different for different operators, an analysis is required in each case. Let us examine the smearing filter for the L R T operator in the cases of regular and irregular offset sampling. Figure 6.3 shows three situations for the smearing filter for the L R T . Figure 6.3a shows the smearing filter at a frequency small enough such that there is no aliasing. As we move to higher frequencies or larger Ah, aliasing appears (Figure 6.3b). The smearing filter causes every real event with a particular slope to produce three different slopes in the adjoint model. Deconvolution of Chapter 6. Aliasing 162 0 0.2 0.4 0.6 0.8 1 P X10"3 Figure 6.4: Resolution and smearing niters (linear R T ) a) resolution filter from smearing filter b) (regular offset) - c) resolution filter from smearing filter d) (irregular offset) this filter is ill posed and aliased events will remain. Figure 6.3c shows the filter resulting from the same p values and frequency, but from irregular offset sampling. The span in offset and the number of offset traces are the same, but the offset spacing is not constant. The aliasing disappears, although the filter becomes more irregular. Figure 6.4 shows an example of the resolution filter, equation (6.18), that connects the true model with the inverted model. Figure 6.4a shows the resolution filter obtained from the smearing filter in Figure 6.4b for a evenly sampled offset case. Figure 6.4c is the resolution filter obtained from Figure 6.4d for irregular offset sampling. Similarly to what happens with the F T , if we use irregular sampling, the aliased events are attenuated. However, the consequences are more interesting because in the Chapter 6. Aliasing 163 (a) _ 5 Q I I I i I I 1 I i 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 « X10"6 Figure 6.5: Smearing filter LHL (parabolic R T ) a) / = 10Hz regular offset - b) / = 25Hz regular offset -c) / = 25Hz irregular offset R T the non-orthogonality is always present and deconvolution is required whether the sampling in offset is regular or not. In the R T , therefore, advantages of regular sampling are not nearly as important as in the F T case. For the parabolic R T , unfortunately, irregular offset sampling is not as efficient to attenuate aliasing as it is for the L R T . The reason is that aliasing for the P R T arises as a consequence of local summations in phase, and irregular sampling produces little difference for the short offset distances involved in these local summations. Figure 6.5a shows the smearing filter without aliasing. For large frequencies the filter appears wildly oscillatory for large values of q and aliasing appears as we see in Figure 6.5b. A n event with a particular curvature given by q will produce a main peak on the adjoint, but also Chapter 6. Aliasing 164 smaller aliased events (secondary lobes) will appear at other q values. As a result the adjoint model will contain mixed information from the true model. Figure 6.5c shows the filter resulting from the same q values and frequency, but from irregular offset. The lateral lobes have decreased a little and again the filter has become more irregular. As the frequency becomes larger, the decreasing of the lobes is less important and the distortion due to aliasing remains. Aliasing produces artifacts in the model space and a loss of resolution because the high frequencies are smeared in the Radon domain. 6.5 A l i a s at tenuation by the t ime domain sparse R T There is a difference between aliased and true events that allows for alias attenuation for both the L R T and P R T . Aliased events are dispersive in the time domain because of their frequency dependency, whereas true events are localized in the time domain, because they have all their frequency components located at the same slope or curvature. Using this characteristic we can attenuate aliasing by using a regularization method leading to sparse models in the time domain (Marfurt and Nemeth, 2001). A time domain sparse R T will fit the data with a few large elements (non aliased components) rather than many small elements (alias artifacts). Another way to understand how a sparse constraint in the time domain can attenuate aliasing is to consider crossing events in the frequency offset domain. These crossings are a problem for the frequency domain R T because of the uncertainty in the relative contributions from the aliased and primary events. Because every frequency is solved separately, there is no a priori information about the relative amplitudes. By imposing a sparseness constraint in the time domain, the spectrum is assumed to be smooth, and therefore the uncertainty disappears by using the spectrum information from neighboring points at other frequencies. Chapter 6. Aliasing 165 200 400 Time ( d ) 0.25 0.2 | 0.15 a I 0.1 0.05 ililiiJiiliJl' lLl-^- l i^^Jll^IJIt^ -1 -0.5 0.5 Freq F i g u r e 6.6: (a) Regular ly sampled t ime series w i t h f\ = 0.35Hz, f2 = 0.7Hz, (b) Irregularly sampled t ime series, (c) Fourier transform of t ime series i n (a), (d) Fourier transform of t ime series i n (b) 6.6 Examples L e t us ana lyze the s p e c t r u m for a s imple t i m e series, composed o f t h e s u m of two h a r m o n i c s ignals w i t h frequencies g iven by f\ = 0.35Hz and f2 = 0.7Hz. I d i g i t i z e d the t i m e series at regular in tervals (F igu re 6.6a) and at in tervals d r a w n f rom a n o r m a l d i s t r i b u t i o n (F igu re 6.6b). In order to o b t a i n an e s t ima t ion of the s p e c t r u m b e y o n d N y q u i s t I use a F o u r i e r t r ans fo rm ob ta ined by least squares invers ion ins tead of the F F T . T h e c o m p u t e d s p e c t r u m f rom the regular s a m p l i n g is a l iased (F igu re 6.6c), as the 0.7 H z frequency appears at 0.7Hz a n d at fa = 2fN — f2 = 2 x 0.5 - 0.7 = 0.3Hz. A l s o , a n a l i ased vers ion of / x = 0.35Hz appears at fa = 2fN — fx = 2 x 0.5 — 0.35 = O.QhHz. T h u s , we see two Chapter 6. Aliasing 166 true frequencies and two aliased frequencies. The spectrum of the irregularly sampled data, on the other hand, shows only the two true frequencies (Figure 6.6d). Usually we look at the spectrum only until the Nyquist frequency, so that only two frequencies, fi = 0.3bHZ and fa = 0.3Hz, appear in the standard F F T . The second example shows the same phenomenon in the linear R T . I have created a synthetic data set with two linear events. The first event at 0.4s simulates the air wave with a slope p = 3.10~ 3s/m. The second event at 0.8s has a slope p — 4.10~ 4 s/m. Figure 6.7a shows the original regularly spaced data set. For this example fmax = 70Hz and dx = 10m, so that PNyquist given by Equation (6.14) is p = 0.0014s/m, which is less than p = 0.003s/m from the air wave. As a consequence the low velocity event will be mapped to the corresponding p value (0.003 s/m), but first it will appear at the aliased location pa = 0.003 - PNyquist-Because the aliased locations are frequency dependent we cannot localize the aliased event in the time-Radon domain, Figure 6.8a, but we can still visualize the problem in the frequency-Radon domain (Figure 6.7c). It is clear in this figure that the aliased events are outside of the chosen Radon subspace for low frequencies, but as we move to larger frequencies the alias appear in the R T space because the separation between pauased and p t r u e decreases. Note that in the spectrum of Figure 6.7b, the energy following a curve proportional to 1 / / at the left side is the aliased event for p t r u e = 0.003s/m, and the energy following a curve at the right side is the aliased event for the high velocity event with Ptrue — 4 x 1 0 - 4 s / m . The elements of the model are very localized in the frequency domain because I have used a sparse (frequency domain) L R T . If we apply a mute to remove the air wave, the higher frequency components of the aliased events will remain and the event will still appear after reconstruction with the inverse Radon transform, Figure 6.8c. Now let us see how the same procedure works for the same model but sampled at non Chapter 6. Aliasing 167 regularly spaced offsets, Figure 6.7b. I have used an offset axis with a mean dx = 10m as before, but now with a perturbation normally distributed. The irregular sampling does not produce aliasing, as can be seen in Figure 6.7d, and the same muting in the Radon domain, Figure 6.8b, separates the events very well (Figure 6.8d). It is clear from these plots that the mute in the Radon domain for the irregular sampled data will produce a better separation between high and low velocity events, and so we can extend our Radon space without limit, never seeing aliased events. A similar alias attenuation can be obtained by applying sparseness in the time domain. Figure 6.9a shows the same aliased data as before and Figure 6.9b shows its / — k spectrum, where the aliasing for the high slope event is clearly seen. The non sparse R T of these data maps the two events to its correct slopes, but the aliased artifacts are spread throughout the space. These can be seen in the f —p domain (Figure 6.9d), where each event shows its main appearance and its dispersive alias with a location dependent on 1/ / . By muting slopes smaller than 0.002s/ra in the Radon domain, and then applying the inverse transform, we recover the larger slope event, and part of the alias of the small slope event. Therefore, the separation is not successful (Figure 6.9e). Figure 6.9f shows the sparse R T applied in the time domain. The aliasing has disappeared (Figure 6.9g) and now the same muting results in a successful separation. Let us see a similar example for the parabolic R T . To keep the discussion focused on aliasing, without considering N M O or time stretching, I have used parabolic events that should be mapped, in theory, to a single wavelet per event. The synthetic data contain two parabolic events, one at 0.4s with a velocity of 1500ra/s, and the other at 0.8s with a velocity of 2500m/s (Figure 6.10a). The / — k spectrum shows that the data are heavily aliased (Figure 6.10b). The standard non sparse R T produces the main events and their alias, that are dispersed throughout the model space (Figures 6.10c and d). Even though the alias events cannot be seen in the figures because they have low amplitudes, the Chapter 6. Aliasing 168 failure when trying the separation of the events by muting proves that the events are not localized. In contrast, the time domain sparse R T allows a complete separation of these events. The aliased events are highly diminished (Figures 6.10f and g). Muting events with q > 3 x 10~7s/m2 separates the two parabolas (Figure 6.10h). 6.7 Conclusions Aliasing has always been understood in terms of the Nyquist limit beyond which events are mixed and information is ambiguous. As most limits tend to do, however, the Nyquist limit appears to dissolve under certain conditions. For example, the Nyquist frequency and Nyquist slope are, indeed, limits for the F T and L R T respectively when regular sampling is enforced, but not when irregular sampling is in effect. In the linear Radon transform, where advantages of regular sampling are no longer so important, irregular sampling can be advantageous and computationally efficient. In other situations like the P R T , irregular sampling does not eliminate aliasing, but the dispersive nature of the aliased events allows for their attenuation using sparse time domain transforms. A time domain method will attenuate aliased events because of the highly dispersive nature of the aliased artifacts in the time domain. Some of the problems with data sampling could also be solved by using interpolation methods, such as prediction methods advocated by Spitz (1991). Wave equation methods themselves provide an implicit interpolation during the migration process (Zhu and Lines, 1997). There are many other situations where the aliasing conditions imposed by our incom-plete understanding of the problem are too strict. To develop precise aliasing conditions is an important practical aspect of seismic research because it leads to improvement in resolution, a decrease of computational load and a decrease in data gathering costs. Chapter 6. Aliasing 169 However, as the dimensionality of the problem increases, the complexity of aliasing phe nomena also increases and presents an interesting challenge for future research. Chapter 6. Aliasing 170 Figure 6.7: Two events sampled with (a) regular and (b) irregular sampling, (c) Radon domain (f-p) from (a), showing the first aliased event (left hand side) for the first high slope event and the first aliased event (right hand side) for the second lower slope event, (d) Radon domain (f-p) from (b), with attenuated alias events Chapter 6. Aliasing 171 p(s/m) p(s/m) 0.001 0.002 0.003 0.001 0.002 0.003 Figure 6.8: Filtering low velocity events in the Radon domain beyond alias from regular sampling data does not remove completely the event but it does almost completely if the data are irregularly sampled (a) Radon domain from regular sampled data (b) Radon domain from irregular sampled data (c) Filtered from (a) and (d) Filtered from (b) Chapter 6. Aliasing 172 offset(m) 200 400 0-2 0.4 p(s/m) 0.002 0.003 (c) R T - n o n - s p a r s e p(s/m) 0.001 0.002 0.003 (a) d a t a (b) f-k s p e c t r u m p(s/m) 0.001 0.002 (d) R T s p e c t r u m - n o n - s p a r s e p(s/m) 0.001 0.002 0.003 oflset(m) 200 400 (e) S e p a r a t e d e v e n t (f) R T - s p a r s e (g) R T s p e c t r u m - s p a r s e (h) S e p a r a t e d e v e n t Figure 6.9: (a) Synthetic aliased data, (b) f-k spectrum ofthe aliased data (c) Non sparse LRT. (d) Spectrum of (c). (e) Recovered data from (c) after muting slopes smaller than 0.002s/m. (f) Sparse L R T obtained with a time domain method, (g) Spectrum of (f). (h) Recovered data from (f) after the same muting. Chapter 6. Aliasing 173 (f) R T - s p a r s e (g) R T s p e c t r u m - s p a r s e (h) S e p a r a t e d e v e n t Figure 6.10: (a) Synthetic aliased data (parabolas), (b) f-k spectrum of the aliased data (c) Non sparse PRT. (d) Spectrum of (c). (e) Recovered data from (c) after muting events with Radon parameter smaller than 3.10 - 7 s/m2. (f) Sparse PRT obtained with a time domain method, (g) Spectrum of (f). (h) Recovered data from (f) after the same muting. Chapter 7 Radon transform applied for migration In this chapter I will merge the Radon transform with a time migration method derived from Kirchhoff migration. We have seen that the Radon transform is a tool to sepa-rate events according to moveout and it is often applied on common midpoint gathers (CMPs) . The moveout in C M P s is produced by differences in the traveltime because of offset variations. Kirchhoff migration has some resemblance with Radon transform in the fact that a stack is performed along a traveltime path, but the moveout is calculated as a function of midpoint displacement and offset, given by the double square root equation (DSR). Equivalent offset migration ( E O M ) is a similar stack along a traveltime path given by the DSR. However, the stack curve is parameterized differently, in such a way that the D S R becomes a single square root equation, similar to the traveltime function in the hyperbolic Radon transform. In this chapter I explore how the E O M method can be combined with the Radon transform to achieve a series of benefits: improved focus-ing in the final section by performing a sort of automatic correction of initial velocity estimates, elimination of coherent and incoherent noise from the data during the migra-tion process, and computation of C S P gathers that give more precise velocities using semblance analysis. I test the combined method with a synthetic 2D finite difference data set and a 2D marine data set. The results suggest important improvements inside reasonable computing time limits. 174 Chapter 7. Radon transform applied for migration 175 Figure 7.1: The prestack volume of seismic data and CSPs mapping. 7.1 EOM The method of Equivalent Offset Migration ( E O M ) maps the energy scattered at every point in the subsurface from the prestack volume (with offset, midpoint and time coor-dinates) to a plane of hyperbolic events called a common scattering point gather (CSP) (Bancroft et al., 1998) (Figure 7.1). A C S P gather at a given location has hyperbolic events corresponding to the constructive summation of energy originated at scatter points below that location, but also has non-coherent energy scattered at other points. These gathers can be used to improve velocity estimation by semblance analysis. A Kirchhoff-N M O process on the CSPs produces a stack-migrated section, where the Kirchhoff-NMO process consists of phase filtering, amplitude scaling, N M O correction and stack. The advantages of this method with respect to prestack Kirchhoff migration resides in faster computation and in the fact that velocities obtained from C S P gathers are more precise than those obtained from C M P s , because of a larger fold. Perhaps most important is the fact that CSPs produces root mean squares (RMS) velocities, close to Chapter 7. Radon transform applied for migration 176 migration velocities, whereas C M P s produce stacking velocities, largely affected by dips. From my point of view, another advantage is that the information inside a large volume of data can be mapped to a plane, facilitating further processing. A reduction of the data volume while preserving the moveout information allows one to apply processing tools that would be very expensive if applied to the whole prestack volume. The position of each seismic sample in a 2D seismic line is given by three coordinates, half offset, midpoint and time (h,x,t). These three coordinates define the prestack vol-ume that contains all the data. A point inside the earth scatters energy with a delay or traveltime that defines a reflection surface in the volume. This surface is called the Cheops Pyramid (Claerbout, 1985), and its mathematical formulation for an homoge-neous medium is the double square root equation (DSR) If the medium is not homogeneous, but the velocities change smoothly, the use of R M S velocities makes the D S R equation (7.1) a good approximation for the traveltimes. In Kirchhoff time migration, the position of a reflector as a function of depth (expressed in terms of two way traveltime) is obtained by constructive interference of the energy scattered by the points on the reflector. The energy recorded at seismic traces inside the prestack volume is stacked following the traveltimes obtained by the D S R . Time Kirchhoff migration takes every sample from the prestack volume and moves it back to the position of the scattering point that originated it. The D S R equation is the link that connects the scatter point in the migrated output to the sample in the data. A jco filter for 3D data or y/juj for 2D data, where to is the frequency and j — completes the migration process, and a final stack gives the final migrated section. In the E O M method (Bancroft et al., 1998) the same final goal is achieved but the traveltime equation is expressed in terms of different parameters. See Figure 7.2 to (7.1) Chapter 7. Radon transform applied for migration 177 Figure 7.2: Relation between two way traveltime equation and the equivalent offset traveltime equation. facilitate the following explanation. The source (S) and receiver (R) are assumed to be co-located at the point E , not necessarily close to the midpoint (MP). The distance along the surface between the point E and the scattering point is given by a new parameter (he) named the equivalent offset. This offset is such that the original traveltime (T) remains unchanged from the original geometry, but contains the dependency on x, h, t, and V, required by the kinematics of the process. In other words, the complexity of the D S R is transfered to the computation of one new parameter, and the traveltime can be expressed as a single square root equation. In the next section we will see that, after this re-parameterization, the C S P gather offers a suitable space for hyperbolic Radon filtering. Chapter 7. Radon transform applied for migration 178 By imposing the equality between the double square root (DSR) and single square root (SSR) equations (Figure 7.2), we can obtain a mathematical expression for the equivalent offset (Bancroft et al., 1998), In this equation x = x c m p — x c s p . Hence, a sample with coordinates (xcmp,h,t) maps to a C S P gather at the point (xcsp, h e q u i v , f). Note that the time coordinate does not change. Every C S P gather contains the contribution of the whole data space, but only the energy coming from the same scattering point adds constructively along a hyperbolic summation to produce the migrated output. 7.2 T h e E O M - R T m e t h o d The hyperbolic Radon transform (HRT) , presented in Chapter 4, can produce a parsimo-nious decomposition ofthe scattered energy in the CSPs, because this energy aligns along hyperbolic curves. I have already discussed that, by using inversion constraints, aperture extension and filtering of non hyperbolic events can be obtained by the application of the sparse H R T . Therefore, we have two processes that can be linked together: E O M and R T . E O M maps a volume of data to planes where the signal aligns along hyperbolic curves, and R T decomposes hyperbolic curves in a sparse model space, where desired and undesired events can be separated. By applying R T to C S P gathers, it is possible to obtain a number of benefits. a) The Radon transform provides information about velocities that is useful for up-dating the initial velocity information. b) The filtering carried out by the R T may improve the quality of the gathers, leading to a better velocity estimation. (7.2) Chapter 7. Radon transform applied- for migration 179 c) Multiples can be filtered on the CSPs, rather than on C M P s , with some advantages to be discussed later. d) If a P R T instead of a H R T is applied, the R T operator can be combined with the filtering + stack operators (Schonewille, 2000), resulting in a very fast method to remove multiples and migrate data simultaneously. There are two different applications for the R T - E O M method, denoising and focusing. 7.2.1 Denoising The R T can be applied to remove or attenuate coherent and incoherent noise in the C S P gathers, before the final process of NMO+stack+filtering. Two common types of coherent noise in CSPs are horizontal artifacts and multiples. The horizontal artifacts arise in CSPs because of the mapping of horizontal reflectors at constant time. For other than ID structures these events become more complicated, but still usually show some horizontal coherence from trace to trace. In the standard E O M , after the C S P is formed, a NMO+stack + filtering process produces the final output. In this process, the horizontal artifacts are attenuated partially by destructive interference. O n the contrary, in the R T these artifacts can be mapped to areas of zero moveout, making it possible to eliminate them by a mute process in the R T , without affecting the data. The other type of coherent noise, multiples, can also be separated from the primaries by prediction and subtraction process. The CSPs are a better scenario to remove mul-tiples than C M P s , because they map to the CSPs with R M S velocities, facilitating the separation of primaries and low velocity multiples from dipping horizons. The attenuation of incoherent noise brings some difficulties. The energy originated at a single scatter-point will be parsimoniously represented in a space spanned by a set of hyperbolic basis functions. The incoherent noise, however will map all over the R T space, and most of it will appear back in the CSP. Some filtering can be obtained Chapter 7. Radon transform applied for migration 180 Figure 7.3: Difference between mapping to one trace (Kirchoff) and several traces (RT). by using the nullspace of the transformation. This nullspace can be controlled by the range of velocities spanned in the R T space, and by the regularization method applied. For example in subspace methods, part of the incoherent noise can be attenuated or removed by limiting the Krilov space to the vectors associated with the largest singular values. This procedure can be risky, however, because there may be loss of desired signal associated with small singular values, like local A V O effects. Some data filtering can also be achieved by data down-weighting, which allows a softening of the data constraint for outliers, bad traces or other data problems. This data down-weighting is particularly important for the CSPs to avoid fitting artifacts generated during the mapping. Chapter 7. Radon transform applied for migration 181 7.2.2 F o c u s i n g The second application that I have tried is improving focusing by computing velocity corrections. To understand the different nature of the E O M output and the R T - E O M output, we must look at the differences between NMO+stack and R T . In the NMO+stack a single hyperbolic path is scanned at each output value for every output trace. If the velocity model departs from the true migration velocity, the resolution at the output decreases and the quality of the migrated image is diminished, because the NMO-stack process with a wrong velocity produces a loss of focusing in the output. The peak of energy maps either too high (low velocity) or too low (high velocity), but in both cases, the amplitude of the output is smaller than when using the right velocity (Figure 7.3). In the R T the process is different. In the first iteration, the R T operator scans a range of velocities, and the best velocity produces the largest output in the adjoint panel, similar to a semblance spectrum. As the inversion proceeds iteratively, more relevant information is added to the velocity space to properly fit the data, while simultaneously preserving the localization of the information in the velocity space. Therefore a high resolution velocity spectrum can be obtained directly from the Radon space. In the R T used in E O M , the velocity space is constructed not by the actual velocities but by the deviations of the R M S velocities from the initial velocity model. Therefore, inaccuracies in the velocity model can be obtained from the R T space and used to correct the prior velocity information. 7.3 I m p l e m e n t a t i o n The description of the implementation can be separated into two sections, one refering to the E O M method, and the other concerning the Radon transform. I have already described the R T implementation in Chapter 4, and the E O M method is described in Chapter 7. Radon transform applied for migration 182 Large Loop -* Loop lor CSP* -* Loop II r traces. Mapping: d(t,x,h)--dc:p(tx;sp,he) Smal Loop 1 Loop lor CSPs Sparse Hyperbolic Radon T----To Amplitude scaling; Itilorpolalion Filtering Stack Figure 7.4: The E O M + Radon transform. detail in Bancroft, et al. (1998). Because of this, I will skip many details and give only a general description. 7.3.1 Equivalent Offset M i g r a t i o n and C S P s I have implemented the E O M method in two loops, one external (large loop) for the mapping from traces to CSPs, and another internal (small loop) for the mapping from CSPs to output traces (Figure 7.4). In the external loop, the C M P s are read in sequence and mapped into the current CSP. The R T is implemented inside the inner loop, act-ing sequentially on each CSP. Finally the Kirchhoff-NMO process ( N M O + amplitude weighting + interpolation + filtering + stack) gives the final output. To reduce the access to the hard disk, I compute several CSPs at the same time by creating a small volume Chapter 7. Radon transform applied for migration 183 of CSPs (3D subarray) that slides along the larger volume of data. The process of mapping from C M P s to CSPs requires, in theory, the computation of the equivalent offset at every sample in (h,x,t). In practice the equivalent offset is computed in a faster way, because the equivalent offset axis is discrete and does not change very fast from one offset bin to the next. Using this property, a faster approach given in Bancroft et al. (1998), consists of computing transition times from one offset bin to another. There are only a few of these transition times per trace. Given a trace at some coordinates (x, h), the time of the first possible arrival is t0 = 2max(x, h)/V ). We can compute the corresponding equivalent offset for this time and assign it to the offset bin at hei. Then we can precompute the transition time to the next offset bin hei+i. Using this approach, a set of time and equivalent offsets is precomputed for every trace, and the whole trace t(x, h, t) is mapped to a set of equivalent offsets. This mapping is very fast and the process must be repeated for every trace and for every CSP. As we see from equation (7.2), a velocity function is required to compute the equivalent offset. W i t h this purpose and following Bancroft et al. (1998), I first apply a velocity analysis to C M P s at regular intervals along the midpoint axis. Then I use this initial velocity model to compute some CSPs at regular intervals. Finally I again apply a velocity analysis on these CSPs. 7.3.2 R a d o n T r a n s f o r m The Radon transform on CSPs can be performed as a hyperbolic time domain R T with time variant kernels and weighting functions, or as a parabolic R T in the frequency domain. The advantages of the time domain H R T over the frequency domain P R T have been discussed in Chapter 4. I summarize them here: 1) A hyperbolic R T localizes events better than a parabolic approach, in particular for deep water marine data with multiples. Chapter 7. Radon transform applied for migration 184 2) The model weights are easier to calculate because only one hyperparameter is required. As a consequence, they are linked together and produce more stable solutions. 3) The data weights allow one to deal with time variant artifacts and outliers. O n the other hand, frequency domain R T is faster, and frequency domain operators pre-serve much better the waveform of the events. However, because the similarity between simple NMO-stack process and the hyperbolic Radon transform makes easier to discuss the connection between them, I will assume the H R T for the following discussion. The implementation is similar to the already explained (Chapter 4). I will just men-tion briefly the main characteristics. • The kernel is stored as a sparse matrix in the R A M to decrease the computing time. • The external iterations required to enforce sparseness can be avoided by using the results from the previous CSPs to compute the model weights. • The Radon C S P space is kept small by traces with velocity dependence along r . • The velocity space is formed by deviations from the initial velocity model. Concerning the processing time, the Radon E O M is more expensive than the E O M . The total time required is similar to the time it takes to apply R T multiple removal and E O M migration. However, if a P R T is applied, a very fast approach is possible by combining the R T operator with the filtering + stack operators (Schonewille, 2000). Velocity Corrections The R T space provides velocity information to correct the original velocity model, but this information presents local variations and some smoothing is required. W i t h this Chapter 7. Radon transform applied for migration 185 Figure 7.5: Computing velocity corrections on the output Radon space. purpose in mind, I apply a weighted summation inside a moving window on the Radon space (Figure 7.5). The window extends over a few points along the time and Radon axes and the weights correspond to the amplitude of every element (normalized by the total summation). The window moves from the top to the bottom, providing a r dependent correction for the original velocity model. This correction must be kept small because of the nonlinearity of the process. In fact the procedure fails when the errors are too large. Because of this, I limit the velocity span at no more than 10 — 20% of the mean velocity. The corrections are applied iteratively a few times until the computed corrections tend to zero. Because changes in a particular layer affects all deeper layers, velocity adjustments must proceed from top to bottom, in small corrections and iteratively. In this fashion, it is possible to perform the velocity Chapter 7. Radon transform applied for migration 186 correction automatically at every CSP. 7.4 Examples Some insight into the possibilities of a combined R a d o n - E O M method can be obtained by applying the method to a synthetic data set. The synthetic data set is generated using finite difference modeling (program "sufdmod2", from Seismic Unix (SU)). About 40 shot gathers, each one with 200 receivers, are simulated according to the model of Figure 7.6(a). The parameters for this model are provided in S U and can be obtained by using the function "unif2". The velocity model could be obtained by semblance analysis but for this test I used the R M S velocities calculated from the true model parameters. This allows me to isolate velocity errors and check for the corrections in the velocity model performed by the algorithm. The synthetic data set is basically free of noise but it has many numerical artifacts resulting from the finite difference modeling. I applied a t2 gain and removed the direct wave from all the shot gathers. There are free surface multiples but they are not very strong because I used absorbing boundary conditions. Figure 7.6(b) is the result of E O M using correct velocities. It is possible to see some multiples that have gone through the migration and stacking processes. Figure 7.6(c) is the E O M result after I introduced some velocity errors between 10-20 % at selected places, on the left, center and right of the model. The effects of the wrong velocities are evident on the migrated section. Figure 7.6 (d) is the result of Radon E O M . Notice how the distortions on the section have been corrected almost completely. In the second example, using the same synthetic data, I tried to extract coherent signal from a very high level of random noise. Figure 7.7(a) shows the E O M result for data highly contaminated by Gaussian noise, with a signal to noise ratio equal to 1. Figure 7.7(b) shows the E O M result after denoising the C S P gathers by using R T . Part of Chapter 7. Radon transform applied for migration 187 the random noise has been filtered out by stopping the C G algorithm before completion. In the third example we will see the improvements in velocity analysis by using the R T . Bancroft et al. (1998) show that CSPs gathers give more precise velocity analysis than C M P s , first because they have larger fold and secondly because they give R M S velocities. By using R T , I expect to improve velocity analysis even more. We can test if this true by comparing the semblance analysis for C M P s , CSPs, and Radon-CSPs, located at three different positions along the model, the left hand side ( C M P 300), middle ( C M P 500), and right hand side ( C M P 900) (Figure 7.8). The solid line represents R M S velocities estimated from the real model parameters. The C M P gathers at 300 and 500m have enough fold to produce good velocity estimates. We can see that C M P s , CSPs and Radon-CSPs give similar information. On the contrary, the gather at 900m has very low fold and no information can be obtained from the C M P by semblance analysis. The C S P gathers, on the other hand, give very precise information. This is because CSPs obtain information not only from the C M P at the particular location, but also from neighboring C M P s . The previous example has shown that velocity analysis of C S P gathers gives better velocity information than analysis of C M P gathers for low fold data, but it has not shown any important advantage for analysis of Radon-CSP gathers, other than perhaps a slightly better resolution. However, the results become more interesting when substantial noise is added to the data (Figure 7.9). In this situation, semblance analysis on the C M P s is not useful at all (Figure 7.9 top row), whereas CSPs give better results (Figure 7.9, center row). Radon E O M (bottom row) improves the velocity analysis for the C S P by enhancing considerably the resolution of the velocity spectrum. Very interesting is the fact that the best result is for the low fold C M P at 900m. This reflects the fact that information is obtained from a whole volume, and the noise stacks out very well. The next example deals with the same sub-salt data set from the Gulf of Mexico, Chapter 7. Radon transform applied for migration 188 used previously. This data set is not appropriate for time migrations because of the large velocity contrast on the top of the salt dome. However, the strong surface multiples in these data are interesting to test multiple elimination in the CSPs. The primaries are very weak and most of the high energy reflexions that appear after 3 seconds are water bottom multiples. I used the R T - E O M method for multiple-attenuation + migration, with some success. Figure (7.10) shows the results for (a) simple stack, (b) migrated with E O M , (c) migrated with R a d o n - E O M with subtraction of multiples, (d) the same but Radon transform is computed in the frequency domain. The R a d o n - E O M migration has been partially successful in attenuating multiples at about 4 s. The same method with the P R T (in the frequency domain) has been less successful to separate the multiples, mainly because the parabolic approximation in the CSPs did not work so well. However, the waveform is better preserved because frequency domain operators behave better, in this aspect, than time domain operators. Figure (7.11) shows the intermediate process of muting horizontal artifacts and sub-tracting multiples. To show more clearly the improvements I have created a poor sampling case by eliminating three of every four shots in the data set, such that the number of traces per C M P decreases from 92 to 23. The offset sampling interval increases from 175 feet to 700 feet. Figure (7.11a) shows the C M P gather. Note that, because of the poor sampling, if the R T were applied directly on this gather aliasing would be a serious problem. Figure (7.11b) shows the C S P gather at the same location as the C M P 1000. The sampling now is better, such that the R T can be applied. The solution to the aliasing problem results because of the gathering of more information from other C M P s . How-ever, it arises also because when constructing the CSPs, an average of the energy is taken in a window around the true offset. This average has a similar effect to the anti aliasing filter commonly applied in Kirchhoff algorithms. Several horizontal artifacts appear in the CSPs. To remove these artifacts I use a high velocity area (almost null moveout) in Chapter 7. Radon transform applied for migration 189 the Radon space. The nature of this area is different from the rest of the R T space in the sense that there is no velocity variation with time. The horizontal artifacts do not have such variation. Figure (7.11c) shows the same C S P after muting artifacts, subtracting multiples and aperture extension. The semblance analysis shown at the bottom of every gather reveals that the C M P and C S P have similar velocity information, with strong surface multiples hiding the primary velocity trend. The C S P after Radon reveals much more clearly the primaries, whose velocity trend can be seen in the semblance analysis. A n interesting characteristic of C S P gathers is that they contain global information. A C S P at a given location contains the C M P at that location but can be constructed without it. This characteristic may be irrelevant for very well sampled data, like the Mississippi Canyon data, but may be very important to obtain velocity information in poorly sampled data. For example, imagine an unfortunate loss of data on the reduced version of the Mississippi Canyon data set. Suppose that the C M P 1000 has disappeared mysteriously from the tape. Can we still obtain information on this location from the other C M P s ? Figure 7.12a shows the C S P at location 1000 created without using the C M P 1000. Let us now compare with C M P 1000 and its velocity spectrum that are shown in Figures (7.11)a and d. Figure 7.12b shows the same C S P after multiple subtraction. Figure 7.12c shows the same gather after forward and inverse R T , with muting in the R T space of the horizontal artifacts. The respective velocity spectra are shown in the same figure, as a measure of information of the gathers. As a final comment to the last example, notice that the C S P gather is not the same as the C M P gather. This is the reason why I cannot present a residuals plot. We can not just replace a missing C M P by a CSP, but both gathers provide velocity information. C M P s provide stacking velocities, whereas CSPs give R M S velocities. The example shows that CSPs contain global information obtained from traces with different C M P numbers, and because of that they are more robust to missing information. A C S P becomes a Chapter 7. Radon transform applied for migration 190 C M P when the aperture on the gathering process is set to one midpoint location. 7.5 Discussion The combined R T - E O M process is completely automatic, no user intervention is required. As the C M P s are read sequentially along the line, the CSPs are formed and the hyper-bolic R T models are calculated. Low velocity multiples are removed by modeling and subtraction, the incoherent noise and horizontal artifacts are removed by muting in the R T space. Finally a Kirchhoff-NMO-stack process produces the final output. Some changes may be required for this sequence. These changes concern efficiency, but also can change the results substantially. In particular, I think that applying the H R T in hyperbolic windows must lead to a superior accuracy in the reconstruction. This results from the fact that, when inverting integral operators, global features are adjusted first (large singular values of the kernel) and local details later (small singular values). When the C G algorithm is truncated after a few iterations, as is required when working with large gathers, many of the details of the original data are lost. If the gather is decomposed into small windows, then few iterations are enough to fit details, and results may improve considerably. I mentioned already that the multiples are removed by subtraction, an efficient approach because as much as possible of the original gather is preserved. The horizontal artifacts, however, have been eliminated by filtering (muting) in the R T space. This is a disadvantage because part of the data can be lost in the nullspace of the transformation. These artifacts could be eliminated also by subtraction, but I found this procedure very difficult because modeling them precisely is extremely difficult due to their complicated nature, truncations and variable offsets. One practical problem difficult to solve concerns the amplitude scaling along offsets in the CSPs. This problem is generated by the variation of fold with offset and time. For Chapter 7. Radon transform applied for migration 191 each CSP, I use an array to keep the precise account for the fold at every output sample. After the mapping I divide every sample in the CSPs by its corresponding fold. This procedure helps to produce more regular CSPs, but often it over-weights samples near zero offset. Other empirical procedures that I tried did not give better results. This is a serious problem for the R T because to achieve a precise fit of these artificial A V O effects the R T space becomes more complicated and less sparse. As a result, the separation of primaries and multiples becomes more difficult. Finally, another difficult problem involves the automatic velocity corrections when strong multiples contaminate the data. Multiples pull velocity corrections towards their velocity trend. A possible solution for this problem is to apply first the multiple removal strategy, and then the velocity corrections. 7.6 Conclusions Prestack time migration is a good first estimate of the depth model because of its ro-bustness and speed. In this chapter I have combined the Radon transform with a time migration method, E O M . This combined procedure pursues two main objectives, noise attenuation and velocity corrections. Both coherent and incoherent noise can be filtered. Some artifacts produced on CSPs can be removed as well, without affecting the signal. In the standard E O M method, velocities can be estimated iteratively using CSPs, but in the end they are interpolated to produce a velocity model at every point of the section. Hence, we must always expect departures from the velocities that optimally focus the energy of the migrated section. Radon E O M is an attempt to improve velocity estimation from CSPs by automatically reducing velocity errors. The simple tests carried out in this study suggest that the R T - E O M method could provide a very good estimate of velocities and a robust prestack time migration at the Chapter 7. Radon transform applied for migration 192 expense of a reasonable increase in computing time. The method presented here is somewhat more complex than the conventional E O M method but appears to hold promise for the future. The combined R T - E O M method may produce superior results to a R T -Kirchhoff migration approach in the following cases: 1) When the multiples overlap primaries in the velocity space because of high stacking velocities due to dips. 2) When incomplete data produce heavy artifacts and aliasing in the R T . There are several problems that affect the method, and still more research is required to produce optimal results. Many of these problems are common to other migration methods, and I believe it is a matter of time before they are overcome. Chapter 7. Radon transform applied for migration 193 Figure 7 .6: (a) Model, (b) E O M using correct RMS velocities, (c) E O M using wrong velocities, (d) Radon E O M using wrong velocities. Chapter 7. Radon transform applied for migration CDP(m) CDP(m) 200 400 600 800 200 400 600 (a) (b) Figure 7.7: Signal to noise ratio =1 (a) E O M (b) Radon E O M . Chapter 7. Radon transform applied for migration 195 Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (nVsec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Semblance plot for CMP 300 Semblance plot for CMP 500 Semblance plot for CMP 900 Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Semblance plot for CSP 300 Semblance plot for CSP 500 Semblance plot for CSP 900 Velocity (m/sec) 1500 2000 2500 3000 3500 Of Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Semblance plot for CSP 300 Semblance plot for CSP 500 Semblance plot for CSP 900 Figure 7.8': Semblance analysis on clean gathers located at 300m, 500m and 9000m. first row: C M P gathers, second row CSP gathers, third row, CSP gathers + RT (c) Chapter 7. Radon transform applied for migration 196 Velocity (m/sec) 1000 1500 2000 2500 3000 Velocity (m/sec) 1000 1500 2000 2500 3000 Velocity (m/sec) 1000 1500 2000 2500 3000 Semblance plot for CMP 300 Semblance plot for CMP 500 Semblance plot for CMP 900 Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Semblance plot for CSP 300 Semblance plot for CSP 500 Semblance plot for CSP 900 Velocity (m/sec) 1500 2000 2500 3000 3500 0+ Velocity (m/sec) 1500 2000 2500 3000 3500 Velocity (m/sec) 1500 2000 2500 3000 3500 Semblance plot for CSP 300 Semblance plot for CSP 500 Semblance plot for CSP 900 Figure 7.9: Semblance analysis on noisy gathers located at 300m, 500m and 9000m. first row: C M P gathers, second row CSP gathers, third row, CSP gathers + RT (c) Chapter 7. Radon transform applied for migration 197 cdp cdp 1000 2000 3000 1000 2000 3000 HRT(t ime)-EOM HRT(freq)-EOM Figure 7.10: (a) Stack, (b) E O M , (c) R A D O N - E O M , (d) Radon (in frequency) - E O M Figure 7.11: (a) C M P gather at location 1000. (b) CSP gather at location 1000. (c) CSP gather after filtering horizontal artifacts in the RT domain and subtracting multiples, (d-e-f) Semblance analysis from (a-b-c) respectively. Figure 7.12: (a) CSP gather at location 1000 obtained without the C M P gather 1000. (b) The same gather after multiple subtraction, (c) The same gather after multiple subtraction and filtering artifacts, (d-e-f) Semblance analysis from (a-b-c) respectively. C h a p t e r 8 D i s c u s s i o n a n d C o n c l u s i o n Todo pasa y todo queda, pero lo nuestro es pasar, pasar haciendo caminos, caminos sobre la mar. Caminante, son tus huellas el camino y nada mas; caminante, no hay camino, se hace camino al andar. Al andar se hace camino y al volver la vista atras se ve la senda que nunca se ha de volver a pisar. Caminante no hay camino sino estelas en la mar... Antonio Machado I have presented different implementations and applications of the sparse R T . I will not discuss the details here since every chapter contains its own conclusions. Let me, however, summarize briefly the central ideas presented in this thesis. I have, initially, outlined the problems that appear in the R T as a consequence of discrete sampling and limited aperture, and how the sparseness constraint helps to deal with some of them. I have analyzed and proposed different numerical methods, or im-plementations, to compute the R T subject to a sparseness constraint and have discussed the advantages and possibilities of computing the R T in both the time and frequency domains. I have also presented new methods for separating signal and coherent noise in different scenarios which are a result of designing convenient modeling operators and applying the C G methodology proposed by Claerbout (1992). The flexibility that is achieved in the R T when the transform space is designed with nonuniform sampling and/or different operators has allowed the design of non conventional R T algorithms for the transformation of complicated data. The application of this methodology to migra-tion shows some of the possibilities when migration and the R T are combined in a single 200 Chapter 8. Discussion and Conclusion 201 process. In the beginning of this thesis I posed two questions: • Does a sparseness constraint solve the discretization and aperture problems for the R T ? The examples in this thesis have hopefully convinced the reader that, often, it does. • When should a sparse R T be applied? The sparse R T will, in general, produce better results than a simple R T when the data have limited aperture, gaps, or multiples. It is probably not needed when the sampling and aperture are well above the minimum requirements imposed by the data and applied operators. Another important question is, how will the R T evolve in the future? I believe that the sparse R T will replace the standard R T completely since, not only is it a generalization of the standard R T , but is also superior in performance. To achieve this, it is important that the algorithms be efficient, robust and have a smooth transition between non sparse and sparse solutions. It is likely, in my opinion, that the applications of the R T as a filtering or interpola-tor tool will disappear at some point in the future from the seismic processing garden. Eventually, multiples will be useful for migration and inversion, or they will be removed with wave equation methods. Interpolation, perhaps, will be performed with more so-phisticated seismic operators based on the wave equation, or not be required at all, as migration algorithms might not require regular sampled data any longer. However, in the near future, there will still many useful R T contributions to seismic processing. Most importantly, there are still several unresolved problems: • A V O preservation. When the algorithms are unable to yield small residuals, poor misfit leads to the destruction of A V O effects. The introduction of an extra parame-ter in the inversion may improve the data fit and open a broad range of possibilities. Chapter 8. Discussion and Conclusion 202 • Frequency domain algorithms are fast but produce several artifacts when the sparse-ness constraint is enforced. Imposing a time domain constraint but still working in the frequency domain is a difficult challenge. • Time domain algorithms are still too slow. Some possible solutions to this problem could involve multi-grid methods, working with hyperbolic windows, or simply data with regular geometries where the inverse operators could be pre computed. In this sense, the R T - E O M algorithm presented in this thesis is appropriate for this kind of scheme. • A R T capable of parsimoniously transforming events with non hyperbolic moveouts is required in many geological areas, in particular in the presence of anisotropy. Also it could be a necessity for deep water marine data and faulted structure. • Better and more efficient algorithms are required. Perhaps more importantly, effi-cient methods for finding hyperparameters or better stopping criteria, for both the external and internal iterations (when IRLS is used) should be very useful. I hope that this thesis has shown some of the problems of the R T , some possible solutions to them and has opened new routes for research. Furthermore, I hope it has shown some of the beautiful landmarks along the path running between signal processing and inversion. References Alliney, S. and Ruzinsky, S. A . , 1994. A n algorithm for the minimization of mixed l\ and I2 norms with application to Bayesian estimation, I E E E Trans. Signal Processing, 42, 618-627. Bancroft, J . , Geiger H . and Margrave, G . , 1998, Equivalent Offset Migration: Geophysics, 63, N O . 6, P. 2042-2053. Bancroft, J . , 1999, A practical understanding of Pre- and Poststack Migrations: S E G , Course Notes Series, N O . 9. Beylkin Gregory, 1987, Discrete Radon Transform: I E E E transactions on acoustics, speech, and signal processing, 35, 162-172. Bickel S H , 2000, Focusing aspects of the hyperbolic Radon tranform: Geophysics, 65, 2, 652-655. Cary P., 1998, The simplest discrete Radon transform: 68th A n n . Internat. M t g Soc. Explo.Geophys., Expanded Abstracts, 1999-2002. Chen S. S., Donoho D . L . , Saunders M . A . 1998 Atomic Decomposition by Basis Pursuit, S I A M Journal on Scientific Computing. Volume 20, Number 1, pp. 33-61 Chong E . and Zahk, S., 1996, A n Introduction to Optimization: Wiley-Interscience Series in Discrete Mathematics and Optimization. Claerbout, J . F . , 1992, Earth sounding analysis, Processing versus inversion: Blackwell Scientific Publications, Inc. 203 REFERENCES 204 Claerbout J . F . , 1985, Imaging the Earth: Blackwell Scientific publications. Dunne, J . and Beresford, G . , 1997, Processing through to stack in the r — p domain, Exploration Geophysics, 28, 341 - 347. Foster, D . J . , and Mosher, C . C , 1992, Suppression of multiple reflections using the Radon transform: Geophysics, 57, 386-395. Haber, E . , 1997, Numerical strategies for the solution of inverse problems. Ph .D . thesis, University of British Columbia. Hale, D . . , 1991, A Non aliased method for time domain D M O : Geophysics, June 1991. Hampson, D . , 1986, Inverse velocity stacking for multiple elimination: Journal of the Canadian S E G , 22, 44-55. Hanke, M . , and Hansen, P. 1993, Regularization methods for large scale problems: Sur-vey. Math. Ind., 3, 253-315. Hansen, P. 1998, Rank-Defficient and Discrete Ill-Posed Problems: Numerical aspects of Linear Inversion. S I A M monographs on mathematical modeling and computation. Hansen P. C . 1999, Numerical aspects of Deconvolution: http://www.imm.dtu.dk/ pch Herrmann P., Mojesky T , and Magesan M . , 1999, Amplitude prserving Radon demulti-ple: beyond sampling and aperture limitations: C S E G 1999 National Convention, Abstracts, 73-74. Huber, P. J . , 1981. Robust Statistics: John Wiley & Sons, New York. Hugonnet, P. and Canadas, G . , 1995, Aliasing in the parabolic Radon transform: 65th A n n . Internat. Mtg. , Eur. Assoc. Explo. Geophys., Expanded Abstracts, 1366-1369. REFERENCES 205 Kabir M . M . Nurul, and Verschuur D. J . . , 1995, Restoration of missing offsets by parabolic Radon transform: Geophys. Prosp., 43, 347-368. Karmarkar, N . 1984, A new polynomial-time algorithm for linear programming. Combi-natorica, 4(4), 373-395. Rostov, C . , 1990, Toeplitz structure in slant-stack inversion: 60th A n n . Internat. Mtg Soc. Explo. Geophys., Expanded Abstracts, 1618-1621. Lumley, D . , 1997, Estimating a pseudounitary operator for velocity stack inversion: Stan-ford Exploration Project Report, 82, 63-78. Lumley, D . , Nichols, D . and Rekdal, T . , 1995, Amplitude-Preserved Multiple Suppres-sion: 65th Ann. Internat. Mtg Soc. Explo. Geophys., Expanded Abstracts, 1460-1463. Maeland Einar, 1998, Focusing aspects of the parabolic Radon transform: Geophysics, 63, 5, 1708-1715. Marfurt, K . and Nemeth, T . , 2001, (paper in review process), Antialias discrete Radon transforms: Geophysics. Marfurt, K . J . , Schneider, R. J . , and Mueller, M . C , 1996, Pitfalls of using conventional and discrete Radon transforms on poorly sampled data: Geophysics, 61 , 5, 1467-1482. Nichols, D . , 1997, A simple example of a null space and how to modify it: Stanford Exploration Project Report, 82, 185-192. Olson T i m and DeStefano Joe, 1994, Wavelet Localization of the Radon Transform: I E E E , Transactions on signal processing, 42, 8, 2055-2067. REFERENCES 206 Radon J . , 1917, Uber die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten, Berichte Sachsische Acadamie der Wissenschaften, Leipzig, Math.-Phys. KI. , vol. 69, pp. 262-267. Saad, Y . , 1996, Iterative methods for sparse linear systems: P W S publishing Company. Sacchi, M . and Porsani, M . , 1999, Fast high resolution parabolic R T : 69th Ann. Internat. M t g Soc. Explo. Geophys., Expanded Abstracts, 1477-1480. Sacchi, M . and Ulrych, T . , 1997, Compression and Inversion of velocity-stack operators: C S E G 1997 National Convention, Abstracts, 77-79. Sacchi M . , 1996, A bidiagonalization procedure for the inversion of time-variant velocity stack operator: C D S S T report, 73-92. Sacchi, M . , 1995, Ph .D. thesis. Aperture compensated Radon and Fourier transforms, The University of British Columbia. Sacchi, M . and Ulrych, T . , 1995, High-resolution velocity gathers and offset space recon-struction: Geophysics, 60, 4, 1169-1177. Sacchi M . and Ulrych T . , 1995a, Improving Resolution of Radon operators using a model re-weighted least squares procedure: Journal of seismic exploration, 4, 315-328. Scales, J . , Gersztenkorn, A . and Treitel, S., 1988, Fast lp solution of large, sparse, linear systems: Application to seismic travel time tomography, Journal of Computational Physics, 75, 314-333. Sheriff, Robert E . , 1991, Encyclopedic dictionary of exploration geophysics , Encyclopedic dictionary of exploration geophysics: Soc. Of Expl . Geophys. REFERENCES 207 Schonewille, M . and Duijndam, A . , 2001, Parabolic Radon transform, sampling and efficiency, 66, 2, 667-678. Schonewille, M . , 2000, Ph .D. thesis. Fourier reconstruction of irregularly sampled seismic data. Delft Unversity of Technology. Shewchuk, J . , 1994. A n introduction to the Conjugate Gradient Method Without the A g -onizing Pain, School of Computer Science. Carnegie Mellon University. Pittsburgh. F T P : / / warp. cs. emu. edu Spitz, S., 1991, Seismic trace interpolation in the F - X domain: Geophysics, Soc. of Expl . Geophys., 56, 785-794. Stewart, R., 1991, Exploration Seismic Tomography: Fundamentals: S E G , Course Notes Series, N O . 3, S. N . Domenico, Editor. Strang, G . , 1986, Introduction to applied mathematics: Wellesley-Cambridge Press. Tarantola, A . , 1987, Inverse Problem Theory: Elsevier Publishing. Tatham, R., 1989, Tau-p Filtering, in Stoffa, P. L . , McCowan, D . W . , Tatham, R. H . , Diebold, J . , Tygel, M . , Hubral, P., Wenzel, F . , Ziolkowski, A . , Fokkema, J . T . , Koster, K . , Confurius, A . , Boom, R. V . , Tau-p - A plane Wave Approach to the Analysis of Seismic Data, Modern approaches in geophysics, Paul Stoffa Editor, 35-70: Kluwer Academic Publishers. Thorson, R. and Claerbout, J . , 1985, Velocity-stack and slant-stack stochastic inversion: Geophysics, 50, 4, 2727-2741. Trad D. , Sacchi M . and Ulrych T . , 2001a, A hybrid linear-hyperbolic Radon transform. Journal of Seismic Exploration, 9: 303-318. REFERENCES 208 Trad D . , Ulrych T . , and Sacchi M . , 2001b, Accurate interpolation with high resolution time variant Radon transforms, accepted for publication in Geophysics. Trad D . and Ulrych T . , 2000, Improving Equivalent Offset Migration by Radon Tech-niques, 2000: 70th Ann. Internat. Mtg Soc. Explo. Geophys., Expanded Abstracts, 818-821. Trad, D . and Ulrych, T . , 1999, Beyond aliasing in the Radon transform: SBGf, Rio'99 Conference. Trad D. , Ulrych T . , and Sacchi M . , 1999, The Radon transform, sampling aliasing and computational aspects: C S E G 1999 National Convention, Abstracts, 75-77. Treitel, S., Gutowski, P. R. and Wagner, D . E . , 1982, Plane-wave decomposition of seismograms: Geophysics, 47, 1375-1401. Turner, G . , 1990, Aliasing in the rp transform and the removal of spatial aliased coherent noise: Geophysics, 55, 1496-1503. Ulrych, T . J . , Sacchi, M . D . and Woodbury, A . , 2001, A Bayes tour of inversion: A tutorial: Geophysics, Soc. of Expl . Geophys., 66, 55-69. Yilmaz, O. , 1989, Velocity-stack processing: Geophys. Prosp., 37, 357-382. Yilmaz, O . , 1987, Seismic data processing: Investigations in Geophysics, S E G , Vol 2. Wisecup, D . , 1998, Unambiguous signal recovery above the Nyquist using random-sample-interval imaging: Geophysics, 63, No 2, 763-771. Wiggins, J . W . and Abma, R., 1990, Using the parabolic Radon transform as a moveout filter: 52th A n n . Internat. Mtg. , Eur. Assoc. Explo. Geophys., Expanded Abstracts, 69-70. REFERENCES 209 Zhou, B . and Greenhalgh, S.A., 1994, Linear and parabolic r-p transforms revisited: Geophysics, 59, 1133-1141. Zhu, J . and Lines, L . R., 1997, Implicit interpolation in reverse-time migration: Geo-physics, Soc. of Expl . Geophys., 62, 906-917. Appendix A Table of Symbols S Y M B O L name or description comments R T Radon transform L R T linear Radon transform P R T parabolic Radon transform H R T hyperbolic Radon transform E R T elliptical Radon transform P H R T pseudo hyperbolic Radon transform C M P common midpoint gather the set of traces with the same midpoint between source and detector C D P common depth point gather similar as C M P but the common midpoint is at depth (reflection point) CSP common scattering point gather the set of traces whose energy is produced by scattering at the same point E O M equivalent offset migration N M O normal moveout the variation of the reflection arrival time due to offset distance D M O dip moveout the variation in the reflection arrival time due to dip of the reflector A V O amplitude variation with offset the variation of the amplitude of a seismic reflection with source-geophone distance 210 Appendix A. Table of Symbols 211 S Y M B O L name or description comments DSR double square root equation traveltime as a function of two spatial coordinates (midpoint and offset) SSR single square root equation traveltime as a function of one spatial coordinate (equivalent offset for E O M ) RMS velocity root mean square velocity G R ground roll surface wave energy that travels along or near the surface P slowness slope parameter in L R T Q Radon parameter curvature parameter in P R T and H R T h offset distance between source and receiver np or nq number of p or q traces nh number of offsets he equivalent offset offset parameter used in E O M •Zref reference depth T two way travel time at zero offset V velocity in general it represents RMS velocity U angular frequency (radian/seconds) f temporal frequency in Hertz (Hz) t time in seconds nt number of time samples h-x wavenumber (cycles/meter) spatial frequency 1Some of these definitions are taken from Sheriff (1991). Appendix A. Table of Symbols 212 S Y M B O L name or description comments II • \\m £ m norm I II" n power of £ m norm A, S, 7 , e hyperparameters CT standard deviation C T 2 variance e angle of incidence angle between a ray and the perpendicular to a particular reflector d data time-offset domain m model data after transformation L forward RT operator operator that maps model to data L* adjoint R T operator operator that maps data to model w m model weight matrix data weight matrix model covariance matrix c d data covariance matrix A " hermitian of a matrix A if A is an complex operator then it represents also A* (adjoint of A) A r transpose of a matrix A if A is an real operator then it represents also A* (adjoint of A ) D diagonal matrix R circulant matrix a matrix whose columns are shifted versions of a vector I identity matrix SPD symmetric positive definite Appendix A. Table of Symbols 213 SYMBOL name or description comments pdf probability density function P(a|b) probability of a given b T or FT Fourier transform DFT discrete Fourier transform FFT fast Fourier transform CG conjugate gradient CGLS conjugate gradient least squares PCGNR preconditioned CG for the normal residuals CGFFT CG with FFT matrix vector multiplication peformed by FFT WTCGLS weighted CGLS GCV cross validation function stopping criteria for iterative algorithms SVD singular value decomposition SVD truncated singular value decomposition GSVD general singular value decomposition IRLS iterative reweighted least squares TFD time-frequency domain time domain algorithm with frequency domain operators
Thesis/Dissertation
eng
Oceanography
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
Implementations and applications of the sparse Radon transform
Text