Estimation of Surface-free Data by Curvelet-domainMatched Filtering and Sparse InversionbyMufeed H. AlMatarB.S.c. Computer Engineering, King Fahad University of Petroleum and Minerals,Dhahran, 2005A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE STUDIES(Geophysics)The University Of British Columbia(Vancouver)April 2011c Mufeed H. AlMatar, 2011AbstractA recent robust multiple-elimination technique, based on the underlying principlethat relates primary impulse response to total upgoing wavefield, tries to changethe paradigm that sees surface-related multiples as noise that needs to be removedfrom the data prior to imaging. This technique, estimation of primaries by sparseinversion (EPSI), van Groenestijn and Verschuur; Lin and Herrmann, proposes aninversion procedure during which the source function and surface-free impulse re-sponse are directly calculated from the upgoing wavefield using an alternating op-timization procedure.EPSI hinges on a delicate interplay between surface-related multiples and pri-maries. Finite aperture and other imperfections may violate this relationship. Inthis thesis, we investigate how to make EPSI more robust by incorporating curvelet-domain matching in its formulation. Compared to surface-related multiple removal(SRME), where curvelet-domain matching was used successfully, incorporatingthis step has the additional advantage that matches multiples to multiples ratherthan predicated multiples to total data as in SRME.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Theme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.4 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . 42 Primary-multiple separation techniques . . . . . . . . . . . . . . . . 82.1 Primary-multiple separation techniques . . . . . . . . . . . . . . 82.2 Surface-related multiples model . . . . . . . . . . . . . . . . . . 92.2.1 Surface-related multiples formulation for normal-incidenceplane wave in a horizontally stratified medium . . . . . . 92.2.2 Data matrix . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.3 Two-dimensional formulation . . . . . . . . . . . . . . . 112.3 Iterative surface-related multiple removal . . . . . . . . . . . . . 142.4 Curvelet-domain matched-filtering . . . . . . . . . . . . . . . . . 15iii2.5 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 Convexified EPSI and curvelet-domain matched-filtering . . . . . . 193.1 Estimation of primaries by sparse inversion . . . . . . . . . . . . 193.1.1 Convexified EPSI . . . . . . . . . . . . . . . . . . . . . . 203.1.2 Bi-convex optimization . . . . . . . . . . . . . . . . . . . 213.2 EPSI with curvelet-domain matched-filtering . . . . . . . . . . . 233.2.1 Frequency-domain monochromatic matching . . . . . . . 233.2.2 Time-domain matching . . . . . . . . . . . . . . . . . . . 244 Observations and results . . . . . . . . . . . . . . . . . . . . . . . . 264.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.1.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . 264.1.2 Results evaluation . . . . . . . . . . . . . . . . . . . . . 264.2 Factors affecting our results . . . . . . . . . . . . . . . . . . . . . 274.2.1 The smoothing and positivity of the curvelet-domain diag-onal scaling . . . . . . . . . . . . . . . . . . . . . . . . . 274.2.2 When and how frequent should we do the matching? . . . 324.2.3 Optimized result . . . . . . . . . . . . . . . . . . . . . . 415 Conclusion and future works . . . . . . . . . . . . . . . . . . . . . . 435.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.2 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Appendix A Implementations . . . . . . . . . . . . . . . . . . . . . . . 48ivList of TablesTable 4.1 Synthetic time-domain finite difference acoustic forward mod-eling data parameters. . . . . . . . . . . . . . . . . . . . . . . 27vList of FiguresFigure 1.1 An illustration of land seismic survey. A seismic source sendsenergy into the earth. Some of this energy is reflected back tothe surface at rocks layers discontinuities. This energy is thenrecorded by geophones at the surface. . . . . . . . . . . . . . 2Figure 1.2 An illustration that shows a primary and two types of multi-ples. The solid blue line has only one upward reflection, so itis a primary reflection. The green dashed line has more thanone upward reflection and it has its shallowest downward re-flection at the surface, hence, it is a surface-related multiple.Finally, the internal-multiple, the red dashed line, has its shal-lowest reflection at an internal rock layer and it has more thanone upward reflection. . . . . . . . . . . . . . . . . . . . . . 3Figure 1.3 Mesh view of a 2D curvelet. Curvelets are oscillatory in onedirection, but smooth in the other. They resemble a windowedsinusoid in the oscillatory direction and a very smooth Gaus-sian window in the other direction. . . . . . . . . . . . . . . . 5Figure 1.4 Spatial (left) and frequency (right) representations of six curveletswith different angles and scales. Notice how each curvelet islocalized in the spatial domain with a rapid decaying amplitudeand how it has a compact support in the frequency domain.Notice also the 90 rotation between any two correspondencerepresentations in both domains of each curvelet. . . . . . . . 6viFigure 1.5 a) Synthetic shot gather b) Curvelet decomposition of the shotgather into different scales and angles. Five scales are shownhere. Notice how the angles double every other scale and howthe different parts of the shot gather are decomposed at thevarious scales and angles. The coarsest (center) scale showsthe DC and low frequencies. The second scale has 4 anglesand the third and the fourth has 8. The fifth (finest) scale has16 angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6Figure 1.6 Curvelet correlation with wavefronts. A large coefficient isyield when the curvelet and the waterfront have locally thesame direction and frequency contents. Otherwise, a small co-efficient is yielded. . . . . . . . . . . . . . . . . . . . . . . . 7Figure 2.1 Two reflection events: a primary (solid blue line) and a surface-related multiple (dashed red line). Both events have similar ar-rival time but have traveled different paths. The primary eventhas traveled deeper into the subsurface and has seen highervelocity than the multiple event. Therefore, it arrived with asmaller angle at the receiver. . . . . . . . . . . . . . . . . . . 9Figure 2.2 Feedback diagram for surface-related multiples generation. First,surface-free data is generated by convolving the source signa-ture s+ with surface-free impulse response g. Then, surface-free data feedbacks into the earth to generate higher order mul-tiples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11Figure 2.3 a) An illustration that shows the generation of surface-relatedmultiples. b) Source signature s(t). c) Surface-free impulseresponse g(t). d) Surface-free data p0(t). e) First-order multi-ples m1(t). f) Second-order multiples m2(t). g) Total responsep(t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 2.4 After convolving the source gather from total data P with thecommon receivers gather from surface-free impulse responseG for each point in xk, the convolution result is summed up toproduce the predicted surface-related multiple trace at xr. . . . 13viiFigure 2.5 A common-offset section from North Sea field dataset withresults from two different multiples-elimination methods. a)Total data plotted with automatic-gain control. b) Primary es-timation from one term SRME. c) Scaled Bayesian iterativethreshold primaries estimate. Adapted from Herrmann et al.(2008b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Figure 4.1 Zero-offset EPSI synthetic data results with no matching. a)Total upgoing wavefield section ( ˆP). b) Primaries estimation( ˆP0). c) primaries estimation minus the total data ( ˆP ˆP0 = ˆG ˆP), i.e. the multiples estimation d) Final estimated sourcewavelet ( ˆQ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Figure 4.2 Mosaic and 2D plots for the curvelet-domain matched-filteringscaling for one shot (offset = 1260m) for l = 0:46 and h =10;35 and 100. The l and h used to generate each subfigureare stated underneath it. . . . . . . . . . . . . . . . . . . . . 30Figure 4.3 Zero-offset synthetic data results for l = 0:46 and h = 10, 35and 100. The l and h used to generate each subfigure arestated underneath it. a-c) zero-offset primaries estimates. d-f)zero-offset multiples estimates. . . . . . . . . . . . . . . . . . 31Figure 4.4 Mosaic and 2D plots for the curvelet-domain matched-filteringscaling for one shot (offset = 1260m) for h =35 and l = 0.046,0.46 and 46.4. The l and h used to generate each subfigureare stated underneath it. . . . . . . . . . . . . . . . . . . . . 33Figure 4.5 Zero-offset synthetic data results for h = 35 and l = 0.046,0.46 and 46.4. The l and h used to generate each subfigureare stated underneath it. a-c) zero-offset primaries estimates.d-f) zero-offset multiples estimates. . . . . . . . . . . . . . . 34Figure 4.6 a-d) Primaries estimates at different iterations. e-h) Multiplesestimates calculated by subtracting the primaries estimate forthe total data ˆP ˆG ˆQ. i-l) Multiples estimates calculated byconvolving the Green’s function with the total upgoing wave-field ˆG ˆP. . . . . . . . . . . . . . . . . . . . . . . . . . . . 36viiiFigure 4.7 Zero-offset results of EPSI with only one matching step at thefirst alternating optimization loop. a) Zero-offset primaries es-timate. b) Zero-offset multiple estimate calculated by subtract-ing the primaries estimate form the total upgoing wavefield. . 37Figure 4.8 Zero-offset results of EPSI with six matching steps starting af-ter the second alternating optimization loop (even loops only).a) Zero-offset primaries estimate. b) Zero-offset multiple esti-mate calculated by subtracting the primaries estimate form thetotal upgoing wavefield. . . . . . . . . . . . . . . . . . . . . 38Figure 4.9 Zero-offset results of EPSI with only one matching step at thetwenty first alternating optimization loop. a) Zero-offset pri-maries estimate. b) Zero-offset multiple estimate calculatedby subtracting the primaries estimate form the total upgoingwavefield. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 4.10 Zero-offset results of EPSI with only one matching step at thetwelve alternating optimization loop. a) Zero-offset primariesestimate. b) Zero-offset multiple estimate calculated by sub-tracting the primaries estimate form the total upgoing wavefield. 40Figure 4.11 Zero-offset EPSI synthetic data result with curvelet-domainmatched filtering incorporated into loop 12. The number of it-erations h = 60 and the trade off parameter l = 0:46. a) Totalupgoing wavefield section ( ˆP). b) Primaries estimation ( ˆP0).c) primaries estimation minus the total data ( ˆP ˆP0 = ˆG ˆP),i.e. the multiples estimation d) Final estimated source wavelet( ˆQ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42ixAcknowledgmentsI would like to thank my supervisor Felix Herrmann for his continuous support andencouragement, and creative ideas without which this thesis would not be possible.I would like also to thank all slim members for their support in filling the gapsin my Geophysical background.I am very fortunate to have Michael Bostock and Eldad Haber as part of mysupervisory committee. I thank Michael for his two insightful courses and Eldadfor the few lectures I attending during his course last fall.I would like also to thank my management at Saudi Aramco for giving me theopportunity to study at UBC and for their support and funding during my studies.A special thank goes to my wife Sakeel for her encouragement and patientduring my research period especially after the twins came in board.Thank you all.xDedicationTo my parents,my wife Sakeel,my daughters Zahraa and ZainabxiChapter 1IntroductionIn early nineteenth century, searching for oil was very simple, easy and straight-forward. People were looking for oil seepages at the surface near oil springs. Backthen, the demand for oil was limited and oil was mainly used for medical pur-poses (Gadallah and Fisher, 2005). But the situation has changed since then, manyobjects around us are made from oil or made by machines that use energy fromoil and gas. The paper used to print this thesis was manufactured using oil energy.The ink itself is an oil product. The demand for oil has changed dramatically andso do procedures to search for oil.The objective of the seismic survey is to reveal the subsurface geology. It is avery costly and complicated procedure. On land, dynamite (impulsive source) orVibroseis (trucks that propagate energy into the earth over an extended period oftime), sources send energy into the subsurface. Some of this energy is reflected atthe boundaries between rocks layers back to the surface where it is recorded by anarray of geophones. The same procedure applies at sea but the seismic source isnow an air gun and the receivers are hydrophones. Figure 1.1 illustrates the seismicsurvey on land.Most seismic imaging algorithms assume that the reflected energy recorded bygeophones at surface has reflected only once in the subsurface (Verschuur, 2002),as shown by the solid line in Figure 1.2. However, in practice, this assumption isnot true. Some of the energy we send bounce between rocks layers before beingrecorded by geophones as shown by the dashed lines in Figure 1.2. As a result1Figure 1.1: An illustration of land seismic survey. A seismic source sendsenergy into the earth. Some of this energy is reflected back to the sur-face at rocks layers discontinuities. This energy is then recorded bygeophones at the surface.multiple reflections events occur. These events are usually considered as undesiredevents that need to be removed from the seismic data before processing to avoidconfusion during the interpretation of seismic images. The success of multiple-removal or multiple-suppression methods can make the difference between hittingthe target and making profit or missing the target and losing millions of dollars.One way to categorize multiples is by the interface where they have their shal-lowest downward bounce. The first downward bounce for the green dashed linein Figure 1.2 happens at the surface. Hence it is called a surface-related multiplewhile the first downward bounce for the red dashed line happens at an internalinterface and so it is called an internal-multiple (Verschuur, 2002).1.1 ThemeThe main theme of this thesis is to extend a recent robust multiple-eliminationtechnique that estimates primaries by sparse inversion, EPSI (van Groenestijn andVerschuur, 2009; Lin and Herrmann, 2009). This technique proposes an inversionprocedure during which the source function and surface-free impulse response are2Figure 1.2: An illustration that shows a primary and two types of multiples.The solid blue line has only one upward reflection, so it is a primaryreflection. The green dashed line has more than one upward reflectionand it has its shallowest downward reflection at the surface, hence, it isa surface-related multiple. Finally, the internal-multiple, the red dashedline, has its shallowest reflection at an internal rock layer and it has morethan one upward reflection.directly calculated from the upgoing wavefield using an alternating optimizationprocedure. We propose to incorporate a curvelet-domain matched-filtering stepinto this optimization procedure to mitigate possible adverse effects that stem fromamplitude mismatches that vary smoothly as a function of position and dip alongthe predicted wavefronts of the current estimate for the surface-free data. Theimplementation and the different parameters affecting this matched-filtering stepare investigated.1.2 ObjectivesThis thesis aims to incorporate a matched-filtering operator into the formulation ofEPSI. This operator will be able to absorb amplitude errors of non-ideal reflectionsat the surface, finite-aperture, and other unknown effects.31.3 OutlineChapter 2 introduces primary-multiple separation followed by an explanation ofthe theory underlying the model behind surface-related multiple removal, SRME(Verschuur et al., 1992), and primary estimation by sparse inversion, EPSI (vanGroenestijn and Verschuur, 2009). Finally, we discuss the theory behind curvelet-domain matched-filtering.In Chapter 3, we introduce estimation of primaries by sparse inversion, beforeincorporating a matching step into the convexified EPSI procedure.Chapter 4 studies the factors affecting the outcome of incorporating a matched-filtering step into EPSI. We conclude this chapter by presenting our results.We will conclude our research in chapter 5 by a summary and suggested futurework.The last section of this chapter introduces the curvelet transform and highlightsits importance as a sparsifying domain to the success of the matched-filtering pro-cedure.1.4 Theoretical backgroundThe ability of the transform to handle seismic data with conflicting dips and to min-imize the overlap between the primary and multiples coefficients in the transformdomain is significant to the success of the primary-multiple separation algorithms(Herrmann, 2008b). In this section, we show that the curvelet transform (Candeset al., 2006) addresses these issues by efficiently representing seismic data sparselyand by its ability to handle conflicting dips.Curvelet transform was first introduced in Candes and Donoho (2000). Curvelettransform is a localized, multiscale, multi-directional and redundant transform.Curvelets are oscillatory in one direction and smooth in the perpendicular direc-tion. They resemble a windowed sinusoid in the oscillatory direction and a verysmooth Gaussian window in the other direction, see Figure 1.3. Each curvelet isspatially localized in the spatial domain and has a compact support in the frequencydomain. In spatial domain, its amplitude decays rapidly outside a certain region.Curvelets obey the parabolic scaling principle with length width2 (Neelamaniet al., 2008).4The curvelet transform computational complexity of data of size N is O(NlogN)which is not an issue. However, the curvelet transform demands memory; it is anovercomplete transform, i.e. the resulting transform coefficients are more than theimage coefficients. It is about 8 times redundant for the 2D curvelet transform and24 times redundant for the 3D case.Figure 1.4 shows a few curvelets of different scales and angles with a rapidamplitude decay in the spatial domain and a compact support in the frequencydomain and Figure 1.5 shows a decomposition of a shot gather at different scalesand angles.Curvelet correlate very well locally with wavefronts of the same direction andfrequency contents, see Figure 1.6. This permits the curvelet transform to rep-resent seismic data sparsely as a small set of significant coefficients which arelocation, scale, and dip dependent. This sparse representation of seismic data helpsin mapping the different seismic events to different coefficients and that will makethe overlap of the primaries and multiples coefficients minimal and hence leads toa successful separation.5Figure 1.3: Mesh view of a 2D curvelet. Curvelets are oscillatory in onedirection, but smooth in the other. They resemble a windowed sinusoidin the oscillatory direction and a very smooth Gaussian window in theother direction.Figure 1.4: Spatial (left) and frequency (right) representations of sixcurvelets with different angles and scales. Notice how each curveletis localized in the spatial domain with a rapid decaying amplitude andhow it has a compact support in the frequency domain. Notice also the90 rotation between any two correspondence representations in bothdomains of each curvelet.6(a) (b)Figure 1.5: a) Synthetic shot gather b) Curvelet decomposition of the shotgather into different scales and angles. Five scales are shown here. No-tice how the angles double every other scale and how the different partsof the shot gather are decomposed at the various scales and angles. Thecoarsest (center) scale shows the DC and low frequencies. The secondscale has 4 angles and the third and the fourth has 8. The fifth (finest)scale has 16 angles.Figure 1.6: Curvelet correlation with wavefronts. A large coefficient is yieldwhen the curvelet and the waterfront have locally the same direction andfrequency contents. Otherwise, a small coefficient is yielded.7Chapter 2Primary-multiple separationtechniquesIn this chapter, we start briefly by categorizing the primary-multiple separationtechniques. Then, we introduce the model behind surface-related multiple removal(SRME) and estimation of primaries via sparse inversion (EPSI). Finally, we layoutthe theory behind the curvelet-domain matched-filtering.2.1 Primary-multiple separation techniquesMultiples removal techniques can be categorized in two main classes: filtering methods based on move-out and dip discrimination. data-driven methods based on periodicity and predictability.Filtering methods rely on the fact that multiple events traveled through differentrock layers. Hence, they traveled at an apparent different velocity or reflected atdifferent geological structures than the primary events, see Figure 2.1. These meth-ods use a transform, e.g. F-K or Radon transform, to map the seismic data events,primaries and multiplies, into different regions, i.e. mapping different frequencycontents or different dip angles into different regions. Then, they eliminate themultiples events and transform the remaining events, the primaries, back to theoriginal time-offset domain (Verschuur, 2002).8Figure 2.1: Two reflection events: a primary (solid blue line) and a surface-related multiple (dashed red line). Both events have similar arrival timebut have traveled different paths. The primary event has traveled deeperinto the subsurface and has seen higher velocity than the multiple event.Therefore, it arrived with a smaller angle at the receiver.Data-driven methods exploit the relation between primaries and multiples eventsto suppress the repetitive patterns of multiples. These methods involve two steps: aprediction and a subtraction step. First, multiples are predicted and then subtractedthrough a matching procedure from total data, i.e. primaries, surface-related mul-tiples and internal multiples. Assumptions are needed for both steps. For examplethe energy after the subtraction step is assumed to be minimum (Verschuur, 2002).2.2 Surface-related multiples modelThis section discusses the model behind SRME and EPSI for 1D and 2D situation.We discuss theoretical details for SRME and EPSI in Section 2.3 and Section 3.1,respectively.92.2.1 Surface-related multiples formulation for normal-incidenceplane wave in a horizontally stratified mediumConsider a horizontal plane wave propagating in a laterally invariant earth. Thiswave reflects at the subsurface geology to produce surface-free data p0, whichcontains primaries and internal multiples but no surface-related multiples (bluesolid and red dashed lines in Figure 2.2a). At the surface, this surface-free dataacts as a new source of events that will feedback into the earth to produce first-order multiples (green dashed lines). The first-order multiples will then feedbackinto the earth again to produce higher order multiples. The feedback diagram inFigure 2.2b illustrates this process.We can express the surface-free data mathematically asp0(t)= g(t) s(t);where p0 is the surface-free data, g is the surface-free impulse response, denotesconvolution, and s(t) is the source signature. The feedback diagram in Figure 2.2permits us to formulate the total response with the surface related multiples p asp(t)= g(t) [s(t) p(t)]= p0(t) g(t) p(t): (2.1)Figure 2.3 illustrates the generation of surface-related multiples with the sourcesignature shown in Figure 2.3b, for a model with only one-reflection layer, Fig-ure 2.3c. The generated surface-free data p0, which consists of only one primary isshown in Figure 2.3d. This surface-free data feedbacks to the earth and convolveswith the surface-free impulse response g to produce first-order multiples m1, Fig-ure 2.3e. The same process produces second-order multiples m2 (Figure 2.3f). Fi-nally, the total response (surface-free data and surface-related multiples) is shownin Figure 2.3g.2.2.2 Data matrixOur seismic data in the 2D case is three-dimensional: time, source position andreceiver position. By considering the earth as a linear time-invariant system, wecan represent an experiment, i.e. a shot gather (one source and multiple receivers),10(a)(b)Figure 2.2: Feedback diagram for surface-related multiples generation. First,surface-free data is generated by convolving the source signature s+with surface-free impulse response g. Then, surface-free data feedbacksinto the earth to generate higher order multiples.as a superposition of many fully independent monochromatic experiments in theFourier domain. After taking the shot gather into the frequency domain and se-lecting a single frequency, we obtain a vector of measurements for the differentreceivers for that frequency. Repeating the process for all shot gathers produces adata matrix, indicated by bold capital P (Berkhout, 1985), whose columns repre-sent common source point gathers and rows represent common receivers gathers.2.2.3 Two-dimensional formulationConvolving surface-free impulse response G with total upgoing wavefield P gen-erates surface-related multiples. Figure 2.4 illustrates the surface-related multiplesprediction for 2D case interpreted as wave field extrapolation. We start with asource positioned at xs and a set of receivers (geophones) positioned at xk. Our11(a)(b) (c)(d) (e)(f) (g)Figure 2.3: a) An illustration that shows the generation of surface-relatedmultiples. b) Source signature s(t). c) Surface-free impulse responseg(t). d) Surface-free data p0(t). e) First-order multiples m1(t). f)Second-order multiples m2(t). g) Total response p(t).12Figure 2.4: After convolving the source gather from total data P with thecommon receivers gather from surface-free impulse response G for eachpoint in xk, the convolution result is summed up to produce the predictedsurface-related multiple trace at xr.goal is to calculate the wavefield at a receiver positioned at xr, which represent thesurface-related multiples. To do that, we need the Green’s function to extrapolatethe path for a source positioned at xr and whose energy is measured at xk.The surface-free impulse response G can act as the Green’s function to predictall possible surface-related multiples. First, we convolve the two wavefields, theshot gather from the total data P and the common receivers gather from the surface-free impulse response G, for each point in xk. Then, we sum the convolution resultto get the wavefield at xr and that gives the surface-related multiples.In practice, both the data P and the surface free impulse response G are trans-formed to the frequency domain and then the convolution and the summation arerepeated over all source-receiver combinations to obtain the predicted surface re-lated multiples M0. If we adapt a similar data matrix notation as in Berkhout(1985), then the multiple prediction can be formulated as a matrix-matrix multipli-cation:ˆM0 = ˆG ˆP; (2.2)where the minus sign represents the reflection at the surface and the hat is used torepresents a monochromatic matrix.13Furthermore, we can use the same approach in Equation 2.2 to generalize Equa-tion 2.1 for the 2D case to getˆP = ˆP0 + ˆM0; (2.3)where ˆP0 is the surface-free data and ˆM0 is the surface-related multiples. Expand-ing Equation 2.3 yieldsˆP = ˆG ˆQ+ ˆGR ˆP; (2.4)where ˆQ is the source signature and R is the surface reflectivity, which is ap-proximated to be I, i.e.ˆP = ˆG( ˆQ ˆP): (2.5)2.3 Iterative surface-related multiple removalIn this section, we aim to give a brief description of the iterative Surface-RelatedMultiple Removal (SRME) procedure (Verschuur et al., 1992). We start from thesurface-related multiples model, Equation 2.4. Now, let A = Q 1R be the sur-face operator. Next, we rewrite the surface-related multiples as ˆM0 = ˆGR ˆP =ˆG ˆQ ˆQ 1R ˆP = ˆP0AP. Substitution this expression in Equation 2.3 yieldsˆP = ˆP0 + ˆP0AP: (2.6)By assuming the source function Q to be omnidirectional, i.e. ˆQ = Iˆq with ˆqthe discretized Fourier representation of the source function, the above equationcan be solved iteratively as,ˆP(n+1)0 = ˆP A(n+1) ˆP(n)0 P; (2.7)where n is the iteration index. The iterative procedure starts with approximating thesurface-free data with total data, i.e. ˆP(0)0 = P. Next, it estimates A in a matchingprocess.This iterative procedure (Equation 2.7) matches predicted multiples to totaldata, which may lead to throwing out multiples energy. It also requires regularsampling.142.4 Curvelet-domain matched-filteringMany seismic processing techniques depend on matching seismic wavefield am-plitudes whether for adaptive subtraction of multiples (Verschuur et al., 1992; Her-rmann et al., 2008b) or ground roll (Yarham and Herrmann, 2008) predictions fromtotal data or for restoration of migration amplitudes through image to migrated-demigrated-image matching (Herrmann et al., 2008a). For a successful matching,we need to control a possible overfitting that may inadvertently lead to primaryenergy removal. Furthermore, we need to handle conflicting dips and apply theseparation after the matching stably. We will see in this section that curvelet sat-isfies all these conditions by its ability to decompose seismic data sparsely intoposition, scale and dip dependant coefficients and by imposing smoothness on itscoefficients in the phase space (Herrmann, 2008a).In Herrmann (2008a), the curvelet-domain matched-filtering is based on twoassumptions: A global Fourier matching removes the stationary differences between thetwo matched wavefields. The Fourier matching removes only amplitude-spectra mismatches and global kinematic errors and fails to remove errorsthat vary spatial as a function of dip. In SRME (Verschuur et al., 1992),these stationary differences come from the convolution of the source signa-ture during the multiples prediction step. The remaining non-stationary differences, which come from the non-idealcircumstances in the prediction, e.g. the absence of 3D effects, vary smoothlyas a function of position and dip along wavefronts.After compensating for the stationary differences via global Fourier matchingprocedure, i.e. no kinematic errors, we can approximate the non-stationary differ-ences mathematically by the operatorY, whose action on a d-dimensional functionf is given by(Yf)(x)=ZRde jk xa(x;z) ˆf(z)dz; (2.8)where x;z are the spatial coordinate and wavenumber vectors, ˆf(z) is the Fouriertransform coefficients of f(x), and a(x;z) is a space- and spatial-frequency depen-dent filter (Herrmann et al., 2008a).15In practice, the operator Y, which after discretizing becomes a full positive-definite matrix, acts on 2D images, e.g. shot gathers or common receivers gathers.We usually have two images to match, g and f . In the context of SRME, f wouldbe predicted multiples and g would be total data. Then, the mismatched can bemodeled as a matrix-vector multiplication, i.e.,g =Yf: (2.9)Following Herrmann et al. (2008a), we approximate the action of Y operator by apositive diagonal scaling in the curvelet domain(Yf)(x) C DYC f(x); (2.10)where C is the forward curvelet transform and C is its adjoint, which is also itspseudo inverse since we are using curvelet-domain via wrapping (Candes et al.,2006), which is a tight frame with C C = I. The diagonal matrix DY approximatesthe symbol a(x;z) evaluated at the curvelet centers.After substituting Equation 2.10 into Equation 2.9, the problem of finding theaction of the operator Y reduces to finding the elements of the curvelet diago-nal scaling matrix DY, which we calculate by solving the following least-squaresproblemargminz12jjg C diag(C f)zjj22; (2.11)where z represents the diagonal elements of DY.Many issues complicate the estimation of the curvelet diagonal scaling of Equa-tion 2.11. The forward model is underdetermined, which leads to mapping manyscaling coefficients to the same image under the adjoint operator C . There is alsothe possibility of overfitting that may lead to primary energy removal or incorrectamplitude corrections. To solve these issues, a smoothness constraint is imposedon the diagonal scaling to giveargminz12jjg C diag(C f)zjj22 +ljjLzjj22; (2.12)16withL =[LT1 LT2 LTq LTscale]T; (2.13)where l is the trade off parameter and L is the sharpening operator, which pe-nalizes fluctuations among neighboring coefficients in the diagonal scaling z. Thesharpening operator penalizes the scaling coefficients of each wedge (Figure 1.5)along the first and second direction, L1 and L2. It also penalizes coefficients be-tween adjacent wedges at the same scale Lq, and between wedges of adjacent scalesLscale, see Herrmann et al. (2008a); Shahidi and Herrmann (2009).To ensure positivity of the diagonal scaling, it is suggested in Herrmann et al.(2008a) to replace z with ez. This is not implemented in the curvelet-domainmatched-filtering used to produce our result and we leave its implementation tofuture works.2.5 MotivationTo illustrate the potential benefit of adding curvelet-domain matching, let us con-sider a common-offset section from North Sea field dataset. Figure 2.5 shows theresults of primaries estimation from two different multiples-elimination methods,namely Surface-Related Multiple Elimination (SRME, Figure 2.5b), and Bayesianthreshold (Figure 2.5c, Saab et al. (2007)), which uses curvelet-domain matched-filtering. Comparing the primary estimations from SRME and Bayesian thresh-old reveals improved continuity and amplitude preservation for the primaries (nearthe lower two arrows in each plot around 2.6 and 3.1s) of the Bayesian thresholdmethod. Furthermore, the strong residual of the first- and second-order water bot-tom multiples in the shallow part (near the top two arrows in each plot around 0.75and 1.2s) are better suppressed.17(a) (b)(c)Figure 2.5: A common-offset section from North Sea field dataset with re-sults from two different multiples-elimination methods. a) Total dataplotted with automatic-gain control. b) Primary estimation from oneterm SRME. c) Scaled Bayesian iterative threshold primaries estimate.Adapted from Herrmann et al. (2008b).18Chapter 3Convexified EPSI andcurvelet-domainmatched-filteringIn this chapter, we layout the theory behind estimation of primaries by sparse in-version (EPSI), before incorporating a matching step into the convexified EPSIprocedure.Will incorporating a matching step into EPSI optimization problem give betterresults than using EPSI alone, knowing that using EPSI one would match multiplesto multiples instead of matching predicted multiplies to the total data as in SRME?The remainder of this thesis is dedicated to answering this question.3.1 Estimation of primaries by sparse inversionEPSI is based on the same model as SRME. However, it estimates primaries inan inversion process rather than in a prediction and a subtraction process. UnlikeSRME, EPSI is not sensitive to limited sampling because it can reconstruct missingoffset during the inversion process.EPSI tries to change the paradigm that sees the surface-related multiples asnoise, which needs to be removed from the data prior to imaging. It proposes anoperator, which is based on the underlying principle that relates the primary im-19pulse response to total upgoing wavefield that includes the source signature andsurface related multiples. EPSI monochromatic relationship is expressed mathe-matically as:ˆG|{z}surface-free impulse responsedowngoing wavefieldz }| { ˆQ+R ˆP upgoing wavefieldz}|{ˆP ; (3.1)where the hat indicates a monochromatic representation of a wavefield arrangedinto a matrix, see Section 2.2.3 for details. The matrix product of any two-hattedwavefields corresponds to a non-stationary convolution in the time domain.The original EPSI (van Groenestijn and Verschuur, 2009) assumes the impulseresponse G to be sparse in the time domain and the source function Q is omni-directional, i.e. a constant-source for all shots or ˆQ = Iˆq with ˆq the discretizedFourier representation of the source function. Sparsity is enforced on the updatesof the surface-free impulse response by a zero norm. An increasing time windowin each iteration is placed over these updates in which the biggest events per traceare selected. Increasing the window size improves the convergence. A Fouriersmoothness (shortness) is also enforced on the source function.We obtain a reasonable estimate for primary impulse response G by choosingthe largest t (sparsity level) elements of the gradient and setting the rest to zeros.We keep repeating this process till we obtain the desired image of primaries.3.1.1 Convexified EPSIMany issues affect the practical adoption of the original EPSI formulation. It isdifficult to estimate many of its inversion parameters, e.g. the sparsity level of theprimaries in each iteration or the precise knowledge of the time window containingthem. The original EPSI also relies on the zero norm which makes it unstable.Lin and Herrmann (2009) address these issues via a ‘0 to ‘1 convexification,where the one norm on the surface-free impulse response is minimized alternatelywith the two-norm on the source function. Introducing ‘1 convexification will notonly preserve sparsity on the Green’s function but will also stabilize the problem. Itwill also eliminate most the inversion parameters. The time window, which is usedto prevent the inversion from being trapped in a local minima, is no longer needed20since convex problems have no local minima. The sparsity level of primaries isreplaced with a more practical and easy to estimate parameter s, which is relatedto the noise level of measured data, see Equation 3.3.3.1.2 Bi-convex optimizationEPSI optimization involves two alternating optimization problems that solve forthe unknown surface-free Green’s function G and the unknown source function Q.To be consistent with the notation of optimization problems, indicated by lowerbold face letters e.g. p = vec(P) in our notation, we rewrite Equation 3.1 in avectorized form in this bi-linear forward model formulationE[ ˆQ]g :=F t blockdiag ˆQ+R ˆP 1 nf I Ftg = p; (3.2)where E[ ˆQ] is a linear operator, which depends on the source signature Q and theupgoing wavefield P that acts on the vectorized surface-free Green’s function g, denotes the conjugate transpose, and denotes the Kronecker product. Thisexpression reformulates the matrix-matrix product into a matrix-vector product.Ft is the Fourier transform, whose action on a vectorized wavefield g is definedas Ftg := ˆg = [ˆg1;ˆg2; ;ˆgn f ] with n f the number of frequencies and F t is theinverse Fourier transform operator that brings the vectorized wavefield back to thetime domain. The blockdiag lays out each frequency slice (1 n f) in a vectorizedform and p is the vectorized upgoing wavefield. For now, we will assume R = I.Equation (3.2) is linear in both the surface-free Green’s function and the Fouriercoefficients of the source function. This permits us to invert for these two un-knowns by alternatingly solving two optimization problems. First, we solve forfree impulse response g subject to data fitting, i.e. the one norm optimization prob-lem on the unknown surface-free Green’s function,eg = argmingjjgjj1 subject to jjp E[ ˆQ]gjj2 s; (3.3)whereeg is the vectorized estimation (estimations are indicated bye) for the surface-free Green’s function and s is the residual, which is linked to the noise level of thedata. Then, we solve the regularized least-squares problem for the unknown source21function,eˆq = argminq12jjey B[eˆG]ˆqjj22 +lFjjLF ˆqjj22; (3.4)whereey = vec([ˆP+eˆG ˆP]1 n f ) is the data vector with eˆG = vec 1(eg) the current es-timate for the surface-free Green’s function and B[eˆG] := blockdiag(eˆG1 n f I). Inthis expression, we assume the source to be omnidirectional, i.e. Q1 n f = Iq1 n fThe sharpening operator LF in the ‘2 penalty term imposes smoothness on the es-timated Fourier coefficients of the source functioneˆq, which corresponds to enforc-ing decay in the time domain. The trade off parameter lF balances Fourier domainsmoothness versus data misfit to avoid overfitting, which may lead to leakage ofmultiples into the source function.In this bi-convex optimization, we start with estimate for the source functionand then solve for the surface-free Green’s function by an one-norm optimizationproblem, which seeks a sparse vector that after convolution with the downgoingwavefield explains the total upgoing wavefield, see Equation (3.3). Given thissparse estimation for the Green’s function, the regularized least-squares problem issolved that seeks the Fourier coefficients of the source function, see Equation (3.4).A sparsfying TransformTo further enforce the sparsity assumption, we can define a suitable sparsifyingtransform on the estimated surface-free Green’s function. In Lin and Herrmann(2009), they define this sparsifying transform as a Kronecker product between the2D curvelet transform along source-receiver coordinates, and the wavelet trans-form along time coordinate, i.e. S := C W, with C, W the curvelet and wavelettransform operators, respectively. The synthesis operator S then sits between E[ ˆQ]and g in Equation 3.3, and the solution seeks the sparsest set of transform domaincoefficients of the surface-free Green’s function Seg instead of the sparse Green’sfunctioneg.Due to the redundancy of the curvelet transform, about 8 times redundant in2D case, and the high demand for memory of our implementation of Algorithm 1,see Section 3.2, the above sparsifying transform is not used in our matching.223.2 EPSI with curvelet-domain matched-filteringEPSI implicitly assumes that the source is omnidirectional, i.e. constant sourcefor all shot gathers; absence of 3D effects; infinite aperture; and ideal total surfacereflectivity, i.e. R = I. Although in principle EPSI can be extended to 3D, it isbeen shown that this assumption is more problematic to handle in practice as expe-riences with SRME-based method have shown (Lin et al., 2004). Deviation formthese assumptions may lead to detrimental effects on the quality of the estimatedsurface-free data.Following the successful application of curvelet-domain matched-filtering withinthe SRME procedure, we propose to incorporate a matching step into our bi-convexoptimization procedure to mitigate possible adverse effects that stem from ampli-tude mismatches that vary smoothly as a function of position and dip along thepredicted wavefronts of the current estimate for the surface-free data. We includea non-trivial surface reflectivity operator, i.e. R6= I, that will absorb amplitudeerrors of non-ideal reflections at the surface, finite-aperture, and other unknowneffects.We start with the assumption that the operator R is symmetric positive defi-nite and pseudo local (no kinematic shifts). Then, we approximate this space-dip-dependent filter by a simple curvelet domain scalingR C diag(z)C: (3.5)After that, the diagonal curvelet scaling is estimated in the regularized least-squaresproblem Equation 2.11 that matches multiples to multiples and not predicted mul-tiples to total data as in SRME.3.2.1 Frequency-domain monochromatic matchingOur first attempt to incorporate the surface reflectivity operator R into our formula-tion was to assume that this operator is frequency independent, even though it canbe easily extended to include angular frequency dependency. Then we incorporatethis operator in Equation 3.1 to getˆP ˆG ˆQ = ˆGC diag(z)C ˆP; (3.6)23which we can solve by least-squares to get the curvelet scaling z:argminz12jjeu ˆMzjj22 +ljjLzjj22; (3.7)where ˆM :=eˆGC diag(C ˆP) andeu :=vec([ˆP eˆGfˆQi]) with i the index correspondingto the frequency for whicheˆq is maximum.This method solves for the curvelet scaling fast, 200x faster than the methodproposed in the next section, and does not demand large memory to execute. How-ever, it does not give the desired space-dip dependent amplitude corrections corre-sponding to the symbol of the pseudo differential operator a(x;z) in Equation 2.8since it is acting on a frequency slice that has two spatial coordinates and no timecoordinate. To solve this problem, the scaling should act on a spatial-time slice,e.g. a shot gather or a common receivers gather.3.2.2 Time-domain matchingIn this approach, we apply the spatial-dip dependant operator R to each shot recordof the upgoing wavefield P separately, to give the desired spatial-dip dependentamplitude corrections. Then, we take the scaled upgoing wavefield to the Fourierdomain and substitute it into EPSI monochromatic relationship Equation 3.1, to getvec(ˆP eˆGeˆQ)= ˆG ˆMz1 ns; (3.8)where ˆMi = C diag(CPi) with i the shot gather index and ns is the number ofsources. Our operator R now acts on each shot record of the upgoing wavefield Pseparately.This new approach seeks a smooth positive spatial-dip dependent curvelet-domain scaling for all the shot gathers of the total upgoing wavefield simultane-ously, which makes it time and memory demanding approach, especially with thecurvelet redundancy of about 8 times for the 2D case.Incorporation of this expression into our formulation, yields a tri-convex opti-24mization problem that now includes curvelet-domain matching:z1 ns = argminz1 ns12jjeu Nz1 nsjj22 +ljjLz1 nsjj22; (3.9)where N := blockdiag([eˆGeˆM]1 n f I) and eu := vec(ˆP eˆGeˆQ]1 n f ). The trade-offparameter l and the sharpening operator L are of Equation 2.12 respectively andns and n f are the number of sources and the number of frequencies respectively.To solve this tri-convex optimization, we start by assuming R = I. Then,given an estimate for the source signature we solve Equation 3.3 for the sparsesurface-free Green’s function. After that, we substitute our Green’s function esti-mate into Equation 3.4, and solve the least-squares problem to obtain a new sourceestimate.Now that we have an estimate for both G and Q, we can start the matchingstep by assuming R6= I and solve Equation 3.9 for the diagonal curvelet-domainscaling. The new estimate for the surface reflectivity operator R is then used in thenext estimations of the surface-free Green’s function and the source function, seeAlgorithm 1.Algorithm 1 Stabilized Estimation of Primaries by Sparse Inversion with CurveletDomain Matched FilteringResult: Estimate for g and Qchoose noise level s, iteration increment rg0 0, Q0 0, R0 I, k 1whilekp E[ ˆQ;R]gk2 s dogk solve (3.3) using initial guess Qk 1;Rk 1 with SPG‘1 1at least rk itera-tions.Qk solve (3.4) for the source signature using gk;Rk 1 with LSQRz solve (3.9) with LSQRRk C diag(z)Cend whileg gkQ Qk1see van den Berg and Friedlander (2008)25Chapter 4Observations and resultsThis chapter studies the factors affecting the outcome of incorporating a matched-filtering step into EPSI. We conclude this chapter by showing our results.4.1 Methodology4.1.1 Synthetic dataWe consider synthetic time-domain finite-difference acoustic forward modelingdata with 128 shot gathers each with 128 receivers (traces). Each trace has 256time samples and a time sampling interval of 6.4ms; and the distance between anytwo sources or receivers is 20m. Delphi f dacmod was used to generate this data,see Table 4.1 for detailed parameters. We assume the data to be receiver ghost free.Unless otherwise stated, all the shot gathers shown in this section are the middleshot gather, i.e. shot index of 64 (offset = 1260m).4.1.2 Results evaluationWe want to compare the primary multiples separation results of EPSI with andwithout the curvelet-domain matched-filtering step. However, since we do not havethe true primaries to evaluate the separation results, e.g. with signal-to-noise ratio,we rely on prior knowledge of the locations of these primaries to evaluate ourresults.26Table 4.1: Synthetic time-domain finite difference acoustic forward modelingdata parameters.Parameter Value Parameter Valuetime sample interval 6.4ms sources depth 5msamples/trace 256 receivers depth 5mnumber of sources 128 grid size (dx) 5mnumber of receivers 128 grid size (dz) 5mwavelet type Ricker wavelet model size (nx) 509frequency range 0-60Hz model size (nz) 251source type monopoleThe matching step will be evaluated by the smoothness and the positivity ofthe curvelet-domain matched-filtering diagonal scaling.4.2 Factors affecting our resultsIn this section, we study the factors that control the output of the matching process,namely the number of least-squares iterations h and the trade off parameter lwhich control the smoothing and positivity of the curvelet-domain scaling. We willalso answer the question of when and how frequent should we do the matching.Before experimenting with these factors, we present EPSI results without match-ing of the synthetic data described in Table 4.1. Figure 4.1a shows the zero-offsetsection of total upgoing wavefield. After applying EPSI to this data, it convergesafter 22 alternating optimization loops to give an estimate for the primaries (Fig-ure 4.1b) and an estimate for the source wavelet (Figure 4.1d). The multiples es-timate (Figure 4.1c) is calculated by subtracting the primary estimate form totalupgoing wavefield.4.2.1 The smoothing and positivity of the curvelet-domain diagonalscalingIn Section 2.4, the second assumption behind the curvelet-domain matched-filteringstates that the non-stationary differences vary smoothly as a function of locationand dip and then it can be modeled mathematically by a matrix-vector multiplica-27(a) (b) (c)(d)Figure 4.1: Zero-offset EPSI synthetic data results with no matching. a) To-tal upgoing wavefield section ( ˆP). b) Primaries estimation ( ˆP0). c) pri-maries estimation minus the total data ( ˆP ˆP0 = ˆG ˆP), i.e. the multi-ples estimation d) Final estimated source wavelet ( ˆQ)tion of a positive definite matrix on a reference vector, see Equation 2.9.To test the effects of the number of least-squares matching iterations h ofEquation 3.9 and the trade off parameter l on the positivity and smoothness ofthe curvelet-domain scaling vector z1 ns, we carry out two experiments on whichwe fix one of the parameters and allow the other to vary. Then, we comparethe smoothness of the curvelet-domain scaling of the middle shot gather (offset= 1260m) with a mosaic plot, see Figure 1.5 for detailed description. We also show2D-plots of the same curvelet-domain scaling to judge the positivity. After that wepresent the final zero-offset primary multiple separation results, i.e. the estimated28primaries ˆP0 and the calculated multiples ˆM0.Fixing l and varying hIn this experiment, we test the effect of the number of curvelet-domain matched-filtering iteration h on the positivity and smoothness of the curvelet-domain scal-ing vector by fixing the trade off parameter l = 0:46 and varying the number ofiteration h = 10;35 and 100.Figure 4.2a, 4.2c and 4.2e show mosaic plots of the curvelet-domain scalingcoefficients for h = 10;35 and 100, respectively and Figure 4.2b, 4.2d and 4.2fshow the corresponding 2D-plots of the same coefficients. Since the surface reflec-tivity coefficients are negative, approximated by R= I in EPSI without matching,and to avoid the confusion with literature (Herrmann et al., 2008a; Yarham andHerrmann, 2008) where the curvelet-domain scaling vector z1 ns is positive, wewill plot z1 ns instead of z1 ns. The corresponding final primaries and multiplesestimations are shown in Figure 4.3(a)-(f).As the number of iterations increases, the curvelet-domain scaling vector be-comes more smooth and positive as shown in Figure 4.2. The smoothness andpositivity also affect the final EPSI output. In Figure 4.3a, we notice a significantleakage of multiples into the zero offset primaries section, indicated by the smallblack arrows at around 0.48, 0.75, 0.8 and 1s. However, in Figure 4.3b, theseindicated multiples are either eliminated or faded. Figure 4.3c has slightly lessmultiples than Figure 4.3b. For example, the multiple, indicated by the third arrow,has been eliminated. Figure 4.3(d)-(f) show the same effects, so as the number ofiterations increases the multiples become more preserved.Increasing the number of iterations produces a more smooth and positive curvelet-domain scaling vector, which results in a better primary-multiples separation. How-ever, that comes at the expense of execution time. Although we are utilizing theparallel execution and the distributed arrays of MATLAB in our implementation ofAlgorithm 1, one matching iteration still takes about 10 minutes to execute on 5nodes each with 8GB memory shared between 2 workers.29(a) l = 0:46;h = 10 (b) l = 0:46;h = 10(c) l = 0:46;h = 35 (d) l = 0:46;h = 35(e) l = 0:46;h = 100 (f) l = 0:46;h = 100Figure 4.2: Mosaic and 2D plots for the curvelet-domain matched-filteringscaling for one shot (offset = 1260m) for l = 0:46 and h =10;35 and 100. The l and h used to generate each subfigure are statedunderneath it..30(a) l = 0:46;h = 10 (b) l = 0:46;h = 35 (c) l = 0:1;h = 100(d) l = 0:46;h = 10 (e) l = 0:46;h = 35 (f) l = 0:46;h = 100Figure 4.3: Zero-offset synthetic data results for l = 0:46 and h = 10, 35and 100. The l and h used to generate each subfigure are stated under-neath it. a-c) zero-offset primaries estimates. d-f) zero-offset multiplesestimates.31Fixing h and varying lThis experiment aims to test the effect of the trade off parameter l on the positiv-ity and smoothness of the curvelet-domain scaling vector by fixing the number ofmatching iterations h = 35 and allowing the trade off parameter to vary l = 0.046,0.46 and 46.4.The mosaic plots of the curvelet-domain scaling coefficients for l = 0.046,0.46 and 46.4 are shown in Figure 4.4a, 4.4c and 4.4e respectively. Next to eachmosaic plot, a 2D-plot of the same coefficients is shown, Figure 4.4b, 4.4d, and4.4f. The corresponding final primaries and multiples estimations are shown inFigure 4.5(a)-(f).Increasing the trade off parameter l will obviously make the curvelet-domainscaling smoother and more positive, see Figure 4.4, but l also balances the curvelet-domain smoothness versus data misfit to avoid overfitting, which may lead to leak-age of multiples into primaries estimate. The later point is obvious in Figure 4.5,where we notice severe leakage of multiples (Figure 4.5f) into the primaries (Fig-ure 4.5c). Figure 4.5a has slightly less multiples than Figure 4.5b especially at thedeepest parts near the lowest three arrows.A large value of the trade off parameter l and the number of matching iter-ations h are desirable, but if not chosen correctly might lead to severe leakageof multiples into primaries estimate or to longer execution time. Hence, in orderto practically incorporate the matching step, we need to choose these parameterswisely to avoid data overfitting and timely execution.4.2.2 When and how frequent should we do the matching?In this section, we carry out different experiments to determine the best way toincorporate curvelet-domain matched-filtering step into EPSI. We present theseexperiments starting with the worst result. Then, we evaluate the quality of eachresult by comparing it to EPSI result with no matching presented in Figure 4.1. Inall the experiments, the number of iteration h and the trade off parameter l arefixed, h = 30 and l = 0:46.Before we carry out the experiments, we will take a look at different alternat-ing optimization iterations snapshots (Figure 4.6). In particular, we look at the32(a) l = 0:046;h = 35 (b) l = 0:046;h = 35(c) l = 0:46;h = 35 (d) l = 0:46;h = 35(e) l = 46:4;h = 35 (f) l = 46:4;h = 35Figure 4.4: Mosaic and 2D plots for the curvelet-domain matched-filteringscaling for one shot (offset = 1260m) for h =35 and l = 0.046, 0.46 and46.4. The l and h used to generate each subfigure are stated underneathit.33(a) (b) (c)(d) (e) (f)Figure 4.5: Zero-offset synthetic data results for h = 35 and l = 0.046, 0.46and 46.4. The l and h used to generate each subfigure are stated under-neath it. a-c) zero-offset primaries estimates. d-f) zero-offset multiplesestimates.34primaries estimates for different iterations ˆP0 = ˆG ˆQ, Figure 4.6a- 4.6d, the mul-tiples estimations calculated by subtracting the primaries estimates from the totalupgoing wavefield MS0 = ˆP ˆG ˆQ, Figure 4.6e- 4.6h, and finally at the multiplesestimates generated by convolving the Green’s functions with the total upgoingwavefield ˆMC0 ˆG ˆP, Figure 4.6i- 4.6l.Referring to Figure 4.6, we notice the importance of the first few iterations;they define the way that the upcoming iterations follow. After these few iterations,we notice only small changes in the inversion process, i.e. small changes in ˆP0, ˆMS0and MC0 . For example, there is a slight change between the primaries estimate afterthe sixth iteration (Figure 4.6c) and after the twenty-first iteration (Figure 4.6d)compared to a significant change between the primaries estimate after the first iter-ation (Figure 4.6a) and after the third iteration (Figure 4.6b).Notice also that after these few iterations, the two multiples estimates, ˆMS0 andˆMC0 start to become closer to each other. In curvelet-domain matched-filtering,we are trying to estimate the surface-reflectivity by taking advantage of EPSI bymatching multiples to multiples instead of predicted multiples to total data like inSRME. Hence, it is desirable to have ˆMS0 close to ˆMC0 , but how close? That is whatwe will try to find out in the following experiments.Only one matching at the first alternating optimization loopIn this experiment, we only solve one curvelet-domain matched-filtering optimiza-tion problem, Equation 3.9, right after we get the first estimation of the Green’sand the source function. Then, we use the estimated curvelet-domain scaling in allthe upcoming iterations.Figure 4.7a- 4.7b show the zero-offset primaries and multiples estimates, re-spectively. There is obvious leakage of multiples into the primaries estimate, indi-cated by the small black arrows. The primaries estimate after the first alternatingloop is still unstable and will go under a lot of changes during the first few loops,which might lead to this leakage. On the other hand, we notice a leakage of pri-maries into the multiples estimate, which might come form the fact that the twomultiples ˆMS0 and ˆMC0 are far apart from each other when we do the matching.Therefore, incorporating a matched-filtering step into EPSI in the first few iter-35(a) Iteration = 1 (b) Iteration = 3 (c) Iteration = 6 (d) Iteration = 21(e) Iteration = 1 (f) Iteration = 3 (g) Iteration = 6 (h) Iteration = 21(i) Iteration = 1 (j) Iteration = 3 (k) Iteration = 6 (l) Iteration = 21Figure 4.6: a-d) Primaries estimates at different iterations. e-h) Multiples es-timates calculated by subtracting the primaries estimate for the total dataˆP ˆG ˆQ. i-l) Multiples estimates calculated by convolving the Green’sfunction with the total upgoing wavefield ˆG ˆP.36(a) (b)Figure 4.7: Zero-offset results of EPSI with only one matching step at thefirst alternating optimization loop. a) Zero-offset primaries estimate.b) Zero-offset multiple estimate calculated by subtracting the primariesestimate form the total upgoing wavefield.ations produces an undesirable leakage of events between primaries and multiplesestimates.Multiple matching starting from the second alternating optimization loopIn this experiment, we solve six curvelet-domain matched-filtering optimizationproblems at even loops starting from the second loop, i.e. in loops indexed 2, 4, ...,12.The zero-offset primaries and multiples estimates are shown in Figure 4.8a-4.8b, respectively. We immediately notice that the primaries are better preserved37(a) (b)Figure 4.8: Zero-offset results of EPSI with six matching steps starting afterthe second alternating optimization loop (even loops only). a) Zero-offset primaries estimate. b) Zero-offset multiple estimate calculatedby subtracting the primaries estimate form the total upgoing wavefield.and there is much less primaries leaking into the multiples. However, there are ob-viously still multiples leaking into the primaries estimate, indicated by the arrows.This leakage might be coming form the first two matching steps at the second andfourth loop where the inversion process still goes under significant changes.This result is better than the previous one. However, solving more than onematched-filtering optimization problem is impractical. Solving one matched-filteringoptimization problem takes longer than solving the same problem without thematching step.38(a) (b)Figure 4.9: Zero-offset results of EPSI with only one matching step at thetwenty first alternating optimization loop. a) Zero-offset primaries es-timate. b) Zero-offset multiple estimate calculated by subtracting theprimaries estimate form the total upgoing wavefield.Only one matching near the last alternating optimization loopIn this experiment, we only solve one matched-filtering optimization problem nearthe last loop, loop 21, just before the inversion problem converges.The zero-offset primaries and multiples estimates are shown in Figure 4.9a-4.9b, respectively. This result is very close to the original result with no curvelet-domain matched-filtering but it suffers from weak multiples events leaking intothe primaries estimate, indicated by the black arrows. In the last few loops, theinversion problem goes through small changes and the two multiples estimates ˆMS0and ˆMC0 are very close and there is no much that the matching step can correct.39(a) (b)Figure 4.10: Zero-offset results of EPSI with only one matching step at thetwelve alternating optimization loop. a) Zero-offset primaries esti-mate. b) Zero-offset multiple estimate calculated by subtracting theprimaries estimate form the total upgoing wavefield.Only one matching at the middle alternating optimization loopIn this experiment, we solve one matched-filtering optimization problem near themiddle loop, loop 12. At this stage in the inversion problem, the dramatic changesof the first few iterations have already stabilized and we still have many upcomingiterations before the inversion problem exits.Figure 4.10a- 4.10b show the zero-offset primaries and multiples estimates,respectively. This result is better than all previous results including the one withoutcurvelet-domain matched-filtering one. We notice less multiples leakage into theprimaries estimate especially where the two black arrows indicate.404.2.3 Optimized resultWe will conclude this section by first presenting an optimized result, a result wherel and h are chosen wisely and the matched-filtering optimization is solved onlyonce in the middle loop, and then compare it to EPSI result without the curvelet-domain matched-filtering.The optimized result is shown in Figure 4.11. Figure 4.11a- 4.11c show thezero-offset section of the total upgoing wavefield, the estimate for the primariesand the calculated multiples, respectively. The estimated source wavelet is shownin Figure 4.11d.This result is better than EPSI without the curvelet-domain matched-filteringstep presented in Figure 4.1. By comparing this result primaries estimate (Fig-ure 4.11b) to the one without matching estimated primaries (Figure 4.1b), we no-tice that the multiple event indicated by the top arrow is now weaker and moreinterestingly the other multiple indicated by the bottom arrow is no longer there.The number of iterations (h = 60) and the trade off parameter (l = 0:49) arechosen to avoid overfitting the data and to solve the least-squares optimizationproblem in a reasonable time and at the same time does not sacrifice the positivityand smoothness of the curvelet-domain scaling vector. The matching step is incor-porate at the twelve loop away from the instability of the few first iteration to avoidleakage between the primaries and the multiples.41(a) (b) (c)(d)Figure 4.11: Zero-offset EPSI synthetic data result with curvelet-domainmatched filtering incorporated into loop 12. The number of iterationsh = 60 and the trade off parameter l = 0:46. a) Total upgoing wave-field section ( ˆP). b) Primaries estimation ( ˆP0). c) primaries estimationminus the total data ( ˆP ˆP0 = ˆG ˆP), i.e. the multiples estimation d)Final estimated source wavelet ( ˆQ)42Chapter 5Conclusion and future works5.1 ConclusionThe aim of this thesis is to incorporate a curvelet-domain matched-filtering stepinto the formulation of EPSI and investigate the factors that affect the outcome ofthis matching step on the outcome of EPSI.We started by leveraging the curvelet transform ability to represent seismic dataand images sparsely as a function of location, scale and dip to end up with a solidmatched-filtering formulation that has the ability to control a possible overfittingthat might lead to primary energy removal.Then we incorporated this formulation into EPSI formulation to produce a tri-convex optimization problem in which after we invert for the Green’s and sourcefunction, we solve our matching-filtering problem.After that we looked at the different factors that affect the outcome of ourinversion namely, the number of iteration h, the value of the trade off parameterl and when and how frequent should we solve the matched-filtering problem. Wealso saw that the number of iterations and the value of the regularization parametershould be chosen wisely to balance the smoothness and positivity of the curvelet-domain scaling versus data overfitting and long execution time.Then, we concluded by showing that solving the curvelet-domain matched-filtering step around the middle alternating optimization loop gives the best resultssince in that stage the inversion problem have already stabilized and the dramatic43changes of the first few iterations have already done.5.2 Future worksIn this section we try to improve our method by making it faster. We suggesttwo ways to do that either by changing the problem formulation or by imposingfrequency regularization over the curvelet-domain scaling vector.New formulationThe problem with the current formulation Equation 3.9 is that we have to solve forthe curvelet domain scaling for the whole data at the same time which requires agreat deal of time and memory. We start by assuming R = I and then change theformulation a bit to getP GQ = C diag(z)CGP: (5.1)This formulation is solved for each shot gather separately which makes it fast andeasy to parallelized. However, this formulation does not have a physical meaninglike matching for the reflectivity coefficients (Equation 3.9).Frequency regularizationIn Shahidi and Herrmann (2009), they show that by imposing a smoothness overthe curvelet-domain scaling vector in the frequency domain, faster convergence isrealized. By incorporating this in our formulation we could shorten the executiontime of Equation 3.9.44BibliographyA. Berkhout. Seismic migration : imaging of acoustic energy by wave fieldextrapolation. Elsevier, 1985. !pagesA. Berkhout. Multiple removal based on the feedback model. The Leading Edge,18:127–131, 1999. !pagesE. Candes, L. Demanet, D. Donoho, and L. Ying. Fast discrete curvelettransforms. Multiscale Modeling and Simulation, 5:861–899, 2006. !pagesE. J. Candes and D. L. Donoho. Curvelets: a surprisingly effective nonadaptiverepresentation of objects with edges. In in Curve and Surface Fitting:Saint-Malo. University Press, 2000. !pagesM. R. Gadallah and R. L. Fisher. Applied Seismology: A Comprehensive Guide toSeismic Theory and Application. PennWell Corporation, 2005. !pagesG. Hennenfent and F. J. Herrmannn. Seismic denoising with nonunformlysampled curvelets. Computing in Science and Engineering, 8:16–25, 2006. !pagesF. J. Herrmann. Curvelet-domain matched filtering. In SEG Technical ProgramExpanded Abstracts, volume 27, pages 3643–3647. SEG, November 2008a. !pagesF. J. Herrmann. Curvelet-domain matched filtering. In SEG Technical ProgramExpanded Abstracts, volume 27, pages 3643–3647. SEG, November 2008b. !pagesF. J. Herrmann and G. Hennenfent. Non-parametric seismic data recovery withcurvelet frames. Geophysical Journal International, 173:233–248, 2008. !pages45F. J. Herrmann, P. P. Moghaddam, and C. C. Stolk. Sparsity- andcontinuity-promoting seismic image recovery with curvelet frames. Appliedand Computational Harmonic Analysis, 24:150–173, 2008a. URLhttp://slim.eos.ubc.ca/Publications/Public/Journals/HerrmannACHA07.pdf.Accepted for publication in the Journal of Applied and ComputationalHarmonic Analysis. !pagesF. J. Herrmann, D. Wang, and D. J. Verschuur. Adaptive curvelet-domainprimary-multiple separation. Geophysics, 73(3):A17–A21, 2008b.doi:10.1190/1.2904986. URL http://slim.eos.ubc.ca/Publications/Public/Journals/herrmann08acd/paper html/paper.html. !pagesD. Lin, J. Young, Y. Huang, and M. Hartmann. 3d srme application in the gulf ofmexico. SEG Technical Program Expanded Abstracts, 23(1):1257–1260, 2004.doi:10.1190/1.1851098. URL http://link.aip.org/link/?SGA/23/1257/1. !pagesT. Lin and F. J. Herrmann. Unified compressive sensing framework forsimultaneous acquisition with primary estimation. In SEG, 2009. !pagesR. Neelamani, A. I. Baumstein, D. G. Gillard, M. T. Hadidi, and W. L. Soroka.Coherent and random noise attenuation using the curvelet transform. TheLeading Edge, 27:240–248, 2008. !pagesC. C. Paige and M. A. Saunders. Lsqr: An algorithm for sparse linear equationsand sparse least squares. ACM Trans. Math. Softw., 8:43–71, March 1982.ISSN 0098-3500. doi:http://doi.acm.org/10.1145/355984.355989. URLhttp://doi.acm.org/10.1145/355984.355989. !pagesR. Saab, D. Wang, O. Yilmaz, and F. Herrmann. Curvelet-based primary-multipleseparation from a bayesian perspective. In SEG International Exposition and77th Annual Meeting, 2007. URL http://slim.eos.ubc.ca/Publications/Public/Conferences/SEG/2007/saab07seg.pdf.!pagesR. Shahidi and F. J. Herrmann. Curvelet-domain matched filtering withfrequency-domain regularization. In SEG Technical Program ExpandedAbstracts, volume 28, page 3645. SEG, 2009. !pagesE. van den Berg and M. P. Friedlander. Probing the pareto frontier for basispursuit solutions. SIAM Journal on Scientific Computing, 31(2):890–912, 2008.doi:10.1137/080714488. URL http://link.aip.org/link/?SCE/31/890. !pages46G. J. A. van Groenestijn and D. J. Verschuur. Estimating primaries by sparseinversion and application to near-offset data reconstruction. Geophysics, 74(3):A23–A28, 2009. doi:10.1190/1.3111115. URLhttp://link.aip.org/link/?GPY/74/A23/1. !pagesD. J. Verschuur, A. J. Berkhout, and C. P. A. Wapenaar. Adaptive surface-relatedmultiple elimination. Geophysics, 57:1166–1177, 1992. !pagesE. Verschuur. Seismic multiple removal techniques: Past, present and future.EAGE Publications BV, 2002. !pagesC. Yarham and F. J. Herrmann. Bayesian ground-roll seperation bycurvelet-domain sparsity promotion. In SEG Technical Program ExpandedAbstracts, volume 27, pages 2576–2580. SEG, November 2008. !pages47Appendix AImplementationsThere are many implementations of the curvelet-domain matched-filtering namely,a MATLAB, a Python and a partially MPI C++ implementation. Although the lastimplementation is the fastest and has more features, we chose to extend the MAT-LAB implementation for these reasons. MATLAB does not only make developingand testing algorithms fast and easy, but a recent extension to MATLAB, namelyParallel Computing Toolbox, makes solving memory and time intensive algorithmspossible. It provides MPI free programming of distributed arrays, different parts ofone large array reside at different computational nodes, and distributed algorithmexecution, different parts of an algorithm is executed at different computationalnodes at the same time. The latter feature makes our algorithm, see Section 3.2.2,possible. For example, the diagonal scaling in Equation 3.9 for our small syntheticdata has more than 32 million elements.SPG‘1, a gradient-based solver that uses spectral step length and Pareto rootfinding (van den Berg and Friedlander, 2008), is used for solving the ‘1 optimiza-tion problem in Equation 3.3 and for the ‘2 optimization problem in Equation 3.4,the least-squares solver described in Paige and Saunders (1982) is used.We have chosen the curvelet-domain transform via wrapping (Candes et al.,2006) as our curvelet transform with 6 scales and 16 angles at the finest scale. Forthis choice of implementation, the pseudo-inverse equals the adjoint.48
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimation of surface-free data by curvelet-domain...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimation of surface-free data by curvelet-domain matched filtering and sparse inversion AlMatar, Mufeed H. 2011-05-09
pdf
Page Metadata
Item Metadata
Title | Estimation of surface-free data by curvelet-domain matched filtering and sparse inversion |
Creator |
AlMatar, Mufeed H. |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | A recent robust multiple-elimination technique, based on the underlying principle that relates primary impulse response to total upgoing wavefield, tries to change the paradigm that sees surface-related multiples as noise that needs to be removed from the data prior to imaging. This technique, estimation of primaries by sparse inversion (EPSI), (van Groenestijn and Verschuur, 2009; Lin and Herrmann, 2009), proposes an inversion procedure during which the source function and surface-free impulse response are directly calculated from the upgoing wavefield using an alternating optimization procedure. EPSI hinges on a delicate interplay between surface-related multiples and pri- maries. Finite aperture and other imperfections may violate this relationship. In this thesis, we investigate how to make EPSI more robust by incorporating curvelet- domain matching in its formulation. Compared to surface-related multiple removal (SRME), where curvelet-domain matching was used successfully, incorporating this step has the additional advantage that matches multiples to multiples rather than predicated multiples to total data as in SRME. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-05-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution 3.0 Unported |
DOI | 10.14288/1.0053426 |
URI | http://hdl.handle.net/2429/34367 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by/3.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2011_fall_almatar_mufeed.pdf [ 2.95MB ]
- Metadata
- JSON: 24-1.0053426.json
- JSON-LD: 24-1.0053426-ld.json
- RDF/XML (Pretty): 24-1.0053426-rdf.xml
- RDF/JSON: 24-1.0053426-rdf.json
- Turtle: 24-1.0053426-turtle.txt
- N-Triples: 24-1.0053426-rdf-ntriples.txt
- Original Record: 24-1.0053426-source.json
- Full Text
- 24-1.0053426-fulltext.txt
- Citation
- 24-1.0053426.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0053426/manifest