PLANE-WAVE DECOMPOSITION AND RECONSTRUCTION OF SPHERICAL-WAVE SEISMOGRAMS AS A LINEAR INVERSE PROBLEM by JOSE JULIAN CABRERA B.Sc. Engineering Geophysics, Institute* Politecnico Nacional (Mexico),1981 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department Of Geophysics And Astronomy) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1983 © Jose Julian Cabrera, 1983 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. Department of Geophysics And Astronomy The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: 5 December 1983 ABSTRACT t The plane-wave decomposition of the vertical displacement component of the spherical-wave field corresponding to a compressional point source is solved as a set of inverse problems. The solution utilizes the power and stability of the Backus & Gilbert (smallest and flattest) model-construction techniques, and achieves computational efficiency through the use of analytical solutions of the integrals which are involved. The theory and algorithms developed in this work allow stable and efficient reconstruction of the spherical-wave field from a relatively sparse set of their plane-wave components. However, the algorithms do not formally conserve the correct amplitudes of the seismograms. Comparison of the algorithms with direct integration of the Hankel transform shows very little or no advantage for the transformation from the time-distance (t-x) domain to the delay time - angle of emergence (r-y) domain if the seismograms are equi-sampled spatially. However, for cases where the observed seismograms are not equally spaced and for the transformation r-y to t-x, the proposed schemes are superior to the direct integration of the Hankel transform. Applicability of the algorithms to reflection seismology is demonstrated via the solution to the problem of trace interpolation and that of separation of converted S modes from other modes presented in common-source gathers. In both cases the application of the algorithms to a set of synthetic reflection seismograms yields satisfactory results. i v TABLE OF CONTENTS ABSTRACT ii LIST OF TABLES viLIST OF FIGURES viii ACKNOWLEDGEMENTS X 1 INTRODUCTION TO THE PLANE-WAVE DECOMPOSITION PROBLEM 1.1 Expansion of a Spherical-Wave Field in Terms of Plane and Cylindrical Waves: Basic Development 1 1.2 On Plane-Wave Decomposition of Digital Data 15 2 BACKUS AND GILBERT FORMULATION OF THE PROBLEM 2.1 Forward Transform: Construction of the Plane-Wave Components of Spherical-Wave Seismograms 20 2.1.1 Smallest model construction 21 2.1.2 Flattest model construction 5 V 2.2 Inverse Transform: Reconstruction of the Spherical-Wave Field from its Plane-Wave Components 28 2.2.1 Smallest model construction 30 2.2.2 Flattest model construction3 EXAMPLES 3.1 Introductory Comments 32 3.2 Separation of Converted S Modes 34 3.3 White Noise and Construction Quality 52 3.4 Angular Sampling and Trace Interpolation 56 3.4.1 Shallow zone of interest 57 3.4.2 Trace interpolation 61 4 COMPUTATIONAL CONSIDERATIONS 66 5 SUMMARY 69 BIBLIOGRAPHY 71 APPENDIX A X2 Value and Observational Errors in Model Construction 73 APPENDIX B Effects of the Weighting Function Q and Standard Deviation Values on Model Construction 76 v i APPENDIX C Inner Product Matrix for the Forward Smallest Model Construction 86 APPENDIX D The Forward Flattest Model: Global Development 87 APPENDIX E Inner Product Matrices for the Inverse Smallest and Flattest Model Constructions 91 vii LIST OF TABLES I. Basic Steps for PWD of a Common-Source Gather 17 Ii. Processing Times for the PWD Algorithms 66 viii LIST OF FIGURES 1.1a Propagating spherical-wave front 4 1.1b Plane waves with different propagation velocities ... 5 1.2 Mode expansion in terms of plane waves 6 1.3 Definition of the wavenumber components 7 1.4 Displacement potential as a superposition of plane waves 9 3.1 Model used for generating the synthetic seimograms .. 35 3.2a P-P seismograms 37 3.2b P-S seismograms 8 3.3 PP+PS seismograms 39 3.4 Signals with large and small moveouts 40 3.5 Predicted times for PWD components 2 3.6 PWS Hankel transform 43 3.7 PWS converted S zone 4 3.8 Converted S modes. Hankel-Hankel 46 3.9 Converted S modes. Smallest-Smallest 47 3.10 Converted S modes. Flattest-Flattest 8 3.11 Converted S zone. Hankel-Flattest 50 3.12 Converted S zone. Smallest-Smallest 1 3.13 PP+PS seismograms. Additive random noise 53 3.14 Converted S modes. Smallest-Smallest. (Noisy seismograms) 54 ix 3.15 Converted S modes. Flattest-Flattest (Noisy seismograms) 55 3.16 PP+PS seismograms. (Near-offset resampling) 58 3.17 PWS smallest model. (Near-offset resampling) 59 3.18 PWS flattest model. (Near-offset resampling) 60 3.19 PP+PS seismograms. Uneven Spacing 63 3.20 Interpolated S modes. Smallest-Smallest 64 3.21 Interpolated S modes. Flattest-Flattest 65 B.1 Modified Bessel functions K0 and K, 77 B.2 The condition number of the inner product matrix .... 80 B.3 Eigenvectors for the (smallest) inner product matrix 81 B.4 The 25th basis function 82 B.5 Effect of standard deviations on the basis functions 85 X ACKNOWLEDGEMENTS I am deeply grateful to my ] /] 1 ~j ^ P (dear, friend) Shlomo Levy for all his valuable advice, knowledge and moral support during the entire development of this thesis. I regard him with profound affection and look forward to a perdurable friendship. It is my pleasure to warmly acknowledge Dr. Ron Clowes, my supervisor, for his wise direction and continuous encouragement. Certainly I have developed much respect for Ron. It is a bit difficult to express, in a language other than my native, my gratitude to Dr. Doug Oldenburg. He has provided me with guidance and counselling, and shared his superb knowledge of Inverse Theory. I hope to continue working with Doug both as student and as "amigo". I wish to acknowledge the generous advice of Dr. George Bluman of the Department of Mathematics at UBC, for his very helpful (and calming) discussions, particularly in relation to integrals of arbitrary and modified Bessel functions. Enthusiastic thanks are also due to Michael Shlax, Kerry Stinson, Kenneth Whittall and Ian Jones for providing me with so much insightful advice. xi I wish I could continue listing all the extraordinary people who have contributed to the nice environment in which I have studied, but space precludes me from doing so. However, I would particularly like to mention Lynda Fisk, Mark Lane, Don White and Gemma Jones. My Master's studies were supported by a graduate scholarship from the Science and Technology National Council (CONACYT) of Mexico, and Operating Grant A7707 from the Natural Sciences and Engineering Research Council of Canada. LA MEMORI A DE MEMITO 1 1. INTRODUCTION TO THE PLANE-WAVE DECOMPOSITION PROBLEM 1.1 Expansion of a Spherical-Wave Field in Terms of Plane and Cylindrical Waves: Basic Development For many years the study of wave fields produced by spherical waves travelling in simple earth models has been facilitated by analysis of the plane waves representing the original spherical waves. Indeed, for a medium consisting of homogeneous layers, reflection, refraction and mode conversion at layer interfaces is simpler to investigate using plane rather than spherical waves. For example, the method of computing synthetic seismograms due to Fuchs and Muller (1971) is based on solving the wave-propagation problem for plane waves and then superimposing the plane-wave solutions to obtain the spherical-wave field (i.e. the synthetic seismograms). The representation of a scalar, time-harmonic, spherical wave field (SWF) as a superposition of plane waves is well documented in the seismic and electromagnetic literature; see for example, Stratton (1941), Bath (1968), Goodman (1968), Born and Wolf (1980), and Aki and Richards (1980). Lucid treatments of the general theory are given by Brekhovskikh (1960), and Devaney and Sherman (1973), and some insightful applications are 2 presented by Asby and Wolf (1971), Muller (1971), and Treitel et al. (1982). Following Aki and Richards (1980; ch.6) we introduce the Weyl plane-wave and the Sommerfeld cylindrical-wave expansions of a SWF. An outline of the mathematical steps necessary for obtaining such representations follows. Consider the problem of a point source at x=0 radiating compressional waves in a homogeneous, isotropic and unbounded medium. Given that the source exhibits time dependence of the form exp[-icjt] (u> is an arbitrary angular frequency), compressional-wave propagation may be described by the (scalar) displacement potential <p, which satisfies d2<j>/d2t - V2V20 = 4TTVp 26(x)expt-iwt] (1-1) where VP represents the P-wave velocity of the medium. The space-time solution to (1-1) is *(x,t) = [ 1/R]exp[io;(R/VP- t)] (1-2) with, AAA x=xi+yj+zk R = /x2 +y2 + z2 Equation (1-1) may also be solved using Fourier transform methods, in which case the wavenumber-time solution reads as tf(k,t) = [4irVp 2/(k2Vp 2-£j2) ]exp[-iut] (1-3) 3 with k = kx i + k3 j + k , The relation between </>(x,t) and <p(k,t) is given by the triple inverse Fourier transform c>(x,t) = [1/R]exp[ iu(R/Vp- t)] = [ 1/2TT2 ]exp[-iut]///[ l/(k2-o;2/VP2) ]exp[ ik- x]dk •too (1-4) -co Equat ion (1-4) represents o spherical wave </>(x,t) travelling with a constant velocity vp as a superposition of an infinite number of homogeneous plane waves, each of which propagates with a velocity w/k (see Figures 1.1a and 1.1b). Because 0<cj/k<o°-, these plane waves do not all travel with the velocity Vp of the medium and hence, equation (1-4) does not represent a true mode expansion of a spherical wave in terms of plane waves. In order to obtain plane waves propagating with the same velocity Vp, one of the integrals in (1-4) is evaluated. The usual choice is the k* integral, which is solved by extending kj to the complex plane and then applying residue theory (Aki and Richards, 1980, p. 195-197; Devaney and Sherman, 1973, p. 770-775). The integration yields <MCJ,X) = [ 1/R]exp[ iwR/Vp ] -too = [1/2ir] //[ 1/ik, ]exp[ ikt z]exp[ ikx x + ik^y]dk^dky (1-5) 4 where k? = ±/w2/V 2"k2-k2 and the sign of kj is chosen positive for z>0 and negative for z<0 (notice that we have eliminated the terms exp{-iwt] and explicitly stated the CJ dependence in the argument of <p) . FIGURE 1.1a Propagating spherical wave front in a medium of constant velocity Vp. For convenience we have shown only the x-z plane. Expression (1-5) is known as the Weyl integral and it represents a spherical wave as a superposition of plane waves travelling with constant velocity Vp (see Figure 1.2). For kt real, i.e. k2 + kv2<aj2/Vp 2, the plane waves are homogeneous and propagate along the direction specified by the wavenumber vector 5 -A. * A A k = k„ i + k., j + kj. k. However, for kt imaginary, these waves are inhomogeneous and travel along directions specified by kxi + kv,j with an exponential amplitude attenuation in the z direction. Z FIGURE 1.1b Propagating plane waves. The integrand in equation (1-4) represents plane waves travelling away from and toward the point x=0 (source location). For a given angular frequency a>, each plane wave has an arbitrary velocity V;''=a>/k; , where k;=|kf |. Thus, for kb<ka. the plane wave represented by the wavenumber vector^k^ travels slower than the plane wave characterized by kb. Transformation of (1-5) to cylindrical coordinates gives the Sommerfeld integral, which reads as <Mcj,r,z) = [ 1/Vr2+z2]exp[i<j(v/r2+z2)/Vp ] CO = /[ 1/iki jexpfikj z]kr J0(kr r)dk. (1-6) o with 6 r y/x2+y2 (1-7a) wsin(7)/Vp (!-7b) ±y/u2/Vf 2-kr 2 = ±CJCOS(7)/VP (1-7C) and 7 denotes the polar angle of the wavenumber vector k (Figure 1.3), and J0 is the zero-order Bessel function of the first kind. FIGURE 1.2 Mode expansion in terms of plane waves. In this case the plane waves travel with the medium velocity V. The restriction on the sign of kt requires that the plane waves propagate away from the point 0. Thus, the homogeneous plane waves shown in this figure have positive vertical wavenumbers kj. The Sommerfeld integral represents a spherical wave as a superposition of cylindrical waves given by krJ0(krr)exptik^z]. These waves are also weighted by the term 1/ik4 and they have the same vertical propagation factor (namely exp[ik4z]) as the x z 7 plane waves. A j FIGURE 1.3 Definition of the wavenumber components. We will be concerned only with the x-z plane and z>G. So far we have dealt with an impulsive point source at x=0. For cases in which the point source has a time dependence F(t) with a spectrum F(CJ), equation (1-6) is modified to 0(o),r,z) = JF(CJ) [ 1/ikj ]exp[ ik„ z]kr J0 (krr )dkr (1-8) o Up to this point we have considered an unbounded homogeneous space. For a medium consisting of a sequence of homogeneous layers, the total displacement potential at an observation point P located at some depth z is given by the 8 contribution from the direct and reflected potentials (recall that these potentials propagate as spherical wave fronts; see Figure 1.4). In this case, then, the corresponding plane waves (in view of 1-5) or cylindrical waves (in view of 1-6) are weighted by plane-wave reflection and transmission coefficients. Additionally, within the itt layer each plane (or cylindrical) wave is characterized by the vertical wavenumber k,r>} = y/cj2/V2 ~ k2 = ucos (7;p )/Vrp (1 - 9a) for P waves. By similar arguments the wavenumber for a S wave, k*; = V/CJVV;5 2 - kr2 = ucos(7;s )/V:S (1-9b) may be found. The horizontal wavenumber is kr= wsin (y* )/Vp, and a superscript letter denotes the wave mode. The displacement potential at P(r,z) is then obtained from (see Figure 1.4) <p = tf>° + 01+tf>2 (1-10) where each potential is expressed as (cf. equation 1-8) 0°(cj,r,z)= /F(o)) [ 1/ik; ]exp[ikj;1z] krJ0(krr)dkr (1 -11 a) o 01 (cj,r,z)= /F(o))A[ 1/ik; ]exp[ikpl/1h1+ikPil (h,-z) ] o krJ0(krr)dkr (1-1 lb) 9 vfv,' Source, vfv* Plr,z) z = 0 O r V,P 7< vP (C) fe) FIGURE 1.4 Displacement potential at the observation point P(r,z) as a superposition of plane waves, (a) The displacement potential at P is given by <p0 + 4>1 +<P2 • (b) A given emergent plane wave with ray parameter p=sin(7' )/Vp. This ray parameter defines a family of reflected and transmitted plane waves, so that in (c) the direct plane wave contributes to 4>°, in (d) the plane wave affected by the reflection coefficient A contributes to the potential 01, and in (e) the plane wave affected by the transmission coefficients B and D and the reflection coefficient C contributes to <p2. 10 <t>2 (u,r , z) = fF(u)BCDexp[ i kpy,h , + i kp 22h2 + i k£ , (h ,-z ) ] [ i/ik*a ]krJ0(krr)dkr . (1-1 1c> The vertical wavenumbers k^, and kp;2 are given by (1-9a) and, because the compressional source is located in the first-layer, ki=ki(1. A and C are plane-wave reflection coefficients, while B and D are plane-wave transmission coefficients. They depend on the elastic properties of the medium and on the ray parameter p (recall that p=cjkr =cjsin (7, )/VP ) . From (1-10) and (1-11) the function defined by V(cj,kr;z) = F(co)exp[ iki,z] + F (o>) Aexp[ i k\t, h , + i k^ , (h ,-z ) ] + F(w)BCDexp[ ik;-1h, + ikp,22h2 + ikP.1 (h,-z) ] (1-12) will be understood as the spectrum of the plane waves defined by kr =cjsin {y* )/Vp . For a given angular frequency w, V(w,kr;z) gives both the plane-wave contribution from the boundaries between layers, and the vertical phase delay that the (homogeneous) plane waves acquired in each layer. If k\x is imaginary, an inhomogeneous plane wave will propagate horizontally in the i1* layer, and the contribution from the bottom interface to this layer will have an exponential attenuating term. Therefore, with reference to Figure 1.4, equation (1-12) gives the response from a direct plane wave (first term .in the left hand side), from a plane wave reflected from the first 11 boundary (second term) and from a plane wave reflected from the second boundary (third term). Notice that these plane waves form part of a system of plane waves sharing the same horizontal wavenumber. In this context, then, V(a>,kr;z) is viewed as the spectrum of a plane-wave seismogram equivalently defined by kr, P or 7f, . Some additional insight into the plane-wave nature of (1-12) may be seen in the time domain. Before inverse-Fourier transforming this equation, it is convenient to make the substitution (see 1-9a) k\- = wcos(7, )/V;p i = 1 ,2 Then, we may rewrite (1-12) as V(u,y\;z) = F(u>)exp[ icjzcos(7p, )/Vp ] + F(u.)Aexp[iw{h,cos(7i )/Vp + (h,-z)cos(7i )/v'} ] + F(w)BCDexp[ icj{h!COs(7')/Vp + 2h2cos(7f2 )/Vp2 + (h1-z)cos(7ei)/V?1}]. (1-13) Notice that we have now used y\ to characterize V(&j,kr;z). For homogeneous plane waves (i.e. for 7.p real), inverse Fourier transform of (1-13) gives 1 2 VCt^w-z) = F(t - zcos(7Pi )/V\) + A F(t - {htCOsO'J/Vj+d^-zJcosfTp/V'}) + BCD F(t - {h,cos(7f, )/V,,1+2h2cos(7?2)/V2 + (h,-z)cos(7?1)/V,}]. (1-14) For a given earth model, this equation demonstrates that the plane wave seismogram V(t,7?,;z) is given by the source function F(t) appearing at delay times dependent on the cosines of the angles of propagation of plane waves within each layer. Further, this source function is weighted by plane-wave reflection and transmission coefficients. Thus, the (homogeneous) plane-wave response for the direct potential <£°(t,r;z) corresponding to the angle is Similarly, the plane-wave response for the displacement 0'(t,r;z) is given by Notice that we have used "plane-wave response" to describe those time functions whose spectra constitute complex weights for plane or cylindrical waves (see equations 1-11). If the observation point P(r,z) is at the surface z=0, the plane-wave response for the direct potential is simply the g0(t,7r,.;z) = F(t-td ) = F(t - zcos(Vi)/V^) (1-15a) gMt,7i;z) = AF(t-tA ) = AF(t - {h,cos(Vi )/Vi + (h,-z)cos(7F, )/V?,}). (1-15b) 1 3 source function with no delay time, i.e. g°(t,7^;z) = F(t) for all angles 7?, . Summarizing and generalizing the results, for a sequence of homogeneous layers with a compressional source F(t) placed at (r=0,z=0), observation point P(r,z>0) and given angular frequency o>, the representation of the total displacement potential at P in terms of a superposition of cylindrical waves is given .by (see 1-10 to 1-13) *(u,r,z) = JV(w,kr ;z) [ 1 / i k\ ]kr J0(krr)dkr . (1-16) o where V(<j,kr;z) is the generalization of (1-13) and it represents the spectra of plane-wave seismograms, each of which may be characterized by either the horizontal wavenumber kr, the ray parameter p, or the emergent angle 7^. These parameters are. related by kr = us in ( y\ )/v\ = cop The inverse Fourier transform of V(w,kr;z) gives plane-wave seismograms characterized by either p or 7F, (but not by kr ) . Because kinematically these time-domain plane-wave seismograms constitute delayed versions of the source function F(t) (see 1-14 and 1-15), the time argument of V will be denoted by the delay time r. Schematically, then, 1 4 V(CJ, 7; z) V(CJ, kr ; z ) V(w,p;z) INVERSE FOURIER TRANSFORM V(r,7;z) V(r,p;z) (for notational convenience we have replaced y\ by 7). The representation of the compressional displacement at the point P(r,z) in terms of cylindrical waves is obtained by applying the gradient operator to (1-16). In particular, the vertical displacement S(w,r,z) is given by S(cj,r,z) = 3/9z{tf>(w,r ,z)} = J3/3z{V(u,kr;z)} [1/ik;]krJ0(krr)dkr (1-17) o In equation (1-13) it is seen that the z-dependence of V(w,kr;z) is given only by the terms exp[iicjzcos(7^ )/Vp] = exp[±ikpz]. Hence (1-17) gives S(o>,r,z) = /U(cj,kr ;z)[ik; ][\/ik\ ]krJ0 (krr )dkr (1-18) o and because of the ± signs in the argument of expfik^z], U(o>,kr;z) and V(cj,kr;z) may differ only in the sign associated to each of their terms. If the observation point is located at the surface z=0, equation (1-18) becomes S(w,r) = /U(aj,kr )J0(krr)krdkr o (1-19) 1 5 Plane-wave decomposition (PWD) is understood as the problem of computing the plane-wave seismograms U(u>,kr) from the spherical-wave observations S(o,r). Muller (1971) recognized equation (1-19) as a zero-order Hankel transform and consequently, presented its formal inversion as U(cj,kr) = fs(cj, r)J0(krr)rdr (1-20) o providing the basic formulation for PWD. 1.2 On Plane-Wave Decomposition of Digital Data The introduction of slant stacking (Schultz and Claerbout, 1978) and plane-wave decomposition (Treitel et al. 1982) into the realm of reflection seismology has revealed a large number of possible applications in exploration seismology. In a recent publication, Treitel et al., (1982) have shown the relationships between slant stacking and plane-wave decomposition. They pointed out that although equation (1-20) is restricted to the vertical component of compressional waves recorded over laterally homogeneous media, this formalism seems to provide ' an acceptable approximation of a larger class of earth models, in particular those consisting of dipping layers. The general procedure followed to perform the plane-wave decomposition of spherical-wave seismograms resulting from a 1 6 common-source gather is illustrated in Table I. The application of PWD in real seismic work requires the discretization of equation (1-20). The form suggested by Treitel et al. (1982) is, U(w,kr) = ArLSdd.r. )J0(krr- )rf (1-21) In this form, equation (1-21) may be applied to a set of seismograms S(o>,r; ) found at evenly spaced intervals Ar = r;,, -r; . In this case, the' integration increment Ar is factored out of the summation and introduced later as a form of global scaling. When uneven seismogram spacing is encountered, we observe that the formal use of equation (1-21) with variable Arf causes deterioration of the decomposition due to inappropriate weighting of certain terms in the summation. In this case it might be better to evaluate (1-20) using an appropriate numerical integration scheme (e.g. interpolating S(<j,r;) and then using 1 -21). An alternative approach to plane-wave decomposition of discretely sampled data has been proposed by Henry et al. (1980). In this approach, PWD is considered as a linear inverse problem and a smallest model for U(a>,kr) is constructed. Henry et al.'s solution is efficient and is applied directly to the forward transformation from time-offset domain to the delay time - angle of emergence domain (i.e. from t-x to r-y; because we are concerned with one dimensional earth models, we make no distinction between the space variables x and r). The same approach may be used for the inverse transformation TABLE I 1 7 BASIC STEPS FOR PWD OF A COMMON-SOURCE GATHER COMMON-SOURCE GATHER Input of N seismic traces S(t,r;). FORWARD FOURIER TRANSFORM Temporal Fourier transformation of each trace in the common-source gather. Calculation of S(uu,r.-). PLANE-WAVE DECOMPOSITION At each angular frequency u>, computation of U(UJ, kr) (or equivalently of U(LU,JT)) for M different angles of emergence, INVERSE FOURIER TRANSFORM (with respect to<^) From U(Lu,kr), computation of M plane-wave seimograms U(^,T). 18 (i.e. from the r-y domain to the t-x domain). Since for some applications PWD may serve as a filtering operation, our objective is to modify Henry et al.'s solution and produce a more flexible and stable algorithm which will perform both the forward (t-x to r-y) and the inverse (r-y to t-x) transformations. Rather than use the inner product given by Henry et al., we introduce an explicit weighting function into (1-19) OO S(u,r) = /[U(w,kr )Q" 1 ][QJ0UR r)kr ]dkr , o for the forward transformation, and into (1-20) U(u,kr) = /[S(u,r)Q- 1 ][QJ0U.r)r]dr o for the inverse transformation. This is an example of a linear quelling (Backus, 1970), and allows the use of the usual definition of the inner product of two functions. With this, in addition to find the smallest model solutions for U(co,kr)Q"1 and for S(cj,r)Q"', we find the flattest model solutions. Finally, the solutions of the forward and inverse transformations are found subject to the x2 criterion (see Appendix A), so that observational errors are accounted for. The algorithms are applied to the problems of separating converted S modes from other modes in a common-source gather, 19 and of trace interpolation. In the first problem we will follow Tatham et al. (1983), and transform the common-source gather (t-x domain) to the plane-wave domain (r-y domain). For reasons which will be outlined later, certain converted modes will occupy a distinct portion of the plane-wave domain. Inverse transformation of only this portion back to t-x space will yield the common-source gather (CSG) representation of the S modes present in the chosen r-y zone. The second problem is solved by utilizing the algorithms to construct additional seismograms at offsets not represented in the original CSG. It is important to emphasize that although the use of (1-19) and (1-20) are restricted to the recorded compressional waves at the surface z=0, we realize that kinematically these equations are still satisfactory for obtaining the plane-wave signature of recorded S waves. The reason for this is understandable from the discussion developing (1-12). In the case of S plane waves, V(u>,kr;z) has vertical propagation terms of the form exp[ icjzcos(7s, )/V* ] and hence, the delay times are still governed by cosine functions. In Chapter 3 we will return to this matter. 20 2. BACKUS AND GILBERT FORMULATION OF THE PROBLEM 2.1 Forward Transform (t-x to r-y): Construction of the Plane-Wave Components of Spherical-Wave Seismograms To apply the Backus-Gilbert (E-G) theory to the problem of plane-wave decompositon, we use equation (1-19) at specified offset r- and angular frequency u, that is S(w,r,- ) = /U(cj,kr )J0(krr; )krdkr (2-1) o We can now solve (2-1) for U(cj,kr) as a set of inverse problems each of which corresponds to a given angular frequency co. In order to expedite the following presentation we introduce the terminology and notations to be used throughout the remainder of this work: (a) S(cj,r: ), the temporal Fourier-transformed elements of the spherical-wave seismograms, at a given angular frequency CJ and offset r- , are termed 'observations' and are denoted by e?. (b) U(w,kr), the temporal Fourier transform (FT) of the plane-wave seismograms at a given angular frequency u>, is termed the 'model' and is denoted by m; m is a continuous function of the horizontal wave number kr.. 21 (c) J0(krr;)kr, the zero-order Bessel functions multiplied by the horizontal wave number, are termed 'kernels' and are denoted by G;. They are continuous functions of the horizontal wave number kr . (d) The inner product of the functions f(k) and g(k) is CO denoted by <f,g> ( i.e., <f,g> = /f(k)g(k)dk ). o In the following two sections we will outline the B-G solution to problems of the form of equation (2-1) to show how smallest and flattest models can be computed. Further treatment of the procedure is found in Parker (1977), and Oldenburg and Samson (1979). 2.1.1 Smallest model construction (forward transform) Consider the problem, e? = <G; ,m> i = 1 , . . . ,N (2-2) where N is the number of observations. Assume that the given observations are contaminated by errors {5e;} with zero mean and estimated standard deviation a,-, i.e. e? = et±6e,- , e* being the true data. Therefore the equation to be solved is, e*±5e,- = <G,- ,m> which upon division by a- becomes, [e;*±6e; ]/o; = <G; /o; ,m> 22 e? = <Gr,m> (2-3) e° are our new observations with associated errors of unit variance, and G; are the new scaled kernels. Given N observations e? and the functional form (2-3), we would like to obtain a continuous model m. This problem is always underdetermined and admits infinitely many solutions but a specific model is obtained by minimizing some norm of the model and using the data equations as constraints. The smallest model corresponds to the requirement that the L2 norm of the constructed model will be smaller than that of any other permissible solution (i.e. all those satisfying (2-3)). Given this requirement, the solution to (2-3) can be expressed as a linear combination of the kernels, that is (Oldenburg and Samson, 1979), The coefficients af are obtained by substituting (2-4) into (2-3), changing the order of summation and integration and solving the system, m = Ia:G; (2-4) e° = Ta (2-5a) that is, a = T-1e° (2-5b) where a is the vector of kernel coefficients, e° is the vector of observations, and 23 T~1 is the inverse of the (NxN) inner product matrix T defined by, r,-= <G,-,G;> (2-6) The formal solution given in equations (2-4) to (2-6) cannot yield a physical solution if the kernels G,- are not square integrable. The current problem of PWD is an example of this occurrence. A way to circumvent this is to use the quelling operation (Backus,1970) which is essentially a mapping of the kernels into a Hilbert space. The method we use is named "quelling by multiplication" (Backus,1970) in which we look for a weighting function Q such that G? = G,-Q is in L2(0,°°) for all i. Once such a Q is specified we rewrite (2-3) as, e? = <GrQ,m/Q> = <G^,m+> (2-7) and continue to find the smallest model m* as outlined in equations (2-4) to (2-6). Subsequently we "unweight" m* and obtain the desired model. The final solution is then given by, m = La; Q2G; (2-8) An important consideration in the choice of the function Q is the ease with which the evaluation of the inner product r;j-= <G*,G/> can proceed. The efficiency of the construction algorithm increases greatly if an analytical expression representing the elements of the inner product matrix"is found. Also, the weighting function Q should lead to an efficient 24 construction of the inverse transform (i.e. from 7-7 to t-x). Indeed, since many inverse solutions are required (one for each frequency), numerical efficiency is gained if Q is chosen so that only a single matrix spectral decomposition is required. If this objective can be achieved, the introduction of the x2 criterion (Appendix A) for noisy data does not decrease the algorithm's efficiency. Our solution to the smallest-model forward construction involves the weighting function Q = kr" °-5Kc0-5 (kr b) , where K0 is a modified Bessel function of zero order, and b is an arbitrary positive real number whose role is demonstrated in Appendix B. With the above choice of weighting the construction proceeds with the following identifications, nr = m/Q = kr °-5K0-°-s (krb)U(cj,kr ) (2-9a) and Gr = GrQ = kr °-5K00-5 (kr b) J0 (kr r: )/o; (2-9b) Hence, from equation (2-8) the frequency "representation of the plane-wave seismograms is given by, U(<j,kr ) = L[a;/o; ]K0(krb)J0(krr; ) (2-10) where N is the number of seismograms in the common-source gather. 25 The evaluation of the inner product <G;Q,GjQ> is described in Appendix C, whereas determination of the coefficients a; to obtain proper x2 value is discussed in Appendix A. 2.1.2 Flattest model construction (forward transform) Consider the problem outlined in the previous section, i.e. given N observations corresponding to N functional relations e° = <G,- ,m>, find a model m which satisfies these relations- In this section, we search for the model which exhibits the least amount of change with respect to the independent variable. The construction of this type of model (commonly referred to as the flattest model) is achieved by the minimization of the norm ||m'||, with m' being the derivative of the model. In the problem of PWD, the flattest model requires that the Fourier transform of the PWS exhibits the least amount of variation with respect to kr at each angular frequency o>, and hence it may yield better lateral continuity in terms of both amplitude and time delay. In order to construct the flattest model we integrate the r.h.s. of equation (2-3) by parts to obtain, OO e? = H,- m| - <H- ,m'> (2-1 1 ) o where, H;(kr) = /G;(u)du 26 Presuming we can evaluate the term R;m| we substract it from the o left hand side to get, e,t = -<H; ,m'> (2-12) Using the technique outlined in section 2.1.1 above we proceed to find the smallest m' model. If the new kernels H; are not square integrable, we introduce a weighting function Q and then solve e;t = <-H;Q,m'/Q> to obtain m'. The solution to this problem is given by (cf. equation 2-8), m' = -Z/3- Q2H,- (2-13a) I - r where the 0,- are obtained from, I = r-'lt (2-13b) and, r..3 = <H,Q,Hj Q> (2-13c) Taking the indefinite integral of equation (2-l3a) we obtain, m(kr) = -Z/3j JQ2 (u)R,- (u)du + C (2-14) Two important considerations should be emphasized at this po i n t: (a) the construction of the flattest model necessitates the additional knowledge of a (boundary) value of m, from which the constant C is found. 27 (b) the choice of the weighting function Q is now burdened by the additional evaluation of the indefinite integral JH.Q2. In the solution to the forward-transform flattest-model construction we make the following identifications: (a) G; = krJ0(krr: )/ot . (b) The new kernels -H;= -krJ,(krr, )/[o;r; ] (see Appendix D). (c) Limiting r; to be greater than zero, we have H;m|=0. On o the other hand, assuming a band limited source function, kr-*<» implies infinitely-attenuated inhomogeneous waves. Hence we consider H;m| =0 (see Appendix D), which means that the new observations e-f are the same as ef. (d) The weighting function Q we have chosen is K,°-5(krb) where is a modified Bessel function of first order and b is an arbitrary positive real number (see Appendix B). We notice that because kr=wsin(7)/V, a different weighting function Q is used for each angular frequency u>. (e) The constant of integration C is equal to zero (see Appendix D). Following the development in Appendix D, the plane-wave seismograms at a given angular frequency are therefore given by, U(u,kr) = -L0; JQ2H,- = |[0j /o; r,]{bkr J, (kr r; )K0(krb) + r. kr J0 (krr,- )K, (krb)} /[r;2+b2] (2-15) 28 2.2 Inverse Transform (r-y to t-x): Reconstruction of the Spherical-Wave Field from its Plane-Wave Components The inverse problem to be solved is expressed in equation (1-20). kr has been replaced by r as the independent variable, and the model and observations have traded places so that the former now represents the temporal FT of the spherical-wave seismograms S(w,r) whereas the latter consists of the FT of the plane-wave seismograms U(&>, kr ).. Indeed, there is no difference between constructing plane-wave and spherical-wave seismograms from each other in the way formulated in the previous section. However, because we have formulated the problem in terms of the horizontal wavenumber kr rather than in terms of the ray parameter p, the inner product matrix for the inverse transformation explicitly depends on the angular frequency (see Appendix E). To see this and compare the form of the inner product matrices for the forward and inverse transformations, let us consider the following integrals r;F- = jG;(krr; )G- (krr_; )Q2(krb)dkr (2-16a) o and T-j = ;G;(kr.r)Gj (kfjr)Q2(br)dr. (2-16b) r;F- is an element of the inner product matrix for the forward transformation while r,* is an element for the inverse transformation. In terms of the ray parameter p, (2-16a) and (2-!6b) read as 29 f oo r.- = /G; (upr; )Gj (cjprj )Q2 (upb)dup (2-17a) o and r?i = o/G; (up. r)G- (up^rjQMbrJdr. (2-l7b) These integrals have the same form. If in (2-l7a) we set b=rc =positive constant, and in (2-17b) b=cjpc , with p„=positive constant we find rf- = /G; (cjpr; )G0 (cjprj )Q2 (uprc )dwp (2-18a) and r?j = 1/CJ /G; (£jp; r)Gj(cjpJr)Q2(upcr)dcjr. (2-l8b) Therefore, integration of (2-l8a) with respect to cjp and (2-l8b) with respect to ur will give matrices r* and rA which can have multiplicative factors dependent on u. In both transformations, spectral decomposition of a single inner product matrix is done only once. For the inverse transformation the definition of b=a>pc in the argument of Q means that, as in the forward problem, a different weighting function is used at each angular frequency CO. 30 2.2.1 Smallest model construction (inverse transform) The procedure here is parallel to that outlined in the section dealing with the corresponding forward transform. Identifying, G; = rJ0(rkr; )/a; Q = r-°-5K0°-s (rb) and using equation (2-8), the smallest model for the spherical-wave seismograms at a given u> is, S(w,r) = Z[a;/c; ]K0 ( rb) J0 ( r kr; ) with b = wsin(c)/V (2-19) where M is the number of plane-wave seismograms. 2.2.2 Flattest model construction (inverse transform) The development here is similar to that of the corresponding forward transform. At the stage of the solution of the smallest model m', the substitution b = cjsin(c)/V is made (Appendix E). Assuming S (a>, r=°°) =0, the integration constant C in equation (2-14) is zero. Hence the frequency representation of the vertical component of the spherical-wave seismograms is given by, 31 b = S(u,r) = Z[/3,/a, kr;]{brJ, (kr;r)K0(rb) i - t + kr;rJ0(kr;r)K1 (rb)}/[ k2+b2 ]} us in(c)/V (2-20) where M is the number of plane-wave seismograms, 32 3. EXAMPLES 3.1 Introductory Comments In the forward transform our goal is to obtain plane-wave seismograms U(r,y) from spherical-wave seimograms S(t,x). To achieve this, we follow the steps illustrated in Table I of section 1.2 (p. 17). Conversely, in the inverse transform our objective is to compute the vertical-displacement spherical-wave seismograms S(t,x) from plane-wave seismograms U{r,y). To achieve this we still follow those steps given in Table I though keeping in mind that the input data are plane-wave seismograms U(r,7). We refer to formation of the plane-wave or spherical-wave seismograms via the discrete form of equation (1-20) or equation (1—19) as the Hankel transform, and construction involving the B-G inversion as the smallest or flattest model. In the examples to follow we will obtain plane-wave seismograms for values of the angle of emergence y between 0° and 90° (that is to say 0<kr<cj/V). Similarly, the reconstruction of spherical-wave seismograms will include plane waves within that range only. Our motivation for doing so is based on the following: (a) We do not consider surface waves in our synthetic seismograms. 33 (b) We consider the receivers to be placed sufficiently far away (i.e. many wavelenghts) from the reflecting interfaces, so that inhomogeneous waves arriving at the .surface z=0 are significantly attenuated (recall that inhomogeneous plane waves have exponential attenuation). (c) The interpretation of the signature of inhomogeneous plane waves is not as straightforward as that of homogeneous plane waves (see 1-12). Finally, it seems that there is not general agreement on how much contribution the inhomogeneous plane waves contribute to the SWF (Carter, 1975). The inverse transform constructions obtained by the B-G variants employ kr as the independent variable and hence, for equi-angularly spaced plane-wave seismograms, the considered observations are not evenly spaced. On the other hand, upon change of variable from kr to y in equation (1-19), the Hankel inverse yields a relation for which equally spaced angular observations are considered. For the Backus and Gilbert model constructions the examples shown have standard deviation values a- defined as percentages of the maximum amplitude spectral value in the input spherical-wave seismograms (forward transformation) or in the input plane-wave seismograms (inverse transformation). Lastly, the phases appearing in the synthetic seismograms are marked by a series of letters denoting the modes in which a 34 particular arrival travelled through the different layers. For example, the sequence PSSP corresponds to a ray which left the source and travelled through the first layer as a P wave, converted to S mode at the first interface and then,travelled both down and up in the second layer as an S mode; finally, it was transmitted back into the first layer as a P wave, the mode recorded by the receiver. 3.2 Separation of Converted S Modes The separation of converted S modes from other modes in common mid-point gathers acquired in a marine environment has been discussed in recent publications by Ryu (1982), and Tatham et al. (1983). In both papers this separation is achieved by using the relatively higher moveout (compared to primary P arrivals) associated with the converted S modes. In this section we will follow the method discussed by Tatham et al. (1983) and will use the construction algorithms presented previously to separate the converted S modes from other modes in common-source gathers. For the purpose of this example we have used the reflectivity method (Fuchs and Muller, 1971) to generate a set of synthetic seismograms. These seismograms correspond to the model shown in Figure 3.1 with geophone offsets ranging from 100 m to 4000 m at 100 m intervals, time sampling interval 8 ms, and compressional-source function represented by a Ricker wavelet with a centre frequency at 16 Hz. For simplicity we have 35 limited our attention to synthetic seismograms which include primaries, head waves and converted modes. MODEL USED FOR GENERATING THE SYNTHETIC SEISMOGRAMS BY THE REFLECTIVITY METHOD SOURCE 500. HORIZONTAL RANGE (METERS) 1000. 1500. 2000. 2500. 3000. 3500, IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII Vp,= 1500 M/5 "si 866 M/S VP2-E500 M/S TS2 1443 M/S VP3=6500 M/S Vs3 = 3753 M/S (HALF SPACE) 200. 400. .. 600. .. 800. m m —i m UJ CO FIGURE 3.1 Homogeneous layer model displaying the geophone layout used to generate the synthetic seismograms (geophone spacing is 100m). Density was assumed constant throughout the model, whereas high velocity contrasts were used to enhance mode conversion. 36 Figure 3.2a shows those phases arriving at vertical-component receivers as P waves, whereas Figure 3.2b displays those phases arriving as S modes. Trace normalization was used to clearly show individual phases. The true amplitude seismograms corresponding to Figure 3.2a were added to those of Figure 3.2b and the result, which represents our field observations, is plotted in Figure 3.3. Our model corresponds to land data and hence it includes a number of converted modes which would not appear in the marine case. In separating these modes from the others in a CSG/ we distinguish two groups. In the first, we include all converted S phases which arrive at the receivers as P waves. This group includes those phases described in Ryu (1982), and Tatham et al. (1983) plus some other phases which may be separable on the basis of their moveout (i.e. phases exhibiting a large moveout will be mapped into the part of the plane-wave domain associated with large emergence angles and vice versa; see Figure 3.4). The second group includes converted modes which arrive at the receivers as S modes. These phases are separable in the plane-wave domain if the decomposition is carried out with a surface-layer velocity V, corresponding to the compressional modes. This separation is possible because the phases included in this group arrive at the surface with a true ray parameter sin^M/V, which is mapped into sin(7)/V^. Therefore we have sin(7)=V,fsin(75 )/Vs, , which suggests that converted S modes of the second group will be mapped to emergence angles which are considerably larger than their true emergence angles. In order to aid the identification 37 P-P SEISMOGRAMS TRACE NUMBER 15 10 15 20 25 30 35 40 Q . 0 ET 0.2 5L IX 0.4 3.2 3.4 3.6 3.8 FIGURE 3.2a Synthetic seismograms (normalized by trace) corresponding to the modes which left the source as P waves and arrived at the receivers as P waves. Three reflected phases are marked on this figure; head waves have not been labelled. 38 P-S SEISMOGRAMS TRACE NUMBER 3. 14 3.6 3.8 14. 0 FIGURE 3.2b Synthetic seismograms (normalized by trace) corresponding to the modes which left the source as P waves and arrived at the receivers as S waves. Three reflected phases are marked on this figure; head waves have not been labelled. 39 PP+PS SEISMOGRAMS (GLOBAL NORMALIZATION) 10 TRACE NUMBER 15 20 25 30 35 7 I -I | L'1 j. 1 40 . 0 . 2 . 4 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 m n FIGURE 3.3 Synthetic vertical component seismograms (with true relative amplitudes) corresponding to the model shown in Figure 3.1 with a point source of compressional waves. The temporal Fourier transforms of these seismograms represent the input data to PWD. 40 FIGURE 3.4 Large and small moveout signals. At first approximation, the PWD response from a given spherical-wave signal is governed by those paths given by geometrical optics. This is a good approximation for krr>>1 (i.e. for high frequencies and/or large offsets; Aki and Richards, 1980, ch.6). In (a) close receivers contribute to form plane-wave signal at small angles of emergence whereas far receivers do it for large emergence angles. Notice the differences in density of angular information for small and large offsets. In (b) the spherical-wave signal from a deep interface has small moveout throughout the receivers and hence, its plane-wave si-gna-1 will be observed at small angles of emergence. 41 of the observed modes in the plane-wave domain, we display in Figure 3.5 the theoretical trajectories of the modes which are expected for the given model. To start this example, we show the Hankel plane-wave seismograms in Figure 3.6. The converted S modes PPSS, PSSP, PSSS and PS are seen quite clearly in the area 0.2S<T<1.0S and 32°<7<88°. This zone is plotted again in Figure 3.7 for the cases of a) Hankel construction, b) smallest model construction, and c) flattest model construction. In cases b) and c) we have carried out the construction using a^=lO%, o^/til =20%, a8 = 15%, a, =12%, 0IO = 8% , and a- = 5% for the rest of the seismograms (the first values were set relatively high to reduce the diffraction effect associated with the high amplitudes at traces 1-11, t^0.5 s of Figure 3.3). Additionally the constant b in both K0(bkr) and K,(bkr) was set to 5. Comparison of the three panels in Figure 3.7 shows slight differences in event continuity and consistency in the zone of interest. These are attributed to the incorporation of the standard deviations and the weighting functions in the B-G variants. However, the general features in all the constructed models are quite similar. To effect the mode separation and compare the results in a domain in which most of us are somewhat more comfortable, we have applied the inverse transforms to the data sets shown in Figure 3.7. The reconstructed converted S modes obtained by the application of the inverse Hankel transform to the data of Figure 3.7a are shown in Figure 3.8, whereas those reconstructed 42 PREDICTED TIMES FOR PWD COMPONENTS PPPP PPSS ANGLE OF EMERGENCE (DEGREES) 2. 10. 20. 30. 40. 50. 60. 70. 3p , i PSSS 80. 90 p. Q fi™yffl4pssp nrff " - ^ 1.0 1. 5 2.0 2.5 3.0 3.5 4.0 m n FIGURE 3.5 Theoretical trajectories in the plane-wave domain of the modes appearing in Figure 3.3 and identified in Figures 3.2a and 3.2b. These times were computed by the delay time expressions, as those presented in equations (1-15), for the different wave-mode combinations. The vertical axis represents the delay time r. 43 PWS HANKEL TRANSFORM ANGLE OF EMERGENCE (DEGREES). 10 20 30 40 5,0 60 70 8.0 . o . 2 . 14 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 . 2 . 4 . 6 . 8 . 0 m n FIGURE 3.6 Plane-wave seismograms (normalized by trace) corresponding to the data of Figure 3.3. The decomposition was performed -using -the forward Hankel algorithm. The concentration of converted S energy in a distinct portion of the r-y domain is clearly visible. 44 PWS CONVERTED 5 ZONE ANGLE OF EMERGENCE (DEGREES) 2 10 20 3.0 40 50 60 70 80 QQ 0.2 0.4 0.6 0.8 1 . 0 1 . 2 psss PWS CONVERTED S ZONE 2 10 20 30 40 50 60 70 80 TT PS wm/G&gmn&Bs&Ea 0.0 0.2 0.4 0.6 0.8 1 . 0 1 . 2 m n PWS CONVERTED S ZONE 20 30 40 50 60 70 80 FIGURE 3.7 Portion of the plane-wave seismograms containing converted S modes. The decomposition was performed using (a) the forward -Hankel algorithm, (b) the forward smallest-model algorithm, and (c) the forward flattest-model algorithm. The seismograms are normalized by trace. 45 via the smallest model (with the data of Figure 3.7b) and the flattest model (with the data of Figure 3.7c) are shown in Figures 3.9 and 3.10, respectively (both reconstructions had a,-=2% for all i's and c=1). Although all the construction schemes were applied to exactly the same portion of the plane-wave domain, some significant differences in reconstruction quality are observed (compare the converted S modes of Figures 3.8 to 3.10 to those of Figures 3.2a and 3.2b). The converted S modes obtained via the Hankel forward followed by Hankel inverse, exhibit a rather poor signal at large offsets. However both B-G reconstructed t-x seismograms present reasonable results, except for a strong diffraction at near offset and small times which is related to the truncation of the transformed zones in the plane-wave domain. We attribute the difference in reconstruction quality of Figures 3.8 to 3.10 to the following factors: a) the relative attenuation of higher . frequencies by the weighting functions K0(busin(y)/V) and K,(bwsin(7)/V) may reduce aliasing effects; b) the ability of the B-G variants to incorporate noise considerations into the reconstruction adds numerical stability; and c) poor angular sampling coupled with a bad choice of integration rule prevents reliable t-x reconstruction using the Hankel inverse. To further illustrate these points we have reconstructed the t-x seismograms via the flattest model with the data of -Figure 3.7a as input (i.e. Hankel transform forward and flattest model inverse). The result is shown in Figure 3.11. Again the improvement over the Hankel forward - Hankel inverse 46 CONVERTED S MODES HANKEL-HANKEL TRACE NUMBER FIGURE 3.8 Converted S-mode t-x seismograms (normalized by trace) corresponding to the data shown in Figure 3.7a. The reconstruction was performed using the inverse Hankel algorithm. Notice the deterioration of all events at long offsets (high trace numbers). CONVERTED S MODES SMALLEST-SMALLEST FIGURE 3.9 Converted S-mode t-x seismograms (normalized by trace) corresponding to the data shown in Figure 3.7b. The reconstruction was performed using the -inverse smallest-model algorithm. 48 CONVERTED S MODES FLATTEST-FLATTEST FIGURE 3.10 Converted S-mode t-x seismograms (normalized by trace) corresponding to the data shown in Figure 3.7c. The reconstruction was performed using the inverse flattest-model algorithm. 49 (Figure 3.8) is significant; the reconstruction differs only slightly from that of the flattest forward-flattest inverse (Figure 3.10), the clearest difference being evident for the event PPSS. In order to familiarize the reader with the effects of applying different standard deviation weights, we have reconstructed the converted S-mode t-x seismograms via the smallest model with the data of Figure 3.7b as input, setting a,-=l0% for all i's. The result is shown in Figure 3.12. Comparison of this result with that of Figure 3.9 (the models of these figures differ only in the standard deviation values used in the smallest inverse) reveals that these changes in a; did not affect significantly the reconstruction. However, a detailed comparison of wavelet definition in events PS and PSSS (particularly at large offsets) shows a loss of frequency band width for the case represented by Figure 3.12. This loss is due mainly to two factors: a) the inversion at any given frequency with a large o; value results in a model which is constructed with a smaller number of basis functions; and b) for the multiple set of problems to be solved in PWD, our choice of a; as a percentage of the maximum value of the amplitude spectrum of the i'th trace, may cause further spectral discrepancies. To clarify the latter point, suppose that the maximum amplitude spectral value in each of the given seismograms is located at the same frequency, say f0. Observational errors at other frequencies are relatively larger than the errors at f0. 50 CONVERTED S ZONE HANKEL-FLflTTEST TRACE NUMBER 15 10 15 20 25 30 35 HO 3. 14 3 . 6 3 . 8 U . 0 FIGURE 3.11 Converted S-mode t-x seismograms (normalized by trace) corresponding to the data shown in Figure 3.7a. The reconstruction was performed using the inverse flattest-model algorithm. Compare this result to the one shown in Figure 3.8. 51 CONVERTED S ZONE SMALLEST-SMALLEST FIGURE 3.12 Converted S-mode t-x seismograms (normalized by trace) corresponding to the data shown in Figure 3.7b. The reconstruction was performed using the inverse smallest-model algorithm. In this case, we assigned a;=10%- of the maximum spectral amplitude value of each (input) plane-wave seismogram. 52 This causes the expected x2 values for the models corresponding to these frequencies to be reached with a smaller number of basis functions (smaller with respect to the number of observations). In particular, at those frequencies with low energy (in comparison to the largest amplitude spectral peak), models are constructed with only a few basis functions and this may result in an even larger amplitude discrepancy (i.e. narrower band). Further description of problems associated with standard deviation weights for the plane-wave decomposition problem is given in Appendix B. 3.3 White Noise and Construction Quality For this part of the presentation we have added a uniformly distributed random noise series to each of the seismograms shown in Figure 3.3. These noise series are limited to the closed region [-0.5X**'"' ,+0.5Xp,v ], where X1*** is the maximum absolute value of the i'th seismogram. The noisy seismograms (normalized by trace) are shown in Figure 3.13. We have reconstructed the t-x representation of the converted S modes using the B-G variants smallest(forward)-smallest(inverse), and flattest(forward)-flattest(inverse), with the same o;, b and c values as in the previous section. These results are shown in Figures 3.14 and 3.15 respectively. The reconstructions closely resemble those of Figures 3.9 and 3.10 for which no additive noise was present. However there is a 53 PP+PS SEISMOGRAMS ADDITIVE RANDOM NOISE 10 TRACE NUMBER 15 20 25 30 35 HO ^?2»llil!iIf41ifm^ffI4i^i^;i tmmiimimmmmmummmiM mmmmmimmmmwummn mmmmmmmmmmimmim mmmimimimmmuiimmmi mmmmuimmmmtmmmiu mmmmmmmmmmm'Miun tmmmmmmmmmmmmm m\mmmmmmmmmmmsm mmmmmm&mmmmiimm tmmmmmmmmmmimixrm ummimmnmmmmmmmm mnnnmmmmmvimiumiim nmmmmmmimmmitmmm 0 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 3 4 0 2 4 6 8 0 2 4 6 1 . 8 2.0 2 4 6 8 0 2 4 6 8 0 in m n FIGURE 3.13 Noisy t-x seismograms (normalized by trace) generated by the addition of random noise series to the seismograms shown in Figure 3.3. 54 CONVERTED S MODES SMALLEST-SMALLEST (NOISY SEISMOGRAMS) TRACE NUMBER 3. 14 3 . 6 3 . 8 14. 0 FIGURE 3.14 Converted S-mode t-x seismograms obtained by (i) PWD of Figure 3.13 using the forward smallest-model algorithm, and (ii) reconstruction of the t-x seismograms using the inverse smallest-model algorithm. The seismograms are normalized by trace. 55 CONVERTED S MODES FLRTTEST-FLATTEST (NOISY SEISMOGRAMS) TRACE NUMBER 3 . 4 3 . 6 3, 8 4.0 FIGURE 3.15 Converted S-mode t-x seismograms obtained by (i) PWD of Figure 3.13 using the forward flattest-model algorithm, and (ii) reconstruction of the t-x seismograms using the inverse flattest-model algorithm. The seismograms are normalized by trace. 56 minor deterioration of the event PPSS at medium and long offsets, whereas events PSSS, PSSP and PS have been reproduced quite faithfully. We mention that the reconstructed t-x converted S-mode seismograms corresponding to the Hankel(forward)-Hankel(inverse) algorithms are very similar to those reconstructed seismograms using the Hankel(forward)-Hankel(inverse) algorithms for the case of no additive noise. Our experience to date shows that white noise does not severely deteriorate the performance of the algorithms presented in this work in regions where the angular data are reasonably dense, for example, at long offsets for shallow horizontally layered earth models. 3.4 Angular Sampling and Trace Interpolation This section is divided into two parts. In the first, we deal with the problem which arises when information concerning a shallow zone of interest is recorded by an evenly spaced geophone array. The second section concerns the problem of trace interpolation, i.e. a number of gaps exist in a set of otherwise evenly spaced seismograms; we will use the algorithms developed here in order to fill these gaps. 57 3.4.1 Shallow zone of interest Geophone arrays consisting of evenly spaced receivers with spacing large as compared to the depth of the first layer of interest (horizontally layered earth model), will tend to have poor information density in the portion of the plane-wave domain corresponding to small angles of incidence (see Figure 3.4). On the other hand, a fairly large spacing will still yield the desired information density when our interest is confined to the portion of the plane-wave domain corresponding to large angles. Consequently, an array of closely spaced geophones at the near offsets, followed by widely spaced geophones at the far offsets could satisfy our information density requirements, at least kinematically. In Figure 3.16 we show a set of seismograms (normalized by trace) corresponding to the model of Figure 3.1. The first twenty traces are evenly spaced with a geophone separation of 50 m, whereas the remaining forty are 100 m apart. The plane-wave components of this set (normalized by trace) are shown in Figures 3.17 and 3.18 for the smallest and flattest models respectively. As is expected, both constructed models show continuity for that part of the PPPP event which was buried inside the diffraction event of Figure 3.6 (this diffraction was caused by the relatively large amplitudes in the t-x domain, traces 1 to 11 between 0.3s and 0.8s. Such large amplitudes at the end of t-x events give significant contribution of energy at all angles of emergence). Also, the PP phase shows more 58 PP+PS SEISMOGRAMS (NEAR-OFFSET RESAMPLING) TRACE NUMBER •3,2 •3.4 •3.6 •3.8 ; 4 . 0 FIGURE 3.16 t-x seismograms (normalized by trace) corresponding to the model of Figure 3.1 . Receiver spacing is 50m for traces 1-20 and 100m for traces 21-50. Various phases are identified. 59 PHS SMALLEST MODEL (NEAR-OFFSET RESAMPLING) ANGLE OF EMERGENCE (DEGREES) 3.4 3 . 6 3 . 8 4. 0 FIGURE 3.17 Plane-wave seismograms (normalized by trace) corresponding to the data shown in Figure 3.16. The decomposition was performed using the forward smallest-model algorithm. Compare this result to that of Figure 3.6 and notice the continuity of event PPPP for small angles of emergence. 60 PWS FLATTEST MODEL (NEAR-OFFSET RESAMPLING) ANGLE OF EMERGENCE (DEGREES) 2. 10. 20. 30. UO. 50. 60. 70. 80. 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 1. 2 1. 4 1. 6 1. 8 i 2. 0 m 2. 2 m n 2 . 4 2 . 6 2. 8 3. 0 3. 2 3. 4 3 . 6 3 . 8 4 . 0 FIGURE 3.18 Plane-wave seismograms (normalized by trace) corresponding to the data shown in Figure 3.16. The decomposition was performed using the forward flattest-model algorithm. The results are very similar to Figure 3.17. 61 prominently at the smaller angles of emergence. In regards to the converted S modes, the PS-phase signature is expected to have better continuity at smaller angles of emergence but because of the normalization used, Figure 3.17 does not show such an improvement (the PP phase dominates the amplitudes). The remaining converted S modes are not expected to improve since their spherical-wave signature appears at large offsets (see Figure 3.16). 3.4.2 Trace interpolation The problem of trace interpolation viewed in the context of plane-wave decomposition is heavily dependent on the information density in the angular domain. Assuming that a certain portion of the plane-wave domain is well sampled, using our algorithms with a reasonable b value (to increase linear dependence of the kernels), we should be able to interpolate seismograms between geophone locations for the phases which belong to such a portion. The converted S modes for the model of Figure 3.1 afford a reasonable example of this operation. For the purpose of this example we have eliminated traces 6,12,18,24,36 and 37 from the set of seismograms of Figure 3.3 and plotted the result (normalized by trace) in Figure 3.19. We constructed the t-x converted S-mode seismograms using the B-G variants smallest(forward) - smallest(inverse) and flattest(forward) - flattest(inverse). The results (normalized 62 by trace) are shown in Figures 3.20 and 3.21, respectively. Both construction schemes succeed in the interpolation of the modes PSSS, PSSP and PS. However both do rather poorly with respect to the mode PPSS. This is because the mode PPSS has its angular signature in the badly sampled zone corresponding to small angles of emergence.•The small loss of frequency band-width observed in the results obtained is due to the respective weighting function (for details please refer to Appendix B). 63 PP+PS SEISMOGRAMS (UNEVEN SPACING) TRACE NUMBER 10 15 20 25 30 35 110 m O FIGURE 3.19 t-x seismograms (normalized by trace) corresponding to the model of Figure 3.1. It is assumed that due to some data acquisition difficulties, traces 6, 12, 18, 24, 30, 36 and 37 are dead. INTERPOLATED S MODES SMALLEST-SMALLEST FIGURE 3.20 Converted S-mode t-x seismograms obtained by (i) PWD of Figure 3.19 using the forward smallest-model algorithm, and (ii) reconstruction of the t-x seismograms using the inverse smallest-model algorithm. The seismograms are normalized by trace. 65 INTERPOLATED S MODES FLATTEST-FLATTEST TRACE NUMBER 1 5 10 15 20 2 5 30 35 40 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 1. 2 1. 4 1. 6 —i 1. 8 i 2. 0 m </> 2. 2 m n 2. 4 2. 6 2. 8 3. 0 3. 2 3. 4 3 . 6 3 . 8 4. 0 FIGURE 3.21 Converted S-mode t-x seismograms obtained by (i) PWD of Figure 3.19 using the forward smallest-model algorithm, and (ii) reconstruction of the t-x seismograms using the inverse flattest-model algorithm. The seismograms are normalized by trace. 66 4. COMPUTATIONAL CONSIDERATIONS The CPU times (on an Amdahl 470 V/8 computer) for the PWD algorithms are presented in Table II below. The seismograms used to generate this table were digitized at 8 ms and were limited to the band . 5-25 Hz (the number of frequency samples in this band is equivalent to those contained in the band 10-50 Hz on data with 4 ms sampling interval). Furthermore, all calculations (with the exception of Bessel function computations) were made in double precision. TABLE II TIME INPUT OUTPUT CPU TIMES SAMPLES TRACES TRACES HANKEL SMALLEST FLATTEST 1 28 256 512 . 40 40 40 45 45 45 2.416s 4.768s 9.349s 3 .556s 6.770s 12.965s 5.342s 10.192s 19.909s The following points are emphasized: (i) Most of the processing time is elapsed in Bessel function computations. For these, we have used polynomial approximations with errors of 0(10"8) given in Abramowitz and Stegun (1970). Efficiency of these computations can be increased by: (a) decreasing the accuracy of the approximations which will reduce processing times, though it 67 may yield some degradation of the output, and (b) using the derivative relation between J0 and J, in flattest model calculat ions. From the summation represented by equation (2-10) for the smallest or that of equation (2-15) for the flattest model, it is deduced that increasing the number of input or output traces causes a linear increase in the number of Bessel function computations. Similarly, an increase in the number of time or frequency samples results in a linear increase in the number of Bessel function computations. (ii) The CPU time for the singular value decomposition (SVD) of the inner product matrix T behaves like the cube of the number of input traces. Hence, depending on the ratio between the number of time and frequency samples to the number of input traces, a significant percentage of the total run time may be spent in SVD. However, when processing a large number of CSG's with fixed geometry and standard deviation estimates, SVD is executed once. In this case, the total smallest-model CPU time is comparable to that of the Hankel algorithm. (iii) The memory requirements of the B-G algorithms include one array of size NPTSxNTRACE and two of size NTRACExNTRACE, where NTRACE is the number of either input or output traces whichever is larger, and NPTS is the number of samples per trace. In contrast, the Hankel algorithm needs only one array of size NPTSxNTRACE. 68 (iv) The CPU run times and memory requirements corresponding to the inverse transform algorithms are equivalent to those of the forward transform. 69 5. SUMMARY Efficient algorithms for the decomposition of a spherical-wave field into its plane-wave components have been presented. Also, it has been shown that these algorithms allow the reconstruction of the spherical-wave field from a relatively sparse sample of its plane-wave components. The practical viability of the proposed algorithms has been demonstrated using the problem of separation of converted S modes from other modes in a common-source gather, and the problem of trace interpolation. The following points should be noted: (i) The plane-wave seismograms are obtained by using the Backus & Gilbert construction techniques, subjected to the requirement of weighted smallest or flattest model. (ii) The construction schemes allow the handling of errors in the data and hence, permit a certain control on the model structure provided by the basis functions. Caution should be exercised in assigning the standard deviation values to unnormalized observations. (iii) The proposed algorithms are not limited to evenly spaced data and consequently, allow the design of an appropriate geophone array which should produce a more faithful representation of the plane-wave components. 70 (iv) Numerical stability is gained by proper use of the b value appearing in the arguments of the weighting functions. Large b values decrease the degree of linear independence of the kernels and seem to be appropriate for the problem of trace interpolation. (v) Dynamic aspects (e.g: true amplitudes) of the forward and inverse contructed models are not formally handled by the algorithms as developed. (vi) For a given angle of emergence 7, homogeneous plane waves corresponding to high frequencies may be significantly attenuated by the weighting function Q. This function then, is viewed as a potential aliasing suppressor. (vii) For a given angular frequency CJ, large wavenumber components are severely attenuated by the weighting function Q. Indeed, depending on the b value chosen, inhomogeneous waves associated with wavenumbers larger than a certain value are practically excluded from the decomposition. This effect is analogous to formulating PWD as an inverse problem with finite limits of integration. 71 BIBLIOGRAPHY Abramowitz, M. , and Stegun, I., 1970, Handbook of Mathematical Functions: New York, Dover Publications Inc., 1046 p. Aki, K., and Richards, P. G., 1980, Quantitative Seismology. Theory and Methods: San Francisco, W. H. Freeman and Co., v. I , 557 p. Asby, R., and Wolf, E., 1971, Evanescent waves and the' electromagnetic field of a moving charged particle: J. Opt. Soc. Am., v. 61, p. 52-59. Backus, G., 1970, Inference from inadequate and inaccurate data, II: Proc. Nat. Acad, of Sci., v. 65,- p. 281 -287. Bath, M., 1968, Mathematical Aspects of Seismology: Amsterdam, Elsevier Publishing Co., 415 p. Born, M., and Wolf, E., 1980, Principles of Optics: Toronto, Pergamon Press., 808 p. Brekhovskikh, L. M., 1960, Waves in Layered Media: New York, Academic Press., 561 p. Carter, W. H., 1975, Band-limited angular-spectrum approximation to a spherical scalar wave field: J. Opt. Soc. Am., v. 65, p. 1054-1058. Devaney, A. J., and Sherman, G. C, 1973, Plane-wave representations for scalar wave fields: SIAM Rev., v.15 , p. 765-786. Fuchs, K., 1968, The reflection of spherical waves from transition zones with arbitrary depth-dependent elastic moduli and density: J. Phys. of the Earth, v.16, p. 27-41. Fuchs, K., and Muller, G., 1971, Computation of synthetic seismograms with the reflectivity method and comparision with observations: Geophys. J. R. Astr. Soc, v. 23, p. 417-433. Goodman, J. W. , 1968, Introduction to Fourier Optics: San Francisco, McGraw-Hill Co., 287 p. Gradshteyn, I., and Ryzhik, I., 1980, Table of Integrals, Series and Products: Toronto, Academic Press, 1160 p. Henry, M., Orcutt, J., and Parker, R., 1980, A new method for slant stacking refraction data: Geophys. Res. Lett., v. 7, p. 1073-1076. 72 Muller, G., 1971, Direct inversion of seismic observations: Zeitschrift fur Geophysik, v. 37, p. 225-235. Oldenburg, D. W, and Samson, J. C, 1979, Inversion of interferometric data from cylindrically symmetric refractionless plasmas: J. Opt. Soc. Am., v. 69, p. 927-941. Parker, R. L., 1977, Understanding inverse theory: Ann. Rev. Earth Planet. Sci., v. 5, p. 35-64. Ryu, J. V., 1982, Decomposition (DECOM) approach applied to wave field analysis with seismic reflection records: Geophysics, v. 47, p. 869-883. Schultz, P. S., and Claerbout, J. F., 1978, Velocity estimation and downward continuation by wavefront synthesis: Geophysics, v. 43, p. 691-714. Stratton, J. A., 1941, Electromagnetic Theory: New York, McGraw-Hill Co., 615 p. Tatham, R. H., Goolsbee, D. V., Marsel, W. F., and Nelson, H. R., 1983, Seismic shear-wave observations in a physical model experiment: Geophysics, v. 48, p. 688-701. Treitel, S., Gutowski, P., and Wagner, D., 1982, Plane-wave decomposition of seismograms: Geophysics, v. 47, p. 1375-1401. 73 APPENDIX A X2 Value and Observational Errors in Model Construction When solving problems which are associated with inaccurate observations, it is undesirable to construct models which reproduce these data exactly. In this case, it is common to require the calculated observations to fit the data in a manner consistent with the observational errors. In this appendix, we outline the steps required in the construction of models with calculated observations which are related to the observed data by an expected x2 value of approximately N. Firstly, from equation (2-5a) we have, Ta = 1° Expressing T in terms of its spectral components and solving for a we get, a = r-1t° = RA"1RTe° (A-1) where R is an (NxN) matrix whose columns are the eigenvectors of T, A is an (NxN) diagonal matrix whose diagonal consists of the eigenvalues of T arranged in decreasing order, and RT is the transpose of R. -Misfitting the observations e° is readily achieved by winnowing, say, the K smallest eigenvalues with their associated eigenvectors, that is, truncating matrices R, A and RT to size 74 (NxM), (MxM) and (MxN) respectively, with M=N-K. The coefficients ac constructed from the truncated set of spectral components are, a? = R^Rji0 (A-2)' Using ac we can compute the calculated observations, i.e. 1° = Tac (A-3) Secondly, consider the x2 value defined by X2 = Z(eT-eP)2 = ||A?||2 (A-4) with A? = e'-e° and E(x2}-N for N>5. The length of the vector Ae is not changed upon rotation. Therefore we can project e1 and e° onto the eigenvectors of R, i.e. X2 = ||RTt" - RTe°||2 or X2 = £(e? - S?)2 (A-5) Finally, from (A-3) and (A-2) we write, e6 = RARTRMA:1Rjt° Premultiplying this equation by RT and writing down the notation for rotated observations we obtain, ec = ARTRwA;'e0 , (A-6) A from which we realize that ef = e° for i=1,...,M and ef = 0 for i=M+1,...,N. Hence (A-5) gives, X2 =z(i?)2 75 For complex data we use, X2 = Ze^e?*" (A-7) where * indicates the complex conjugate. Starting with M=N-1 we form the summation given in equation (A-7), and keep adding terms until this summation yields the closest value to N (number of observations). The final index M gives the number of eigenvalues and eigenvectors to be retained in the calculation of ac . These coefficents, when used in the construction, yield a model which satisfies the observations in a manner consistent with observational errors. 76 APPENDIX B Effects of the Weighting Function Q and Standard Deviation Values on Model Construction After introducing the standard deviation values o- and weighting function Q, our original problem e,- =<G,- ,m> has been modified to, e; = <G; ,m* > (B-1) where, e; = e; / O ; G* = G;Q/a; m+ = Q"1m Q"1 = 1/Q with Q=k;°-5K00-5 (kr b) for the (forward) smallest problem Q=K,°'5(krb) for the (forward) flattest problem and b an arbitrary positive real number. In this appendix we highlight the effects of Q with a given b value, and of a- on the constructed model m=U(cj,kr). 1. Effects of Q on the size of m. The role of the constant b in attenuating large horizontal wavenumber components is portrayed in Figure B.1. This figure shows the plots of the modified Bessel functions K0(r(7)) and 77 K,(r(7)), where the argument is defined by r(7)=bwsin(7)/V with 0.1°<7<90°, cj=407rrad/s, V=1500 m/s , and b assumes the values 1 and 10. Clear differences on attenuation rate imposed by these functions are observed, such that the size of the constructed model will be significantly affected. In what follows we discuss the smallest model problem, whose weighting function involves K0 . 0.0 90.0 0.0 90.0 AXIS IN DEGREES PXIS IN DEGREES (a) (b) FIGURE B.I Modified Bessel functions (a) K0, and (b) K, for two different b values in the argument r(7) (r(7)=bwsin(y)/V, tj=407rrad/s and V=1500m/s). Large b values severely attenuate large horizontal wavenumbers and decrease the linear independence of the kernels. The minimization of ||m*|| requires that (Backus, 1970), | |Q- 'm| |<M (B-2) 78 where M is a finite real positive number. The size of m for large horizontal wavenumbers depends heavily on the rate of growth of Q" 1 . Using an asymptotic approximation to K0 for large arguments, Q~ 1 behaves as (kr )0-5 (bkr ) °-2 5exp(bkr/2) . Hence, the search for models m=U(ci>,kr) is confined to those whose high wavenumbers are strongly attenuated. This model attenuation is more severe for weighting functions Q with large b values. Homogeneous plane waves have horizontal wavenumbers restricted between 0 and CJ/V (i.e. angles of emergence between 0° and 90°). For typical exploration seismic work, CJ/V is smaller than unity and hence, with small b values, these wave components are not severely attenuated by Q"1 in the model m. On the other hand, inhomogeneous waves have wavenumbers between u/V and °° (i.e. complex angles between 90° and 90°-i<»: Brekhovskikh, 1960). When kr>u)/V but is "reasonably small", these waves are still controlled by the data equation (B-1). However, as kr becomes larger the inhomogeneous waves are increasingly attenuated by the requirement specified in (B-2) so that their amplitudes will decrease exponentially. It is clear that the "transition" value of kr at which constraint (B-2) predominates (B-1) depends on the chosen b value, that is to say | |kr°-5(krb)°-25exp(krb/2)U(cj,kr) | | has to be kept finite as kr becomes large (kr-*•»). The preceding discussion applies to the forward flattest problem as well. In this case, the weighting function Q affects the derivative of m, or equivalently, the rate of change of the 79 contribution of the homogeneous and inhomogeneous plane waves at a given angular frequency GJ. 2. Effects of Q on the inner product matrix T. The elements of the inner product matrix T for the (forward) smallest and flattest problems are inversely proportional to [ (r;2+b2 + r-2) 2-4r-2rJ2 ], where r; and rj are the offsets corresponding to the i'th and j'th geophone locations such that r- <rj for i<j. Consequently, in the absence of normalizaton by the standard deviation estimates, we have ri 1>T2 2>. . .>ruiJ . For b<<r„„j the diagonal elements of T dominate greatly over the off-diagonal ones. That is to say, the weighted kernels G-, Q are almost orthogonal. In this case, the matrix T is numerically well behaved (good condition number, see Figure B.2). Furthermore, the eigenvalues of T are approximately equal to its diagonal elements, and its eigenvectors are almost the unit vectors (specifically, the eigenvector matrix R is approximately equal to the identity matrix I). Therefore, the A w i' th rotated data e; =ER,.e; and the l'th basis function =[ 1 //A- ]?Rj;Gj Q become almost the i ' th observed data e; and the i'th scaled weighted kernel G;Q//A| respectively. An increasing b value will decrease the linear independence of the weighted kernels G;Q with the effect of making the diagonal elements of T less prominent with respect to the off-diagonal ones. In this case, the i'th rotated data e; and the i' th basis function \p- will be formed from a linear combination 80 of more than one element of the observed data e; weighted kernels G;Q respectively. and of the rr UJ m o o O FIGURE B.2 9. 0 2. 0 -3. 0 3. 0 LOG OF b (SMALLEST MODEL) 10 The condition number of the inner product matrix T (forward smallest model) as a function of b. Notice the wide domain of b values for which T is well behaved. To illustrate the previous comments, we have displayed in Figure B.3 some of the eigenvectors corresponding to the forward smallest problem for a set of distances r;=iAr, i=1,2,...,40, Ar=l00 m, with b=0.1 (Panel (a)) and b=50.0 (Panel (b)). In Panel (a) we observe that the eigenvectors are practically unit vectors and hence, each basis function will essentially correspond to a single, scaled weighted kernel. For example, the 25th basis function shown in Figure B.4a (b=0.l) resembles a 81 Panel a 0 . 1 E » 0 11 i—i—i—i—i— Panel b 0 . I E ' 0 1 i i i —i—i—L_I i—1_ -0 . I E » 0 1J l l l I l—l—I—1—1— 1 . ilO. EIGENVECTOR 1 0 . 1E *0 li i i i—i—i—i—i—u EIGENVECTOR 1 0 . 1 E * 0 11 i i—i —i—i—i i i i_ •0. IE "Oil i i I i I—I—i—I—i—I -O.lE-OlJ—I—I—i—i—i—I—i—i— 1 . 40. 1 . 40. EIGENVECTOR 13 EIGENVECTOR 13 0. IE'Oil i i—i—i—i— 0 . 1 E»0 11 c—i—i—i—i—i—i—L. - 0 . 1 E » 0 1J I I i—J—i—i—L—J—J—1Q -O.lE'OlJ—I—I—1_ 40. EIGENVECTOR 25 0 . 1 E •» 0 11_i i_i i i i—i i—i— -0. 1 E•»0 U I I I I—I—u 1 . EIGENVECTOR 25 0 . I E • 0 i i i i i i i—i— EIGENVECTOR 35 -0. 1E •» 01J i—i I—i—I—I—I—I—u_„ 40. 1. 40. EIGENVECTOR 35 FIGURE B.3 (HORIZONTAL AXES: ELEMENT OF EIGENVECTOR) Some of the eigenvectors corresponding to the inner product matrix of the forward smallest-model algorithm. The parameter b is set to (a) b=0.l0, and (b) b=50.0. By assigning a constant value o,-=1 for all i's, the inner product matrix reflects only the geometry of the geophone array. Notice that the large b value has "dispersed" the structure of the eigenvectors. 82 single (zero-order) Bessel function of the first kind with the attenuation specified by the modified Bessel function K0(krb). On the other hand, in Panel (b) of Figure B.3 we observe that the large b value has "dispersed" the eigenvector structures. In this case, then, the 25th basis function (shown in Figure B.4b) contains contributions from a number of scaled weighted kernels, specifically from the 12th to the 40th inclusively. 10 HERTZ 0 . 6 E + 0 2 L__U_L_—i : i i i i—I 0.1E + 03|——J 1 1 1 1 1 L. -0 . 8E + 02. (HORIZONTAL AXES IN DEGREES) FIGURE B.4 The 25th basis function for the forward smallest-model contruction with (a) b=0.l0, and (b) b=50.0. Referring to the 25th eigenvector displayed in Figure B.3, notice that the 25th basis function in (a) is practically a single kernel (namely the 25th), whereas this basis function in (b) is effectively a linear combination of several kernels. 83 3. Effect of a- on the inner product matrix T. In section 2. above we have shown that neglecting standard deviations leads to an inner product matrix Y in which T,,>T22>...>C„. This behaviour of the diagonal elements of T provides a natural ordering for the basis functions such that the long-wavelength structure of the model is controlled by the basis functions associated with the largest eigenvalues, whereas the fine structure originates from those basis functions associated. with the smallest eigenvalues. Dividing each element of T by a- <7j the relations , >T2 2> . . . >rufJ do not necessarily hold. In particular, for the case of small b value and largely varying standard deviations, we do not expect that the eigenvalues of T arranged in decreasing order will correspond to its diagonal elements in their original order. The new order relations will depend on the relationships between the standard deviations. For example, a very small o2 value could make rv>J the largest diagonal element and consequently, the first ordered eigenvalue will essentially correspond to this element. But the most important result of this reordering of eigenvalues versus diagonal elements of T is the consequent reordering of the basis functions i>; . Hence, in the above example, the last weighted kernel (associated with YKtl) could become the first basis function and consequently small-wavelength model structure will stem from This is illustrated in Figure B.5 where we have plotted on Panel (a) some basis functions for the case a,=1, whereas on Panel (b) we have displayed the same basis functions with O,=10% of the maximum (amplitude) spectral value of the 84 i'th trace of the data presented in Figure 3. We summarize by stating that in the process of winnowing basis functions, the standard deviation values will play a major role in determining the type of information to be included in the constructed model. Thus application of a-, values to a given PWD problem should be exercised with caution. 85 Panel a 13 HER7Z Panel b 0. 3E'0 2t_j__L_i_(_L i—i—L—i O.UE-OSt -0. 2E-0 3L_i_i—I I 1 1—J—I—I—I -0.5E»03J—i—I—I—I—I—i—J—I—u 2 . 90. 2. 90. p pf-.ni, , , , . i i i i i 0.2E«03i—i—i—i—i—i 0. 3E->03L_i—j i i i i i i i r 0.2E*03i^-i i—I—i—i—I—i—i—t 25 - 0. 3 E • 0 3J1 I I I I I—I—1_ 2. 0. 2E»03i i I i—i—i—i—i -0. 2E-»03J 25 _i I -0.3E*031 I 1 I I—I—I—I—1— 90. 2. 90. 0. UE*03i i I i i i i—i—I—L _ -0.1E*03J-90. 2. (HORIZONTAL RXES IN DEGREES! FIGURE B.5 Basis functions 1, 13, 25 and 35 corresponding to the forward smallest-model construction. Standard deviations are set to (a) 1 for all input traces, and (b) percentages of the maximum spectral amplitude of each trace (see text for details). Notice that the re-ordering of the basis functions is such that in (a) the first basis function contributes with long-wavelength model structure, whereas in (b) this basis function gives short-wavelength structure. 86 APPENDIX C Inner Product Matrix for the Forward Smallest Model Construction The construction of the forward smallest model requires the evaluation of the elements of the inner product matrix X. Using the kernels from equation (2-9b) we have, r:. = <GiQ,GjQ> = [l/a,-o;] fkrK0(krb)J0(kr r; )J0(kT rj )dkr (C-1) 0 and from Gradshteyn and Ryzhik, 1980, equation 6.578.15 we obta in, r;- = [ l/a; a-]/[ (r?+b2 + ri2)2-4ri2rj2]°-5 (C-2) 87 APPENDIX D The Forward Flattest Model: Global Development This appendix describes the solutions to the set of problems which are encountered in the construction of the flattest-model forward transform. In particular, the following problems are undertaken: 1. Computation of the flattest model kernels H;. 2. Evaluation of the term H; m°f . 3. Calculation of the inner product <R,Q,HjQ> Kr 4. Solution to the indefinite integral JQ2H;. 5. Evaluation of the integration constant of equation (2-14). 1. Computation of the flattest model kernels H;. With the identification, = kr J0 iK rr )/o: we evaluate the indefinite integral, H; = [ ]/k'J0(k'r; )dk' Upon change of variable r' = k'r; we get, H; = [\/a. r,2]/r'J0(r')dr' and from Gradshteyn and Ryzhik, 1980, equation 5.56.2 we obtain, H; = [l/a; r; ]krJ,(krr; ) (D-1) 88 2. Evaluation of the term H;mj . Using (D-1) and the asymptotic behaviour of J,(krrr ) for large argument, the first term of the r.h.s. of equation (2-11) for kr-»°° behaves as, U(oj,kr-°°) lim {[ 1/ff, r,- ]kr[2Akrr, ] ^ } = U(u,kr») lim {2kr,,2/ff, 7r1'2r;3/2} (D-2) This expression goes to zero provided U(a>,kr—°°) goes to zero faster than kr"1/2 for all frequencies (because the modified Bessel functions K0(x) and K,(x) are not defined at x=0 we exclude both o>=0 and 7=0°). Recalling that kr=wsin (y)/V, and considering a band-limited source, kr-*°° represents infinitely-attenuated inhomogeneous waves. In the present work we set • H; mT =0 Then we will be concerned with models U(w,kr) which go to zero faster than kr~1'2 as kr-»°°. 3. Calculation of the inner product matrix. With Hj as given in (D-1) and Q = K,°-5(krb), the entries of the inner product matrix T are given by, r;- = <H,Q,HJQ> = [ 1/r- r- o.o- ]|kr2K, (krb)J, (krr. )J, (krrj )dkr and from Gradshteyn and Ryzhik, 1980, equation 6.578.15 we have, ri-=[\/oios] 4b/[(r2+b2 + r.2)2-4r2r2]1-5 (D-3) 89 Kr 4. Evaluation of the indefinite integral JQ2H-. The solution to the indefinite integral appearing in the l.h.s. of equation (2-15) can be derived from the following (indefinite) integral given in Gradshteyn and Ryzhik, 1980, equation 5.54.1, JyZp(ay)BP(cy)dy = {cyZ ?(ay )B*_, (cy) - ay Z?_, (ay ) BP (cy)}/[ a2-c 2 ] (D-4) where Zv and B? are any Bessel function of order p, and a and c are constants. For p=1 we identify, Z, = J, B, = H,"' H,(1)is the first Hankel function of order one. (D-4) now reads as, JyJ,(ay)H,'1'(cy)dy = {cyJ,(ay)H0< 1'(cy) - ayJ0(ay)H,(1'(cy)}/[a2-c2] (D-5) but from Gradshteyn and Ryzhik, 1980, equations 8.407.1 and 8.407.2 we have, H0(1)(iby) = -[2/TT] iK0(by) (D-6a) H,(1»(iby) = -[2/TT] K,(by) • (D-6b) Using (D-6) and identifying c=ib, (D-5) transforms to, -[2/TT] JyJ, (ay)K, (by)dy = {(ib)yJ,(ay)[-2/ir iK0(by)] - ayJ0(ay)[-2/w K,(by)]} /[a2-(ib)2] 90 that is to say, fyj,(ay)K,(by)dy = - {byJ,(ay)K0(by) + ayJ0(ay)K,(by)}/[a2+b2] (D-7) from which after making y = kr and a = r-( we obtain the r.h.s. of equation (2-15). 5. The integration constant of equation (2-14). We show that with the condition U(u, kr-»°°) = 0 the constant of integration C is equal to zero. From equations (2-14), (D-1) and (D-7) (with y = kr and a = r,- ) we see that, C = -Z[0,-/a;rj ] [l/(r?+b2)]{ 1 im{ bkr J , ( rr kr ) K0 (bkr )} + lim{r; kr J0(rr kr)K, (bkr)} } +U(co,kr^») (D-8) Asymptotic expansions for large arguments of the involved Bessel functions and the modified Bessel functions read as (Abramowitz and Stegun, 1970), J0(r(.kr) ~ [2/7rr; kr ]0-5 {cos(r,.kr) + 0 ( | r - k r | - 1 )} J, (r; kr) ~ [2/jrrf kr ] °-5 {sin(r; kr ) + 0( | r. kr | " 1 ) } K0(bkr) ~ [7r/2bkr]0-5 exp(-bkr){l + 0(|bkr|-1)} K,(bkr) ~ [7r/2bkr]°-5 exp(-bkr){l +0(|bkr|-1)} Therefore, the leading behaviours of the limiting expressions in (D-8) are given by, lim {bkr [2/irv, kr ] °-5 [ 7r/2bkr ] °-5 exp(-bkr)} + lim {rf kr[2/7rr; kr ]°-5[7r/2bkr ]0-5 exp(-bkr)} = 0 (D-9) that is to say, C = 0 91 APPENDIX E Inner Product Matrices for the Inverse Smallest and Flattest Model Constructions The computation of the inner product matrices for the inverse smallest and flattest problems parallels that made for the forward problems. Indeed, by replacing r; by kr; as the parameter and kr by r as the independent variable appearing in the formulation of the forward problems, one readily obtains the inner product matrices for the inverse problems. In particular, the entries of the T matrix for the smallest inverse transform are given by, rVj(ej) = [Wo; 0j]/[(k2+b2+k2)2 - 4kr2k2.]0-5 (E-1) After the substitutions kr(=cjsin (7,- )/V and b=ojpc =a>sin (c ) /V, where c is a constant such that sin(c)>0, and V is the surface P-wave velocity, (E-1) yields r;j(cj) = V2/CJ2 {[ I/a; cjj]/[ (sin2 (7,. )+sin2(c)+sin2 (7j )2 -4sin2(7; )sin2(7i )]0-5}. We write rr(u) = V2/GJ2 F • (E-2) where, rr= [ 1/of a; ]/[ (sin2 (7, )+sin2 (c)+sin2 (7. )2 -4sin2 (7- )-sin2 (75 ) ]0-5 92 Since the frequency-dependent term in (E-2) is factored out of the inner product matrix, the inverse r~ 1 (a>) = CJ2/V2 T" 1 is calculated with minimal computational effort. Similarly, the entries of the inner product matrix for the inverse flattest problem are given by, r,.(«) = [ l/o; o: ]4b/[ (k2+b2+k2)2 - 4kr2k2]1-5 (E-3) Upon the substitutions b=cjsin(c)/V and kr- =u>sin(7r )/V we obtain from (E-3), rv(cj) = V5/u5 r;j where T;j = [ i/o; Oj ]4sin(c)/[ (sin2 (7,- )+sin2 (c)+sin2 (y± ) -4 s i n2 (7 r ) s i n2 (7^ ) ] '-5.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Plane-wave decomposition and reconstruction of spherical-wave...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Plane-wave decomposition and reconstruction of spherical-wave seismograms as a linear inverse problem Cabrera, Jose Julian 1983-05-10
pdf
Page Metadata
Item Metadata
Title | Plane-wave decomposition and reconstruction of spherical-wave seismograms as a linear inverse problem |
Creator |
Cabrera, Jose Julian |
Publisher | University of British Columbia |
Date Issued | 1983 |
Description | The plane-wave decomposition of the vertical displacement component of the spherical-wave field corresponding to a compressional point source is solved as a set of inverse problems. The solution utilizes the power and stability of the Backus & Gilbert (smallest and flattest) model-construction techniques, and achieves computational efficiency through the use of analytical solutions of the integrals which are involved. The theory and algorithms developed in this work allow stable and efficient reconstruction of the spherical-wave field from a relatively sparse set of their plane-wave components. However, the algorithms do not formally conserve the correct amplitudes of the seismograms. Comparison of the algorithms with direct integration of the Hankel transform shows very little or no advantage for the transformation from the time-distance (t-x) domain to the delay time - angle of emergence (τ-γ) domain if the seismograms are equi-sampled spatially. However, for cases where the observed seismograms are not equally spaced and for the transformation τ-γto t-x, the proposed schemes are superior to the direct integration of the Hankel transform. Applicability of the algorithms to reflection seismology is demonstrated via the solution to the problem of trace interpolation and that of separation of converted S modes from other modes presented in common-source gathers. In both cases the application of the algorithms to a set of synthetic reflection seismograms yields satisfactory results. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-05-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | 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. |
DOI | 10.14288/1.0052984 |
URI | http://hdl.handle.net/2429/24556 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1984_A6_7 C33.pdf [ 4.21MB ]
- Metadata
- JSON: 831-1.0052984.json
- JSON-LD: 831-1.0052984-ld.json
- RDF/XML (Pretty): 831-1.0052984-rdf.xml
- RDF/JSON: 831-1.0052984-rdf.json
- Turtle: 831-1.0052984-turtle.txt
- N-Triples: 831-1.0052984-rdf-ntriples.txt
- Original Record: 831-1.0052984-source.json
- Full Text
- 831-1.0052984-fulltext.txt
- Citation
- 831-1.0052984.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.831.1-0052984/manifest