Claerbout - Fourier Transforms And Waves

FOURIER TRANSFORMS AND WAVES: in four lectures Jon F. Clærbout Cecil and Ida Green Professor of Geophysics Stanford University c January 18, 1999

Contents 1 Convolution and Spectra

1

1.1

SAMPLED DATA AND Z-TRANSFORMS . . . . . . . . . . . . . . . . .

1

1.2

FOURIER SUMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3

FOURIER AND Z-TRANSFORM . . . . . . . . . . . . . . . . . . . . . .

8

1.4

CORRELATION AND SPECTRA . . . . . . . . . . . . . . . . . . . . . . 11

2 Discrete Fourier transform

17

2.1

FT AS AN INVERTIBLE MATRIX . . . . . . . . . . . . . . . . . . . . . 17

2.2

INVERTIBLE SLOW FT PROGRAM . . . . . . . . . . . . . . . . . . . . 20

2.3

SYMMETRIES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.4

TWO-DIMENSIONAL FT . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3 Downward continuation of waves

29

3.1

DIPPING WAVES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.2

DOWNWARD CONTINUATION . . . . . . . . . . . . . . . . . . . . . . 32

3.3

A matlab program for downward continuation . . . . . . . . . . . . . . . . 36

3.4

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.6

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.7

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.8

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.9

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

CONTENTS

3.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.14 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.17 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Index

39

Why Geophysics uses Fourier Analysis When earth material properties are constant in any of the cartesian variables it is useful to Fourier transform (FT) that variable.

then

In seismology, the earth does not change with time (the ocean does!) so for the earth, we can generally gain by Fourier transforming the time axis thereby converting time-dependent differential equations (hard) to algebraic equations (easier) in frequency (temporal frequency). In seismology, the earth generally changes rather strongly with depth, so we cannot usefully Fourier transform the depth axis and we are stuck with differential equations in . On the other hand, we can model a layered earth where each layer has material properties that are constant in . Then we get analytic solutions in layers and we need to patch them together. Thirty years ago, computers were so weak that we always Fourier transformed the and coordinates. That meant that their analyses were limited to earth models in which velocity was horizontally layered. Today we still often Fourier transform but not , so we reduce the partial differential equations of physics to ordinary differential equations (ODEs). A big advantage of knowing FT theory is that it enables us to visualize physical behavior without us needing to use a computer. The Fourier transform variables are called frequencies. For each axis have a corresponding frequency . The ’s are spatial frequencies, temporal frequency.

we is the

The frequency is inverse to the wavelength. Question: A seismic wave from the fast earth goes into the slow ocean. The temporal frequency stays the same. What happens to the spatial frequency (inverse spatial wavelength)? In a layered earth, the horizonal spatial frequency is a constant function of depth. We will find this to be Snell’s law. In a spherical coordinate system or a cylindrical coordinate system, Fourier transforms are useless but they are closely related to “spherical harmonic functions” and Bessel transformations which play a role similar to FT. Our goal for these four lectures is to develop Fourier transform insights and use them to take observations made on the earth’s surface and “downward continue” them, to extrapi

CONTENTS

ii

olate them into the earth. This is a central tool in earth imaging.

0.0.1 Impulse response and ODEs When Fourier transforms are applicable, it means the “earth response” now is the same as the earth response later. Switching our point of view from time to space, the applicability of Fourier transformation means that the “impulse response” here is the same as the impulse response there. An impulse is a column vector full of zeros with somewhere a one, say (where the prime means transpose the row into a column.) An impulse response is a column from the matrix

(0.1)

The impulse response is the that comes out when the input is an impulse. In a typical application, the matrix would be about and not the simple example that I show you above. Notice that each column in the matrix contains the same waveform . This waveform is called the “impulse response”. The collection of impulse responses in Equation (0.1) defines the the convolution operation. Not only do the columns of the matrix contain the same impulse response, but each row likewise contains the same thing, and that thing is the backwards impulse response . Suppose were numerically equal to . Then equation (0.1) would be like the differential equation . Equation (0.1) would be a finitedifference representation of a differential equation. Two important ideas are equivalent; either they are both true or they are both false: 1. The columns of the matrix all hold the same impulse response. 2. The differential equation has constant coefficients. The story gets more complicated when we look at the boundaries, the top and bottom few equations. We’ll postpone that.

0.0.2 Z transforms There is another way to think about equation (0.1) which is even more basic. It does not involve physics, differential equations, or impulse responses; it merely involves polynomials.

CONTENTS

iii

(That takes me back to middle school.) Let us define three polynomials. (0.2) (0.3) (0.4) Are you able to multiply ? If you are, then you can examine the coefficient of . You will discover that it is exactly the fifth row of equation (0.1)! Actually it is the sixth row because we started from zero. For each power of in we get one of the rows in equation (0.1). Convolution is defined to be the operation on polynomial coefficients when we multiply polynomials.

0.0.3 Frequency The numerical value of doesn’t matter. It could have any numerical value. We haven’t needed to have any particular value. It happens that real values of lead to what are called Laplace transforms and complex values of lead to Fourier transforms. Let us test some numerical values of . Taking we notice the earliest coefficient in each of the polynomials is strongly emphasized in creating the numerical value of the polynomial, i.e., . Likewise taking , the latest value is strongly emphasized. This undesirable weighting of early or late is avoided if we use the Fourier approach and use numerical values of that fulfill the condition . Other than that forces us to use complex values of , but there are plenty of those. Recall the complex plane where the real axis is horizontal and the imaginary axis is vertical. For Fourier transforms, we are interested in complex numerical values of which have unit magnitude, namely, . Examples are , or . The numerical value

gives what is called the zero frequency. Evaluating , finds the zero-frequency component of . The value gives what is called the “Nyquist frequency”. . The Nyquist frequency is the highest frequency that we can represent with sampled time functions. If our signal were then all the terms in would add together with the same polarity so that signal has a strong frequency component at the Nyquist frequency. How about frequencies inbetween zero and Nyquist? These require us to use complex numbers. Consider , where . The signal could be segregated into its real and imaginary parts. The real part is . Its wavelength is twice as long as that of the Nyquist frequency so its frequency is exactly half. The values for used by Fourier transform are . Now we will steal parts of Jon Claerbout’s books, “Earth Soundings Analysis, Process-

CONTENTS

iv

ing versus Inversion” and “Basic Earth Imaging” which are freely available on the WWW1 . To speed you along though, I trim down those chapters to their most important parts.

1

http://sepwww.stanford.edu/sep/prof/

Chapter 1 Convolution and Spectra Time and space are ordinarily thought of as continuous, but for the purposes of computer analysis we must discretize these axes. This is also called “sampling” or “digitizing.” You might worry that discretization is a practical evil that muddies all later theoretical analysis. Actually, physical concepts have representations that are exact in the world of discrete mathematics.

1.1 SAMPLED DATA AND Z-TRANSFORMS Consider the idealized and simplified signal in Figure 1.1. To analyze such an observed

Figure 1.1: A continuous signal sampled at uniform time intervals. cs-triv1 [ER] signal in a computer, it is necessary to approximate it in some way by a list of numbers. The usual way to do this is to evaluate or observe at a uniform spacing of points in time, call this discretized signal . For Figure 1.1, such a discrete approximation to the continuous function could be denoted by the vector (1.1) Naturally, if time points were closer together, the approximation would be more accurate. What we have done, then, is represent a signal by an abstract -dimensional vector. Another way to represent a signal is as a polynomial, where the coefficients of the polynomial represent the value of at successive times. For example, (1.2) 1

2

CHAPTER 1. CONVOLUTION AND SPECTRA

This polynomial is called a “ -transform.” What is the meaning of here? should not take on some numerical value; it is instead the unit-delay operator. For example, the coefficients of are plotted in Figure 1.2. Figure 1.2 shows

Figure 1.2: The coefficients of are the shifted version of the coefficients of . cs-triv2 [ER] the same waveform as Figure 1.1, but now the waveform has been delayed. So the signal is delayed time units by multiplying by . The delay operator is important in analyzing waves simply because waves take a certain amount of time to move from place to place. Another value of the delay operator is that it may be used to build up more complicated signals from simpler ones. Suppose represents the acoustic pressure function or the seismogram observed after a distant explosion. Then is called the “impulse response.” If another explosion occurred at time units after the first, we would expect the pressure function depicted in Figure 1.3. In terms of -transforms, this pressure function would be expressed as .

Figure 1.3: Response to two explosions. cs-triv3 [ER]

1.1.1 Linear superposition If the first explosion were followed by an implosion of half-strength, we would have . If pulses overlapped one another in time (as would be the case if had degree greater than 10), the waveforms would simply add together in the region of overlap. The supposition that they would just add together without any interaction is called the “linearity” property. In seismology we find that—although the earth is a heterogeneous conglomeration of rocks of different shapes and types—when seismic waves travel through the earth, they do not interfere with one another. They satisfy linear superposition. The plague of nonlinearity arises from large amplitude disturbances. Nonlinearity is a dominating feature in hydrodynamics, where flow velocities are a noticeable fraction of the wave velocity. Nonlinearity is absent from reflection seismology except within a few meters from the source. Nonlinearity does not arise from geometrical complications in the propagation path. An example of two plane waves superposing is shown in Figure 1.4.

1.1. SAMPLED DATA AND Z-TRANSFORMS

3

Figure 1.4: Crossing plane waves superposing viewed on the left as “wiggle traces” and on the right as “raster.” cs-super [ER]

1.1.2 Convolution with Z-transform Now suppose there was an explosion at , a half-strength implosion at , and another, quarter-strength explosion at . This sequence of events determines a “source” time series, . The -transform of the source is . The observed for this sequence of explosions and implosions through the seismometer has a -transform , given by

(1.3) The last equation shows polynomial multiplication as the underlying basis of time-invariant linear-system theory, namely that the output can be expressed as the input times the impulse-response filter . When signal values are insignificant except in a “small” region on the time axis, the signals are called “wavelets.”

1.1.3 Convolution equation and program What do we actually do in a computer when we multiply two -transforms together? The filter would be represented in a computer by the storage in memory of the coefficients . Likewise, for , the numbers would be stored. The polynomial multiplication program should take these inputs and produce the sequence . Let us see how the computation proceeds in a general case, say (1.4) (1.5) Identifying coefficients of successive powers of , we get

4

CHAPTER 1. CONVOLUTION AND SPECTRA

(1.6)

In matrix form this looks like

(1.7)

The following equation, called the “convolution equation,” carries the spirit of the group shown in (1.6): (1.8) To be correct in detail when we associate equation (1.8) with the group (1.6), we should also assert that either the input vanishes before or must be adjusted so that the sum does not extend before . These end conditions are expressed more conveniently by defining in equation (1.8) and eliminating getting (1.9) A convolution program based on equation (1.9) including end effects on both ends, is convolve(). # convolution: Y(Z) = X(Z) * B(Z) # subroutine convolve( nb, bb, nx, xx, yy ) integer nb # number of coefficients in integer nx # number of coefficients in # number of coefficients in real bb(nb) # filter coefficients real xx(nx) # input trace real yy(1) # output trace integer ib, ix, iy, ny ny = nx + nb -1 call null( yy, ny) do ib= 1, nb do ix= 1, nx yy( ix+ib-1) = yy( ix+ib-1) return; end

filter input output will be nx+nb-1

+ xx(ix) * bb(ib)

1.2. FOURIER SUMS

5

This program is written in a language called Ratfor, a “rational” dialect of Fortran. It is similar to the Matlab language. You are not responsible for anything in this program, but, if you are interested, more details in the last chapter of PVI1 , the book that I condensed this from.

1.1.4 Negative time Notice that and need not strictly be polynomials; they may contain both positive and negative powers of , such as (1.10) (1.11) The negative powers of in and show that the data is defined before . The effect of using negative powers of in the filter is different. Inspection of (1.8) shows that the output that occurs at time is a linear combination of current and previous inputs; that is, . If the filter had included a term like , then the output at time would be a linear combination of current and previous inputs and , an input that really has not arrived at time . Such a filter is called a “nonrealizable” filter, because it could not operate in the real world where nothing can respond now to an excitation that has not yet occurred. However, nonrealizable filters are occasionally useful in computer simulations where all the data is prerecorded.

1.2 FOURIER SUMS The world is filled with sines and cosines. The coordinates of a point on a spinning wheel are , where is the angular frequency of revolution and is the phase angle. The purest tones and the purest colors are sinusoidal. The movement of a pendulum is nearly sinusoidal, the approximation going to perfection in the limit of small amplitude motions. The sum of all the tones in any signal is its “spectrum.” Small amplitude signals are widespread in nature, from the vibrations of atoms to the sound vibrations we create and observe in the earth. Sound typically compresses air by a volume fraction of to . In water or solid, the compression is typically to . A mathematical reason why sinusoids are so common in nature is that laws of nature are typically expressible as partial differential equations. Whenever the coefficients of the differentials (which are functions of material properties) are constant in time and space, the equations have exponential and sinusoidal solutions that correspond to waves propagating in all directions. 1

http://sepwww.stanford.edu/sep/prof/

6

CHAPTER 1. CONVOLUTION AND SPECTRA

1.2.1 Superposition of sinusoids Fourier analysis is built from the complex exponential (1.12) A Fourier component of a time signal is a complex number, a sum of real and imaginary parts, say (1.13) which is attached to some frequency. Let be an integer and be a set of frequencies. A signal can be manufactured by adding a collection of complex exponential signals, each complex exponential being scaled by a complex coefficient , namely, (1.14) This manufactures a complex-valued signal. How do we arrange for to be real? We can throw away the imaginary part, which is like adding to its complex conjugate , and then dividing by two: (1.15) In other words, for each positive with amplitude , we add a negative with amplitude (likewise, for every negative ...). The are called the “frequency function,” or the “Fourier transform.” Loosely, the are called the “spectrum,” though in formal mathematics, the word “spectrum” is reserved for the product . The words “amplitude spectrum” universally mean . In practice, the collection of frequencies is almost always evenly spaced. Let integer so that

be an (1.16)

Representing a signal by a sum of sinusoids is technically known as “inverse Fourier transformation.” An example of this is shown in Figure 1.5.

1.2.2 Sampled time and Nyquist frequency In the world of computers, time is generally mapped into integers too, say . This is called “discretizing” or “sampling.” The highest possible frequency expressible on a mesh is , which is the same as . Setting , we see that the maximum frequency is (1.17)

1.2. FOURIER SUMS

7

Figure 1.5: Superposition of two sinusoids. cs-cosines [NR]

Time is commonly given in either seconds or sample units, which are the same when . In applications, frequency is usually expressed in cycles per second, which is the same as Hertz, abbreviated Hz. In computer work, frequency is usually specified in cycles per sample. In theoretical work, frequency is usually expressed in radians where the relation between radians and cycles is . We use radians because, otherwise, equations are filled with ’s. When time is given in sample units, the maximum frequency has a name: it is the “Nyquist frequency,” which is radians or cycle per sample.

1.2.3 Fourier sum In the previous section we superposed uniformly spaced frequencies. Now we will superpose delayed impulses. The frequency function of a delayed impulse at time delay is . Adding some pulses yields the “Fourier sum”: (1.18)

The Fourier sum transforms the signal to the frequency function . Time will often be denoted by , even though its units are sample units instead of physical units. Thus we often see in equations like (1.18) instead of , resulting in an implied .

CHAPTER 1. CONVOLUTION AND SPECTRA

8

1.3 FOURIER AND Z-TRANSFORM The frequency function of a pulse at time is occurs so often in applied work that it has a name:

. The factor

(1.19) With this , the pulse at time is compactly represented as . The variable makes Fourier transforms look like polynomials, the subject of a literature called “ -transforms.” The -transform is a variant form of the Fourier transform that is particularly useful for time-discretized (sampled) functions. From the definition (1.19), we have lencies, equation (1.18) becomes

,

, etc. Using these equiva-

(1.20)

1.3.1 Unit circle In this chapter, is a real variable, so variable. It has unit magnitude because ranges on the unit circle .

. As

is a complex ranges on the real axis,

1.3.2 Differentiator A particularly interesting factor is , because the filter is like a time derivative. The time-derivative filter destroys zero frequency in the input signal. The zero frequency is with a -transform . To see that the filter destroys zero frequency, notice that . More formally, consider output made from the filter and any input . Since vanishes at , then likewise must vanish at . Vanishing at is vanishing at frequency because from (1.19). Now we can recognize that multiplication of two functions of or of is the equivalent of convolving the associated time functions. Multiplication in the frequency domain is convolution in the time domain. A popular mathematical abbreviation for the convolution operator is an asterisk: equation (1.8), for example, could be denoted by . I do not disagree with asterisk notation, but I prefer the equivalent expression , which simultaneously exhibits the time domain and the frequency domain. The filter

is often called a “differentiator.” It is displayed in Figure 1.6.

1.3. FOURIER AND Z-TRANSFORM

9

Figure 1.6: A discrete representation of the first-derivative operator. The filter is plotted on the left, and on the right is an amplitude response, i.e., versus . cs-ddt [NR]

1.3.3 Gaussian examples The filter is a running average of two adjacent time points. Applying this filter times yields the filter . The coefficients of the filter are generally known as Pascal’s triangle. For large the coefficients tend to a mathematical limit known as a Gaussian function, , where and are constants that we will not determine here. We will not prove it here, but this Gaussian-shaped signal has a Fourier transform that also has a Gaussian shape, . The Gaussian shape is often called a “bell shape.” Figure 1.7 shows an example for . Note that, except for the rounded ends, the bell shape seems a good fit to a triangle function. Curiously, the filter

Figure 1.7: A Gaussian approximated by many powers of

. cs-gauss [NR]

also tends to the same Gaussian but with a different . A mathematical theorem says that almost any polynomial raised to the -th power yields a Gaussian. In seismology we generally fail to observe the zero frequency. Thus the idealized seismic pulse cannot be a Gaussian. An analytic waveform of longstanding popularity in seismology is the second derivative of a Gaussian, also known as a “Ricker wavelet.” Starting from the Gaussian and multiplying be produces this old, favorite wavelet, shown in Figure 1.8.

CHAPTER 1. CONVOLUTION AND SPECTRA

10

Figure 1.8: Ricker wavelet. cs-ricker [NR]

1.3.4 Inverse Z-transform Fourier analysis is widely used in mathematics, physics, and engineering as a Fourier integral transformation pair: (1.21) (1.22) These integrals correspond to the sums we are working with here except for some minor details. Books in electrical engineering redefine as . That is like switching to . Instead, we have chosen the sign convention of physics, which is better for wavepropagation studies (as explained in IEI). The infinite limits on the integrals result from expressing the Nyquist frequency in radians/second as . Thus, as tends to zero, the Fourier sum tends to the integral. When we reach equation (??) we will see that if a scaling divisor of is introduced into either (1.21) or (1.22), then will equal . The -transform is always easy to make, but the Fourier integral could be difficult to perform, which is paradoxical, because the transforms are really the same. To make a -transform, we merely attach powers of to successive data points. When we have , we can refer to it either as a time function or a frequency function. If we graph the polynomial coefficients, then it is a time function. It is a frequency function if we evaluate and graph the polynomial for various frequencies .

EXERCISES: 1 Let . Graph the coefficients of the powers of . Graph the coefficients of . 2 As moves from zero to positive frequencies, where is around the unit circle, clockwise or counterclockwise?

as a function of

and which way does it rotate

1.4. CORRELATION AND SPECTRA

11

3 Identify locations on the unit circle of the following frequencies: (1) the zero frequency, (2) the Nyquist frequency, (3) negative frequencies, and (4) a frequency sampled at 10 points per wavelength. 4 Sketch the amplitude spectrum of Figure 1.8 from 0 to

.

1.4 CORRELATION AND SPECTRA The spectrum of a signal is a positive function of frequency that says how much of each tone is present. The Fourier transform of a spectrum yields an interesting function called an “autocorrelation,” which measures the similarity of a signal to itself shifted.

1.4.1 Spectra in terms of Z-transforms Let us look at spectra in terms of -transforms. Let a spectrum be denoted

, where (1.23)

Expressing this in terms of a three-point -transform, we have (1.24) (1.25) (1.26) It is interesting to multiply out the polynomial coefficients of :

with

in order to examine the

(1.27) The coefficient

of

is given by (1.28)

Equation (1.28) is the autocorrelation formula. The autocorrelation value at lag is . It is a measure of the similarity of with itself shifted units in time. In the most frequently occurring case, is real; then, by inspection of (1.28), we see that the autocorrelation coefficients are real, and .

CHAPTER 1. CONVOLUTION AND SPECTRA

12

Specializing to a real time series gives (1.29) (1.30) (1.31) (1.32) cosine transform of

(1.33)

This proves a classic theorem that for real-valued signals can be simply stated as follows: For any real signal, the cosine transform of the autocorrelation equals the magnitude squared of the Fourier transform.

1.4.2 Two ways to compute a spectrum There are two computationally distinct methods by which we can compute a spectrum: (1) compute all the coefficients from (1.28) and then form the cosine sum (1.32) for each ; and alternately, (2) evaluate for some value of Z on the unit circle, and multiply the resulting number by its complex conjugate. Repeat for many values of on the unit circle. When there are more than about twenty lags, method (2) is cheaper, because the fast Fourier transform discussed in chapter 2 can be used.

1.4.3 Common signals Figure 1.9 shows some common signals and their autocorrelations. Figure 1.10 shows the cosine transforms of the autocorrelations. Cosine transform takes us from time to frequency and it also takes us from frequency to time. Thus, transform pairs in Figure 1.10 are sometimes more comprehensible if you interchange time and frequency. The various signals are given names in the figures, and a description of each follows: cos The theoretical spectrum of a sinusoid is an impulse, but the sinusoid was truncated (multiplied by a rectangle function). The autocorrelation is a sinusoid under a triangle, and its spectrum is a broadened impulse (which can be shown to be a narrow sinc-squared function). sinc The sinc function is . Its autocorrelation is another sinc function, and its spectrum is a rectangle function. Here the rectangle is corrupted slightly by “Gibbs sidelobes,” which result from the time truncation of the original sinc. wide box A wide rectangle function has a wide triangle function for an autocorrelation and a narrow sinc-squared spectrum.

1.4. CORRELATION AND SPECTRA

13

Figure 1.9: Common signals and one side of their autocorrelations. cs-autocor [ER]

Figure 1.10: Autocorrelations and their cosine transforms, i.e., the (energy) spectra of the common signals. cs-spectra [ER]

14

CHAPTER 1. CONVOLUTION AND SPECTRA

narrow box A narrow rectangle has a wide sinc-squared spectrum. twin Two pulses. 2 boxes Two separated narrow boxes have the spectrum of one of them, but this spectrum is modulated (multiplied) by a sinusoidal function of frequency, where the modulation frequency measures the time separation of the narrow boxes. (An oscillation seen in the frequency domain is sometimes called a “quefrency.”) comb Fine-toothed-comb functions are like rectangle functions with a lower Nyquist frequency. Coarse-toothed-comb functions have a spectrum which is a fine-toothed comb. exponential The autocorrelation of a transient exponential function is a double-sided exponential function. The spectrum (energy) is a Cauchy function, . The curious thing about the Cauchy function is that the amplitude spectrum diminishes inversely with frequency to the first power; hence, over an infinite frequency axis, the function has infinite integral. The sharp edge at the onset of the transient exponential has much high-frequency energy. Gauss The autocorrelation of a Gaussian function is another Gaussian, and the spectrum is also a Gaussian. random Random numbers have an autocorrelation that is an impulse surrounded by some short grass. The spectrum is positive random numbers. smoothed random Smoothed random numbers are much the same as random numbers, but their spectral bandwidth is limited.

1.4.4 Spectra of complex-valued signals The spectrum of a signal is the magnitude squared of the Fourier transform of the function. Consider the real signal that is a delayed impulse. Its -transform is simply ; so the real part is , and the imaginary part is . The real part is thus an even function of frequency and the imaginary part an odd function of frequency. This is also true of and any sum of powers (weighted by real numbers), and thus it is true of any time function. For any real signal, therefore, the Fourier transform has an even real part RE and an imaginary odd part IO. Taking the squared magnitude gives (RE+ IO)(RE IO)= (RE + (IO . The square of an even function is obviously even, and the square of an odd function is also even. Thus, because the spectrum of a real-time function is even, its values at plus frequencies are the same as its values at minus frequencies. In other words, no special meaning should be attached to negative frequencies. This is not so of complex-valued signals. Although most signals which arise in applications are real signals, a discussion of correlation and spectra is not mathematically complete without considering complex-valued signals. Furthermore, complex-valued signals arise in many different contexts. In seismology, they arise in imaging studies when the space axis is Fourier transformed, i.e., when

1.4. CORRELATION AND SPECTRA

15

a two-dimensional function is Fourier transformed over space to . More generally, complex-valued signals arise where rotation occurs. For example, consider two vector-component wind-speed indicators: one pointing north, recording , and the other pointing west, recording . Now, if we make a complex-valued time series , the magnitude and phase angle of the complex numbers have an obvious physical interpretation: corresponds to rotation in one direction (counterclockwise), and to rotation in the other direction. To see why, suppose and . Then . The Fourier transform is (1.34) The integrand oscillates and averages out to zero, except for the frequency frequency function is a pulse at :

. So the (1.35)

Conversely, if were , then the frequency function would be a pulse at meaning that the wind velocity vector is rotating the other way.

,

1.4.5 Time-domain conjugate A complex-valued signal such as can be imagined as a corkscrew, where the real and imaginary parts are plotted on the - and -axes, and time runs down the axis of the screw. The complex conjugate of this signal reverses the -axis and gives the screw an opposite handedness. In -transform notation, the time-domain conjugate is written (1.36) Now consider the complex conjugate of a frequency function. In -transform notation this is written (1.37) To see that it makes a difference in which domain we take a conjugate, contrast the two equations (1.36) and (1.37). The function is a spectrum, whereas the function is called an “envelope function.” You might be tempted to think that it is not.

, but that is true only if

is real, and often

1.4.6 Spectral transfer function Filters are often used to change the spectra of given data. With input , filters and output , we have and the Fourier conjugate . Multiplying these two relations together, we get

,

(1.38)

16

CHAPTER 1. CONVOLUTION AND SPECTRA

which says that the spectrum of the input times the spectrum of the filter equals the spectrum of the output. Filters are often characterized by the shape of their spectra; this shape is the same as the spectral ratio of the output over the input: (1.39)

EXERCISES: 1 Suppose a wavelet is made up of complex numbers. Is the autocorrelation relation true? Is real or complex? Is real or complex?

Chapter 2 Discrete Fourier transform Happily, Fourier sums are exactly invertible: given the output, the input can be quickly found. Because signals can be transformed to the frequency domain, manipulated there, and then returned to the time domain, convolution and correlation can be done faster. Time derivatives can also be computed with more accuracy in the frequency domain than in the time domain. Signals can be shifted a fraction of the time sample, and they can be shifted back again exactly. In this chapter we will see how many operations we associate with the time domain can often be done better in the frequency domain. We will also examine some two-dimensional Fourier transforms.

2.1 FT AS AN INVERTIBLE MATRIX A Fourier sum may be written (2.1) where the complex value is related to the real frequency by . This Fourier sum is a way of building a continuous function of from discrete signal values in the time domain. In this chapter we will study the computational tricks associated with specifying both time and frequency domains by a set of points. Begin with an example of a signal that is nonzero at four successive instants, . The transform is (2.2) The evaluation of this polynomial can be organized as a matrix times a vector, such as

(2.3)

17

18

CHAPTER 2. DISCRETE FOURIER TRANSFORM

Observe that the top row of the matrix evaluates the polynomial at , a point where also . The second row evaluates , where is some base frequency. The third row evaluates the Fourier transform for , and the bottom row for . The matrix could have more than four rows for more frequencies and more columns for more time points. I have made the matrix square in order to show you next how we can find the inverse matrix. The size of the matrix in (2.3) is . If we choose the base frequency and hence correctly, the inverse matrix will be

(2.4)

Multiplying the matrix of (2.4) with that of (2.3), we first see that the diagonals are +1 as desired. To have the off diagonals vanish, we need various sums, such as and , to vanish. Every element ( , for example, or ) is a unit vector in the complex plane. In order for the sums of the unit vectors to vanish, we must ensure that the vectors pull symmetrically away from the origin. A uniform distribution of directions meets this requirement. In other words, should be the -th root of unity, i.e., (2.5) The lowest frequency is zero, corresponding to the top row of (2.3). The next-to-thelowest frequency we find by setting in (2.5) to . So ; and for (2.4) to be inverse to (2.3), the frequencies required are (2.6)

2.1.1 The Nyquist frequency The highest frequency in equation (2.6), , is almost . This frequency is twice as high as the Nyquist frequency . The Nyquist frequency is normally thought of as the “highest possible” frequency, because , for integer , plots as . The double Nyquist frequency function, , for integer , plots as . So this frequency above the highest frequency is really zero frequency! We need to recall that . Thus, all the frequencies near the upper end of the range (2.6) are really small negative frequencies. Negative frequencies on the interval were moved to interval by the matrix form of Fourier summation.

2.1.2 The comb function Consider a constant function of time. In the frequency domain, it is an impulse at zero frequency. The comb function is defined to be zero at alternate time points. Multiply this

2.1. FT AS AN INVERTIBLE MATRIX

19

constant function by the comb function. The resulting signal contains equal amounts of two frequencies; half is zero frequency, and half is Nyquist frequency. We see this in the second row in Figure 2.1, where the Nyquist energy is in the middle of the frequency axis. In the third row, 3 out of 4 points are zeroed by another comb. We now see something like a new Nyquist frequency at half the Nyquist frequency visible on the second row.

Figure 2.1: A zero-frequency function and its cosine transform. Successive rows show increasingly sparse sampling of the zerofrequency function. dft-comb [NR]

2.1.3 Undersampled field data Figure 2.2 shows a recording of an airgun along with its spectrum. The original data

Figure 2.2: Raw data is shown on the top left, of about a half-second duration. Right shows amplitude spectra (magnitude of FT). In successive rows the data is sampled less densely. dft-undersample [ER] is sampled at an interval of 4 milliseconds, which is 250 times per second. Thus, the Nyquist frequency is 125 Hz. Negative frequencies are not shown, since the amplitude spectrum at negative frequency is identical with that at positive frequency. Think of extending the top row of spectra in Figure 2.2 to range from minus 125 Hz to plus 125 Hz. Imagine the even function of frequency centered at zero frequency—we will soon see it. In the second row of the plot, I decimated the data to 8 ms. This drops the Nyquist

20

CHAPTER 2. DISCRETE FOURIER TRANSFORM

frequency to 62.5 Hz. Energy that was at Hz appears at Hz in the second row spectrum. The appearance of what were formerly small negative frequencies near the Nyquist frequency is called “folding” of the spectrum. In the next row the data is sampled at 16 ms intervals, and in the last row at 32 ms intervals. The 8 ms sampling seems OK, whereas the 32 ms sampling looks poor. Study how the spectrum changes from one row to the next. The spectrum suffers no visible harm in the drop from 4 ms to 8 ms. The 8 ms data could be used to construct the original 4 ms data by transforming the 8 ms data to the frequency domain, replacing values at frequencies above 125/2 Hz by zero, and then inverse transforming to the time domain. (Airguns usually have a higher frequency content than we see here. Some high-frequency energy was removed by the recording geometry, and I also removed some when preparing the data.)

2.2 INVERTIBLE SLOW FT PROGRAM Because Fourier sums are exactly invertible, some other things we often require can be done exactly by doing them in the frequency domain. Typically, signals are real valued. But the programs in this chapter are for complexvalued signals. In order to use these programs, copy the real-valued signal into a complex array, where the signal goes into the real part of the complex numbers; the imaginary parts are then automatically set to zero. There is no universally correct choice of scale factor in Fourier transform: choice of scale is a matter of convenience. Equations (2.3) and (2.4) mimic the -transform, so their scaling factors are convenient for the convolution theorem—that a product in the frequency domain is a convolution in the time domain. Obviously, the scaling factors of equations (2.3) and (2.4) will need to be interchanged for the complementary theorem that a convolution in the frequency domain is a product in the time domain.

2.2.1 The slow FT code The slowft() routine exhibits features found in many physics and engineering programs. For example, the time-domain signal (which I call “tt()”), has nt values subscripted, from tt(1) to tt(nt). The first value of this signal tt(1) is located in real physical time at t0. The time interval between values is dt. The value of tt(it) is at time t0+(it-1)*dt. I do not use “if” as a pointer on the frequency axis because if is a keyword in most programming languages. Instead, I count along the frequency axis with a variable named ie. subroutine slowft( adj, add, t0,dt,tt,nt, integer it,ie, adj, add, nt,

f0,df,

ff,nf) nf

2.3. SYMMETRIES

21

complex cexp, cmplx, tt(nt), real pi2, freq, time, scale, t0,dt, f0,df call adjnull( adj, add, tt,2*nt, pi2 = 2. * 3.14159265; scale = 1./sqrt( 1.*nt) df = (1./dt) / nf f0 = - .5/dt do ie = 1, nf { do it = 1, nt {

ff(nf) ff,2*nf)

freq= f0 + df*(ie-1) time= t0 + dt*(it-1)

if( adj == 0 ) ff(ie)= ff(ie) + tt(it) * cexp(cmplx(0., pi2*freq*time)) * scale else tt(it)= tt(it) + ff(ie) * cexp(cmplx(0.,-pi2*freq*time)) * scale }} return; end

The total frequency band is radians per sample unit or Hz. Dividing the total interval by the number of points nf gives . We could choose the frequencies to run from 0 to radians/sample. That would work well for many applications, but it would be a nuisance for applications such as differentiation in the frequency domain, which require multiplication by including the negative frequencies as well as the positive. So it seems more natural to begin at the most negative frequency and step forward to the most positive frequency.

2.3 SYMMETRIES Next we examine odd/even symmetries to see how they are affected in Fourier transform. The even part of a signal is defined as (2.7) The odd part is (2.8) By adding (2.7) and (2.8), we see that a function is the sum of its even and odd parts: (2.9) Consider a simple, real, even signal such as is an even function of , since Consider the real, odd signal is imaginary and odd, since

. Its transform .

. Its transform .

22

CHAPTER 2. DISCRETE FOURIER TRANSFORM

Likewise, the transform of the imaginary even function is the imaginary even function . Finally, the transform of the imaginary odd function is real and odd. Let and refer to real and imaginary, and to even and odd, and lower-case and upper-case letters to time and frequency functions. A summary of the symmetries of Fourier transform is shown in Figure 2.3.

Figure 2.3: Odd functions swap real and imaginary. Even functions do not get mixed up with complex numbers. dft-reRE [NR]

More elaborate signals can be made by adding together the three-point functions we have considered. Since sums of even functions are even, and so on, the diagram in Figure 2.3 applies to all signals. An arbitrary signal is made from these four parts only, i.e., the function has the form . On transformation of , each of the four individual parts transforms according to the table. Most “industry standard” methods of Fourier transform set the zero frequency as the first element in the vector array holding the transformed signal, as implied by equation (2.3). This is a little inconvenient, as we saw a few pages back. The Nyquist frequency is then the first point past the middle of the even-length array, and the negative frequencies lie beyond. Figure 2.4 shows an example of an even function as it is customarily stored.

Figure 2.4: Even functions as customarily stored by “industry standard” FT programs. dft-even [NR]

2.4. TWO-DIMENSIONAL FT

23

2.3.1 Convolution in the frequency domain Let . The coefficients can be found from the coefficients and by convolution in the time domain or by multiplication in the frequency domain. For the latter, we would evaluate both and at uniform locations around the unit circle, i.e., compute Fourier sums and from and . Then we would form for all , and inverse Fourier transform to . The values come out the same as by the time-domain convolution method, roughly that of our calculation precision (typically fourbyte arithmetic or about one part in ). The only way in which you need to be cautious is to use zero padding greater than the combined lengths of and . An example is shown in Figure 2.5. It is the result of a Fourier-domain computation which shows that the convolution of a rectangle function with itself gives a triangle. Notice that the triangle is clean—there are no unexpected end effects. Figure 2.5: Top shows a rectangle transformed to a sinc. Bottom shows the sinc squared, back transformed to a triangle. dft-box2triangle [NR] Because of the fast method of Fourier transform described next, the frequency-domain calculation is quicker when both and have more than roughly 20 coefficients. If either or has less than roughly 20 coefficients, then the time-domain calculation is quicker.

2.4 TWO-DIMENSIONAL FT Let us review some basic facts about two-dimensional Fourier transform. A two-dimensional function is represented in a computer as numerical values in a matrix, whereas a onedimensional Fourier transform in a computer is an operation on a vector. A 2-D Fourier transform can be computed by a sequence of 1-D Fourier transforms. We can first transform each column vector of the matrix and then each row vector of the matrix. Alternately, we can first do the rows and later do the columns. This is diagrammed as follows:

The diagram has the notational problem that we cannot maintain the usual convention of using a lower-case letter for the domain of physical space and an upper-case letter for

24

CHAPTER 2. DISCRETE FOURIER TRANSFORM

the Fourier domain, because that convention cannot include the mixed objects and . Rather than invent some new notation, it seems best to let the reader rely on the context: the arguments of the function must help name the function. An example of two-dimensional Fourier transforms on typical deep-ocean data is shown in Figure 2.6. In the deep ocean, sediments are fine-grained and deposit slowly in

Figure 2.6: A deep-marine dataset from Alaska (U.S. Geological Survey) and the real part of various Fourier transforms of it. Because of the long traveltime through the water, the time axis does not begin at . dft-plane4 [ER] flat, regular, horizontal beds. The lack of permeable rocks such as sandstone severely reduces the potential for petroleum production from the deep ocean. The fine-grained shales overlay irregular, igneous, basement rocks. In the plot of , the lateral continuity of the sediments is shown by the strong spectrum at low . The igneous rocks show a

2.4. TWO-DIMENSIONAL FT

25

spectrum extending to such large that the deep data may be somewhat spatially aliased (sampled too coarsely). The plot of shows that the data contains no low-frequency energy. The dip of the sea floor shows up in -space as the energy crossing the origin at an angle. Altogether, the two-dimensional Fourier transform of a collection of seismograms involves only twice as much computation as the one-dimensional Fourier transform of each seismogram. This is lucky. Let us write some equations to establish that the asserted procedure does indeed do a 2-D Fourier transform. Say first that any function of and may be expressed as a superposition of sinusoidal functions: (2.10) The double integration can be nested to show that the temporal transforms are done first (inside):

The quantity in brackets is a Fourier transform over done for each and every . Alternately, the nesting could be done with the -integral on the inside. That would imply rows first instead of columns (or vice versa). It is the separability of into a product of exponentials that makes the computation this easy and cheap.

2.4.1 Signs in Fourier transforms In Fourier transforming -, -, and -coordinates, we must choose a sign convention for each coordinate. Of the two alternative sign conventions, electrical engineers have chosen one and physicists another. While both have good reasons for their choices, our circumstances more closely resemble those of physicists, so we will use their convention. For the inverse Fourier transform, our choice is (2.11) For the forward Fourier transform, the space variables carry a negative sign, and time carries a positive sign. Let us see the reasons why electrical engineers have made the opposite choice, and why we go with the physicists. Essentially, engineers transform only the time axis, whereas physicists transform both time and space axes. Both are simplifying their lives by their choice of sign convention, but physicists complicate their time axis in order to simplify their many space axes. The engineering choice minimizes the number of minus signs associated with the time axis, because for engineers, is associated with instead of, as is the case for us and for physicists, with . We confirm this with equation (2.11). Physicists and

26

CHAPTER 2. DISCRETE FOURIER TRANSFORM

geophysicists deal with many more independent variables than time. Besides the obvious three space axes are their mutual combinations, such as midpoint and offset. You might ask, why not make all the signs positive in equation (2.11)? The reason is that in that case waves would not move in a positive direction along the space axes. This would be especially unnatural when the space axis was a radius. Atoms, like geophysical sources, always radiate from a point to infinity, not the other way around. Thus, in equation (2.11) the sign of the spatial frequencies must be opposite that of the temporal frequency. The only good reason I know to choose the engineering convention is that we might compute with an array processor built and microcoded by engineers. Conflict of sign convention is not a problem for the programs that transform complex-valued time functions to complex-valued frequency functions, because there the sign convention is under the user’s control. But sign conflict does make a difference when we use any program that converts real-time functions to complex frequency functions. The way to live in both worlds is to imagine that the frequencies produced by such a program do not range from to as the program description says, but from to . Alternately, we could always take the complex conjugate of the transform, which would swap the sign of the -axis.

2.4.2 Examples of 2-D FT An example of a two-dimensional Fourier transform of a pulse is shown in Figure 2.7.

Figure 2.7: A broadened pulse (left) and the real part of its FT (right). [ER]

dft-ft2dofpulse

Notice the location of the pulse. It is closer to the time axis than the frequency axis. This will affect the real part of the FT in a certain way (see exercises). Notice the broadening of the pulse. It was an impulse smoothed over time (vertically) by convolution with (1,1) and over space (horizontally) with (1,4,6,4,1). This will affect the real part of the FT in another way.

2.4. TWO-DIMENSIONAL FT

27

Another example of a two-dimensional Fourier transform is given in Figure 2.8. This example simulates an impulsive air wave originating at a point on the -axis. We see a wave propagating in each direction from the location of the source of the wave. In Fourier space there are also two lines, one for each wave. Notice that there are other lines which do not go through the origin; these lines are called “spatial aliases.” Each actually goes through the origin of another square plane that is not shown, but which we can imagine alongside the one shown. These other planes are periodic replicas of the one shown.

Figure 2.8: A simulated air wave (left) and the amplitude of its FT (right). [ER]

dft-airwave

EXERCISES: 1 Most time functions are real. Their imaginary part is zero. Show that this means that can be determined from . 2 What would change in Figure 2.7 if the pulse were moved (a) earlier on the -axis, and (b) further on the -axis? What would change in Figure 2.7 if instead the time axis were smoothed with (1,4,6,4,1) and the space axis with (1,1)? 3 What would Figure 2.8 look like on an earth with half the earth velocity? 4 Numerically (or theoretically) compute the two-dimensional spectrum of a plane wave [ ], where the plane wave has a randomly fluctuating amplitude: say, rand is a random number between , and the randomly modulated plane wave is [ ]. 5 Explain the horizontal “layering” in Figure 2.6 in the plot of the “layer” separation? What determines the “layer” slope?

. What determines

28

CHAPTER 2. DISCRETE FOURIER TRANSFORM

Chapter 3 Downward continuation of waves 3.1 DIPPING WAVES We begin with equations to describe a dipping plane wave in a medium of constant velocity. Figure 3.1 shows a ray moving down into the earth at an angle from the vertical. Perpendicular to the ray is a wavefront. By elementary geometry the angle between the

x ray

nt

fro

Figure 3.1: Downgoing ray and wavefront. dwnc-front [NR]

z

wavefront and the earth’s surface is also . The ray increases its length at a speed . The speed that is observable on the earth’s surface is the intercept of the wavefront with the earth’s surface. This speed, namely , is faster than . Likewise, the speed of the intercept of the wavefront and the vertical axis is . A mathematical expression for a straight line like that shown to be the wavefront in Figure 3.1 is (3.1) In this expression is the intercept between the wavefront and the vertical axis. To make the intercept move downward, replace it by the appropriate velocity times time: (3.2) 29

CHAPTER 3. DOWNWARD CONTINUATION OF WAVES

30

Solving for time gives (3.3) Equation (3.3) tells the time that the wavefront will pass any particular location . The expression for a shifted waveform of arbitrary shape is . Using (3.3) to define the time shift gives an expression for a wavefield that is some waveform moving on a ray. moving wavefield

(3.4)

3.1.1 Snell waves In reflection seismic surveys the velocity contrast between shallowest and deepest reflectors ordinarily exceeds a factor of two. Thus depth variation of velocity is almost always included in the analysis of field data. Seismological theory needs to consider waves that are just like plane waves except that they bend to accommodate the velocity stratification . Figure 3.2 shows this in an idealized geometry: waves radiated from the horizontal flight of a supersonic airplane. The airplane passes location at time flying hori-

speed at depth z speed at depth z

1

2

Figure 3.2: Fast airplane radiating a sound wave into the earth. From the figure you can deduce that the horizontal speed of the wavefront is the same at depth as it is at depth . This leads (in isotropic media) to Snell’s law. dwnc-airplane [NR] zontally at a constant speed. Imagine an earth of horizontal plane layers. In this model there is nothing to distinguish any point on the -axis from any other point on the -axis.

3.1. DIPPING WAVES

31

But the seismic velocity varies from layer to layer. There may be reflections, head waves, shear waves, converted waves, anisotropy, and multiple reflections. Whatever the picture is, it moves along with the airplane. A picture of the wavefronts near the airplane moves along with the airplane. The top of the picture and the bottom of the picture both move laterally at the same speed even if the earth velocity increases with depth. If the top and bottom didn’t go at the same speed, the picture would become distorted, contradicting the presumed symmetry of translation. This horizontal speed, or rather its inverse , has several names. In practical work it is called the stepout. In theoretical work it is called the ray parameter. It is very important to note that does not change with depth, even though the seismic velocity does change with depth. In a constant-velocity medium, the angle of a wave does not change with depth. In a stratified medium, does not change with depth. Figure 3.3 illustrates the differential geometry of the wave. The diagram shows that

Figure 3.3: Downgoing fronts and rays in stratified medium horizontal translations of one another. dwnc-frontz [NR]

. The wavefronts are

(3.5) (3.6) These two equations define two (inverse) speeds. The first is a horizontal speed, measured along the earth’s surface, called the horizontal phase velocity. The second is a vertical speed, measurable in a borehole, called the vertical phase velocity. Notice that both these speeds exceed the velocity of wave propagation in the medium. Projection of wave fronts onto coordinate axes gives speeds larger than , whereas projection of rays onto coordinate axes gives speeds smaller than . The inverse of the phase velocities is called the stepout or the slowness. Snell’s law relates the angle of a wave in one layer with the angle in another. The constancy of equation (3.5) in depth is really just the statement of Snell’s law. Indeed, we have just derived Snell’s law. All waves in seismology propagate in a velocity-stratified medium. So they cannot be called plane waves. But we need a name for waves that are near to plane waves. A Snell wave will be defined to be the generalization of a plane

32

CHAPTER 3. DOWNWARD CONTINUATION OF WAVES

wave to a stratified medium . A plane wave that happens to enter a medium of depthvariable velocity gets changed into a Snell wave. While a plane wave has an angle of propagation, a Snell wave has instead a Snell parameter . It is noteworthy that Snell’s parameter is directly observable at the surface, whereas neither nor is directly observable. Since is not only observable, but constant in depth, it is customary to use it to eliminate from equations (3.5) and (3.6): (3.7) (3.8)

3.1.2 Evanescent waves Suppose the velocity increases to infinity at infinite depth. Then equation (3.8) tells us that something strange happens when we reach the depth for which equals . That is the depth at which the ray turns horizontal. Below this critical depth the seismic wavefield damps exponentially with increasing depth. Such waves are called evanescent. For a physical example of an evanescent wave, forget the airplane and think about a moving bicycle. For a bicyclist, the slowness is so large that it dominates for all earth materials. The bicyclist does not radiate a wave, but produces a ground deformation that decreases exponentially into the earth. To radiate a wave, a source must move faster than the material velocity.

3.2 DOWNWARD CONTINUATION Given a vertically upcoming plane wave at the earth’s surface, say , and an assumption that the earth’s velocity is vertically stratified, i.e. , we can presume that the upcoming wave down in the earth is simply time-shifted from what we see on the surface. (This assumes no multiple reflections.) Time shifting can be represented as a linear operator in the time domain by representing it as convolution with an impulse function. In the frequency domain, time shifting is simply multiplying by a complex exponential. This is expressed as (3.9) (3.10)

3.2.1 Continuation of a dipping plane wave. Next consider a plane wave dipping at some angle . It is natural to imagine continuing such a wave back along a ray. Instead, we will continue the wave straight down. This

3.2. DOWNWARD CONTINUATION

33

requires the assumption that the plane wave is a perfect one, namely that the same waveform is observed at all . Imagine two sensors in a vertical well bore. They should record the same signal except for a time shift that depends on the angle of the wave. Notice that the arrival time difference between sensors at two different depths is greatest for vertically propagating waves, and the time difference drops to zero for horizontally propagating waves. So the time shift is where is the angle between the wavefront and the earth’s surface (or the angle between the well bore and the ray). Thus an equation to downward continue the wave is (3.11) (3.12) Equation (3.12) is a downward continuation formula for any angle . We can generalize the method to media where the velocity is a function of depth. Evidently we can apply equation (3.12) for each layer of thickness , and allow the velocity vary with . This is a well known approximation that handles the timing correctly but keeps the amplitude constant (since ) when in real life, the amplitude should vary because of reflection and transmission coefficients. Suffice it to say that in practical earth imaging, this approximation is almost universally satisfactory. In a stratified earth, it is customary to eliminate the angle which is depth variable, and change it to the Snell’s parameter which is constant for all depths. Thus the downward continuation equation for any Snell’s parameter is (3.13) It is natural to wonder where in real life we would encounter a Snell wave that we could downward continue with equation (3.13). The answer is that any wave from real life can be regarded as a sum of waves propagating in all angles. Thus a field data set should first be decomposed into Snell waves of all values of , and then equation (3.13) can be used to downward continue each , and finally the components for each could be added. This process akin to Fourier analysis. We now turn to Fourier analysis as a method of downward continuation which is the same idea but the task of decomposing data into Snell waves becomes the task of decomposing data into sinusoids along the -axis.

3.2.2 Downward continuation with Fourier transform One of the main ideas in Fourier analysis is that an impulse function (a delta function) can be constructed by the superposition of sinusoids (or complex exponentials). In the study of time series this construction is used for the impulse response of a filter. In the study of functions of space, it is used to make a physical point source that can manufacture the downgoing waves that initialize the reflection seismic experiment. Likewise observed upcoming waves can be Fourier transformed over and .

CHAPTER 3. DOWNWARD CONTINUATION OF WAVES

34

Recall that a plane wave carrying an arbitrary waveform is given by equation (3.4). Specializing the arbitrary function to be the real part of the function gives moving cosine wave

(3.14)

Using Fourier integrals on time functions we encounter the Fourier kernel . To use Fourier integrals on the space-axis the spatial angular frequency must be defined. Since we will ultimately encounter many space axes (three for shot, three for geophone, also the midpoint and offset), the convention will be to use a subscript on the letter to denote the axis being Fourier transformed. So is the angular spatial frequency on the -axis and is its Fourier kernel. For each axis and Fourier kernel there is the question of the sign before the . The sign convention used here is the one used in most physics books, namely, the one that agrees with equation (3.14). With this convention, a wave moves in the positive direction along the space axes. Thus the Fourier kernel for -space will be taken to be Fourier kernel

(3.15)

Now for the whistles, bells, and trumpets. Equating (3.14) to the real part of (3.15), physical angles and velocity are related to Fourier components. The Fourier kernel has the form of a plane wave. These relations should be memorized! Angles and Fourier Components

(3.16) A point in -space is a plane wave. The one-dimensional Fourier kernel extracts frequencies. The multi-dimensional Fourier kernel extracts (monochromatic) plane waves. Equally important is what comes next. Insert the angle definitions into the familiar relation . This gives a most important relationship: (3.17) The importance of (3.17) is that it enables us to make the distinction between an arbitrary function and a chaotic function that actually is a wavefield. Imagine any function . Fourier transform it to . Look in the -volume for any nonvanishing values of . You will have a wavefield if and only if all nonvanishing have coordinates that satisfy (3.17). Even better, in practice the -dependence at is usually known, but the -dependence is not. Then the -dependence is found by assuming is a wavefield, so the -dependence is inferred from (3.17). Equation (3.17) also achieves fame as the “dispersion relation of the scalar wave equation,” a topic developed more fully in IEI.

3.2. DOWNWARD CONTINUATION

35

3.2.3 Linking Snell waves to Fourier transforms To link Snell waves to Fourier transforms we merge equations (3.5) and (3.6) with equations (3.16) (3.18) (3.19) The basic downward continuation equation for upcoming waves in Fourier space follows from equation (3.13) by eliminating by using equation (3.18). For analysis of real seismic data we introduce a minus sign because equation (3.19) refers to downgoing waves and observed data is made from up-coming waves. (3.20)

In Fourier space we delay signals by multiplying by , analogously, equation (3.20) says we downward continue signals into the earth by multiplying by . Multiplication in the Fourier domain means convolution in time which can be depicted by the engineering diagram in Figure 3.4.

Figure 3.4: Downward continuation of a downgoing wavefield. dwnc-inout [NR]

Downward continuation is a product relationship in both the -domain and the domain. Thus it is a convolution in both time and . What does the filter look like in the time and space domain? It turns out like a cone, that is, it is roughly an impulse function of . More precisely, it is the Huygens secondary wave source that was exemplified by ocean waves entering a gap through a storm barrier. Adding up the response of multiple gaps in the barrier would be convolution over . A nuisance of using Fourier transforms in migration and modeling is that spaces become periodic. This is demonstrated in Figure 3.5. Anywhere an event exits the frame at a side, top, or bottom boundary, the event immediately emerges on the opposite side. In practice, the unwelcome effect of periodicity is generally ameliorated by padding zeros around the data and the model.

36

CHAPTER 3. DOWNWARD CONTINUATION OF WAVES

Figure 3.5: A reflectivity model on the left and synthetic data using a Fourier method on the right. dwnc-diag [ER]

3.3 A matlab program for downward continuation The matlab program below produces a movie. It shows a wave propating from the surface downward. My plan is to clean it up and make it into a little lab. % Wave Movies for MATLAB (thanks to JianLiang Qian) % : downgoing spherical wave by dispersion relation clear nt=12; nx=48; nz=96; nw=2; nlam=4; dx=2; dz=1; v=1; % allocate memory pp=zeros(nz,nx,nt); q =zeros(nx,1); qx=zeros(nx,1); tqb= zeros(nx,nt); bshift = zeros(1,nt); % column vector % some parameters lambda = nz*dz/nlam; pi2 = 2.*3.141592; dw = v*pi2/lambda; dt = pi2/(nt*dw); x0 = nx*dx/3; z0 = nz*dz/2; kx0 = -0.5*pi2/dx; dkx = pi2/(nx*dx); % Generate the data for the movies for iw = 1:nw % superimpose nw frequencies w= iw*dw; wov = w/v; % frequency / velocity % initial conditions for a collapsing spherical wave qx(1:nx)=exp( 0.-i*wov*sqrt( z0*z0+((1:nx)*dx-x0).*((1:nx)*dx-x0) ) ); bshift = exp( 0.0-i*w*dt*(1:nt) ); for iz = 1:nz % extrapolation in depth q=fftshift(fft(fftshift(qx))); % to Spatial freq. domain for ikx=2:nx kx = kx0+(ikx-1)*dkx; wkx= (kx/wov)*(kx/wov); if( wkx <= 1.0)

3.3. A MATLAB PROGRAM FOR DOWNWARD CONTINUATION

q(ikx)=exp(i*wov*sqrt(1.0-wkx)*dz)*q(ikx); end end qx =fftshift( ifft( fftshift( q))); % back to space domain tqb= qx*bshift; % an update matrix tqb= reshape(tqb,1,nx,nt); % 2d array to 3d array pp(iz,:,:)= pp(iz,:,:) + tqb(1,:,:); % to use 3d array addition end end % % play the movie:(Thanks to Tapan) h=imagesc(real(pp(:,:,1))); colormap gray; set(h,’erasemode’,’none’); for k=1:6*nt index=mod(k,nt)+1; set(h,’cdata’,real(pp(:,:,index))); drawnow; end

37

38

3.4 . 3.5 . 3.6 . 3.7 . 3.8 . 3.9 . 3.10 . 3.11 . 3.12 . 3.13 . 3.14 . 3.15 . 3.16 . 3.17 . 3.18 .

CHAPTER 3. DOWNWARD CONTINUATION OF WAVES

Index -transform, 2, 8 -transform and Fourier transform, 8

and -transform, 8 discrete, 17 two-dimensional, 23, 25, 26 front, 29

airgun, 19 amplitude spectrum, 6 autocorrelation, 11, 12

Gaussian, 9, 14 Gibbs sidelobes, 12

basement rock, 25

Hertz, 7 Hz, 7

Cauchy function, 14 comb, 14, 18 complex-valued signal, 6, 14, 15 conjugate signal in time-domain, 15 convolution, 8 convolve subroutine, 4 corkscrew, 15

impulse response, 2 index, 39 linearity, 2 mesh, 6 negative frequency, 21 nonlinearity, 2 nonrealizable, 5 Nyquist frequency, 7, 10, 18, 19

differentiate, 8 digitizing, 1 dip, 32 double-sided exponential, 14

odd function, 14

envelope, 15 evanescent, 32 even function, 14, 22 exponential, 14 exponential double-sided, 14

Pascal’s triangle, 9 phase velocity, 31 plane wave, 2 polynomial multiplication, 3 precision, 23

filter, 3 filter nonrealizable, 5 fold, 20 Fourier downward continuation, 33 Fourier integral, 10 Fourier sum, 7, 10, 17 Fourier transform, 8 Fourier transform

quefrency, 14 radian, 7, 10 random, 14 ray, 29, 30 ray parameter, 31 realizable, 5 rectangle function, 12 Ricker wavelet, 9 39

40

sampling, 1 scale factor, 20 sign convention, 10, 25 signal complex-valued, 14, 15 sinc, 12 slowft subroutine, 20 slowness, 31 Snell parameter, 32 Snell wave, 31, 33, 35 Snell’s law, 31 spatial alias, 25, 27 spectral ratio, 16 spectrum, 6, 11, 14 spectrum amplitude, 6 stepout, 31 subroutine convolve, convolve, 4 slowft, slow FT, 20 superpose, 2, 7 symmetry, 21 time-domain conjugate, 15 unit-delay operator, 2 wave equation, 34 wavelet, 3 wavelet Ricker, 9 zero frequency, 8, 9 zero pad, 23

INDEX

Jon F. Claerbout (M.I.T., B.S. physics, 1960; M.S. 1963; Ph.D. geophysics, 1967), professor at Stanford University, 1967. Consulted with Chevron (1967-73). Best Presentation Award from the Society of Exploration Geophysicists (SEG) for his paper, Extrapolation of Wave Fields. Honorary member and SEG Fessenden Award “in recognition of his outstanding and original pioneering work in seismic wave analysis.” Founded the Stanford Exploration Project (SEP) in 1973. Elected Fellow of the American Geophysical Union. Published three books: Fundamentals of Geophysical Data Processing, 1976; Imaging the Earth’s Interior 1985, and Processing versus Inversion 1992. Elected to the National Academy of Engineering. Maurice Ewing Medal, SEG’s highest award.

41