1. Q-filtering and absorption-dispersion pairs

This section on Q is based on John Costain and Coruchs’ book “Basic theory of exploration seismology”. The Mathematica programs were written by Costain and has been further developed by Knut Sørsdal who remains responsible for errors of omission or content. We also have some thought from Waters book “Reflection seismology”.

Q stands for quality factor. Formulas that define the dimensionless number Q will be given below, but in simple terms the numerical value of Q is just a direct measure of how much mechanical energy is converted to heat as a seismic disturbance propagates away from the its source. Generally speaking, except for unconsolidated, water-saturated sediments, Q has been commonly thought to be essentially independent of frequency; Q is a constant of the medium in which a disturbance propagates. This kind of absorption of energy is also called intrinsic damping. The energy lost is dissipated and cannot be recovered. There are many reasons why a determination of Q is important. For example, it has been experimentally shown that Q rather than velocity is strongly correlated with permeability (Klimentos and McCann [96]). Thus it should be possible to relate Q anisotropy with anisotropic permeability in partially or completely saturated rocks.

The concepts of absorption and body wave dispersion are introduced in this section. Body wave dispersion means that the different frequency components that make up the seismic wavelet travel at different speeds. Higher frequency components travel faster than lower frequency components. Absorption means that the amplitude of the seismic wavelet decreases along the path of propagation as adjacent grains in a rock matrix rub against each other and mechanical energy is continuously converted to heat. You cannot have absorption without dispersion. The reason for this is simple. Absorption causes the seismic wavelet to spread. It spreads to the right toward increasing time and it spreads to the left toward decreasing time. But it cannot spread into times before the theoretical arrival time of the wavelet. In order to prevent this from happening the different frequencies must travel at different speeds. The lower frequencies travel more slowly thereby keeping the wavelet causal.

The loss in energy has two affects on the seismic wavelet:

1. A gradual decrease in wavelet amplitude along the path of propagation even in a homogeneous, isotropic medium.

2. A change in wavelet shape—that is, the time between two successive peaks or troughs increases. Also, the relative amplitudes of the peaks and troughs change. For higher-Q rocks the time changes are not much, maybe only a fraction of a millisecond, but the change in wavelet shape is impressive even for higher-Q rocks.

The absorption is different for different frequencies so that the wavelet undergoes a change in shape as well as a change in amplitude. A change in shape means that the "time-shifting theorem" is not obeyed. The phenomenon is clearly evident even on a single shot record where it becomes obvious that wavelets recorded at long distances from the sourcepoint are more attenuated and spread out more than those recorded at closer distances. These changes in wavelet shape and amplitude are a result of absorption accompanied by body wave dispersion.

The fact that you cannot have absorption without dispersion means that the phenomenon can be described mathematically by an "absorption-dispersion pair". There has been much discussion in the literature as to what the correct formulation for this "pair" is. It is Ecevitoglu's [66] approach that is described here.

Introduction

The change in amplitude by body wave dispersion is always associated with a change in shape, so that we are not simply scaling the wavelet amplitude—it does not retain the same shape. High-Q rocks (Q = 250, say) are associated with a much smaller energy loss than low-Q rocks (Q = 30, say). "High-quality" rocks distort the wavelet less than "low-quality" rocks. Ecevitoglu [66], using methods of singular value decomposition, found a value of Q = 250 for the rocks of the indurated but relatively unmetamorphosed Paleozoic shelf sequence in southwestern North Carolina. Demirbag [56, 57, 58, 59] found a value of Q = 50 for the relatively unconsolidated sediments of the Atlantic Coastal Plain. (In this context one might say that granites are higher quality rocks than shales.) In this section we are concerned only with the absorption of energy due to the mechanical energy of particle motion that is converted to heat when adjacent grains in a rock matrix rub against each other as the disturbance passes.

The higher frequencies in the wavelet are attenuated more than the lower frequencies but this is not the only reason for the change in wavelet shape. The different frequency components travel at different velocities. Those wavelets that have traveled to more distant receivers are attenuated and changed in shape more than those recorded at near-offset receivers. In other words, the phase shift is not linear with frequency. Clearly then, comparison of the shapes and amplitudes of wavelets recorded at different receiver offsets should contain information about the rock properties traversed by the wavelets. The universally used common-mid-point (CMP) method of acquiring seismic data and then summing traces in a "gather" always sums traces with different source-receiver offsets, and therefore wavelets of different amplitude and shape are summed together even though they have been reflected from the same reflector. This is undesirable but unavoidable; however, a correction for the change in shape and amplitude can be made before summing if the value of Q is known.

Intrinsic damping, like velocity or density, is a measurable rock physical property and if Q can be correlated, like velocity and density, with a particular geologic unit, then using this additional information brings us one step closer to a better geologic interpretation.

There are, of course, other factors that cause changes in wavelet amplitude and shape, for example geometrical spreading of the wavefront, destructive interference with another wavelet, reflection or refraction at an interface just to name a few, and we know from this volume that reflection from a thin bed changes wavelet shape as well as amplitude. Attenuation of energy from the seismic wavelet can also result from layering by a combination of transmission losses and intrabed multiples (O'Doherty and Anstey [127]). Schoenberger and Levin [156] reported that attenuation due to layering accounted for 30 — 50 per cent of the total frequency dependent attenuation estimated from field seismograms.

In our discussion of Q we do not concern ourselves with the physical mechanisms that cause the absorption of energy because the definition of Q does not depend on the details by which the energy is dissipated. We assume here that intrinsic damping can be isolated in the presence of any other amplitude-changing mechanisms, but this is not always an easy task. The extremely clean waveforms associated with many vertical seismic profiles (VSPs), however, clearly show gradual changes in wavelet shape and should allow quite accurate determinations of Q by the mechanism of intrinsic damping. The widespread use of VSPs in itself justifies a closer look at how to measure intrinsic damping. Dispersion affects traveltimes and is known to be a cause of the discrepancy between integrated sonic times and seismic times. This section is devoted to an understanding of body wave dispersion.

Definitions

The most commonly used measure of absorption by intrinsic damping in rocks is the quality factor or specific dissipation constant Q. Other commonly used measures of absorption in rocks due to intrinsic damping are the attenuation coefficient α(f), the internal friction or dissipation factor 1/Q, and the logarithmic decrement δ. The lower the value of Q the more the wavelet is attenuated by intrinsic damping and the more its shape as well as its amplitude changes as the travel distance increases. Again, this is an important effect because in a single gather wavelets from the same reflector recorded at short offsets will have a different shape than those recorded at long offsets. Relatively unconsolidated sediments have low values of Q. Granites have much higher values.

As noted, except for unconsolidated, water-saturated sediments, the quality factor Q is believed to be essentially independent of frequency f and is related to the rate at which the mechanical energy of vibration is converted irreversibly into heat. The above quantities are related as follows:

(1.1)

(1.2)

where V is the non-dispersive constant rock velocity at f = ∞. From (1.2) linear-with-frequency absorption implies a constant Q. The attenuation coefficient α can be expressed in nepers/unit, length (or simply inverse length) or in dB per unit length. The relationship between the two is given by a (dB per unit length) = 8.686 a (nepers per unit length). Also, α(dB/λ) = 8.686 π/Q.

We could also measure attenuation as a percentage of stored energy lost per cycle of a Fourier component , and from our assumption that Q is independent of frequency deduce that this percentage is the same independent of the length of the cycle.

And from here we can sum up and also introduce other measures:

1/Q = specific dissipation constant. It is related to the rate at which the mechanical energy of vibration is converted irreversibly into heat energy and does not depend on the detailed mechanism by which the energy is dissipated.

δ = Logarithmic decrement – the natural logarithm of the ratio of amplitudes of two successive maxima or minima in an exponentially decaying free vibration.

a = damping amplitude coefficient in the expression for a free vibration:

e^-at sin (2ft) (1.3)

Closely related to the damping coefficient is the attenuation coefficient α as a measure of attenuation for plane wave in an infinite medium:

e^-αx sin (2f(t-x/V)) where V = wave velocity (1.4)

∆f/f = Relative bandwidth of a resonance curve between the half-power or 0.707 amplitude points for a solid undergoing forced vibrations is a measure of the sharpness of the resonance curve.

∆E/E = Fraction of strain energy lost per cycle.

Now we can put up a relation between the different parameters on basis of a specific damping capacity b. This damping measure b is related to others by

b = 2π/Q = 2δ = 2 a/f = 2V α/f = 2π ∆f/f = ∆E/E (1.5)

Dispersion

As the higher frequencies are attenuated relative to the lower frequencies, the wavelet changes shape. The change in shape must be such that the wavelet remains causal. In order for this to happen each Fourier frequency component f must travel at a different phase velocity c(f). This dependence of body wave velocity on frequency is called body wave dispersion. From available data it is reasonable to assume that for most rocks over the seismic frequency range of interest q is a linear function of frequency f as defined by 1.2.

The amplitude H₂(f) of a plane seismic wave after having travelled a distance x is

H₂(f) = H₁(f)e-^a(f)x (1.6)

where the amplitude spectra H₁(f) and H₂(f) are measured at points separated by the travel distance x as shown in Figure 1. Taking the log of each side of 1.6 gives

α(f) = (1/x) ln[A₂(f) / A₁(f) ] (1.7)

For most rocks a plot of a(f) versus f appears to be linear with frequency. If a least-squares straight line is fitted to α(f) versus f then α(f) is proportional to f and if x is the (constant) distance between observation points, then

α(f)x=bf (1.8 a)

where b is a constant equal to

b=(πx)/(VQ) (1.8.b)

The relationship between α(f) and f is shown in Figure 2.

Fig.2. Attenuation as a function of frequency

Having defined some well-known relationships we now proceed to derive the expression for the phase velocity c at any frequency f for sampled data. The earth is a filter and it has a frequency response defined by an amplitude term and a phase term. We recall the relationship between the amplitude |H(f)| and phase φ(f) of any function. A complex amplitude spectrum H(f) can be written

H(f) = |H(f)| e ^iφ(f) (1.9)

where H(f) is the complex Fourier spectrum. |H(f)| is the modulus of the spectrum, and φ(f) is the phase angle spectrum. This is simply a definition for any function, causal or not. With h(t) <=> H(f) if h(t) happens to be some wavelet then an estimate of the wavelet amplitude spectrum |H(f)| can be obtained by taking the square root of the power spectrum, which we know to be the Fourier transform of the autocorrelation of the wavelet h(t).

The autocorrelation φ ₁₁(τ) of h(t) is

(1.10)

and the Fourier transform of φ ₁₁(τ) is the power spectrum Ф(f). Thus,

(1.11)

and (1.9) becomes

(1.12)

Now, taking the log of both sides of (1.12) we get

log H(f) = log √Ф(f) + iφ(f) (1.13)

= 1/2 log √Ф(f) + iφ(f) (1.14)

= H_e (f) + i H₀ (f) (1.15)

where H_e(f) = l/21og Ф(f) is even and H₀(f) = φ(f) is odd. We can obtain the phase spectrum of the wavelet by taking the Hilbert transform of 1/2 log Ф(f) = H_e(f). That is,

log H(f) = 1/2 log P(f) + iН 1/2 log P(f) (1.16)

= H_e + i H_o (1.17)

₌ H_e + i H (H_e) (1.18)

where H( ) denotes Hilbert transformation of the argument in parentheses. Having now both the amplitude and phase spectrum we could take the inverse Fourier transform of H(f) to obtain the causal wavelet h(t). The mathematical relationship (Hilbert transformation) between the real and imaginary parts of (8.110) guarantees that the wavelet h(t) will be causal.

Now suppose that h(t) is instead the impulse response of the earth and that h(t)<=> H(f). We don't have access to h(t) but we do know from experiments and observation how the amplitude H (f) of any single Fourier component of frequency f falls off with distance x. An impulse δ(t) leaving the sourcepoint contains all frequencies each with a Fourier amplitude of "1". After a sinusoid of amplitude "1" travels a distance x then from (1.9) its amplitude is given by

(l)e-^α(f)x (1.19)

This is the amplitude part of the complex absorption response H (f) of the earth. Following the notation of (1.12) we can write the general expression for the frequency response H(f) of the earth as

H(f) = (l)e-^α(f)x e ^iφ(f) (1.20)

Now take the log of (1.20) as we did before and get

log H (f) = — α(f)x + iφ(f) (1.21)

But from (1.8.a)

— α(f) x = —bf

So

Log H(f) = — α(f) x + iφ(f) (1.22)

So, - bf is the part of the earth absorption filter that is the analog of 1/2 log Ф(f). When we took the Hilbert transform of log |H(f)| we got a phase spectrum φ(f) with very special properties — the inverse Fourier transform of H(f) will be causal. So (1.22) can be written

Log H(f) = [- α(f)x] + iφ(f) = — bf + iφ(f) (1.23)

This says that the phase spectrum φ(f) of the causal earth absorption filter can be obtained by taking the Hilbert transform of —bf. This phase spectrum is special because of the causal nature of the Hilbert transformation so we assign to it a different symbol B(f):

Because b is a constant the quantity B(f) can be written in several equivalent ways:

B(f) = H(-bf)

= -bH(f)

= bH(-f)

We will arbitrarily adopt the latter form of the expression for B(f) because it looks nice in the final form of the derivation of body wave dispersion (Equation (1.37)). Thus,

Log H(f) = -bf + iB(f) (1.24)

= H_e(f) + iH₀(f) (1.25)

What if the phase term B(f) is ignored and only the amplitude part defined by Equation (1.20) is used? This won't work because the attenuation of frequency components in H(f) results in spreading of the wavelet, as shown in Figure 1.3. The wavelet spreads to the right toward increasing time as well as to the left toward decreasing time. But it cannot spread into times before the wavelet is supposed to arrive based on the known rock velocity and source-receiver distance. In order to prevent this from happening the different frequencies f that make up the spectrum of the wavelet must travel at different speeds. This means that the phase spectrum B(f) cannot be zero.

Figure 3. Absorption without dispersion. The attenuated wavelet without using the B(f) term spreads into negative as well as positive time after inverse Fourier transformation and therefore violates causality.

We can determine H₀ from H_e in (1.25) for continuous data by using the expressions here:

(1.26)

Once we have H₀(f) we can construct H(f) by

H(f) = e- ^{He (f) + Ho (f)}

= e- ^{bf +i H (-bf)}

and then take the inverse Fourier transform to get h(t). That is,

The function h(t) is the response of the Earth to an impulsive input. It is the absorption response of the Earth. The impulse response is causal because of the Hilbert transformation step used to compute the phase spectrum.

Equation (1.26) becomes

(1.27)

Prom (1.23) we substitute H _e(f) = — bf in (1.27), which then becomes

(1.28)

In summary, the Hilbert transform of — bf gives the phase B(f) associated with linear-with-frequency absorption (the causal absorption) of the "Earth filter". Equation (1.28) is a Hilbert transformation. In (1.28), — bf is an even function of f and B(f) is odd. If we take the inverse Fourier transform of H(f) (which we will display) we obtain the absorption response h(t) of the Earth to an impulse.

So far we have just looked at the part B(f) of the total phase spectrum θ(f) associated with causal absorption. There is also a part 2πfτ associated with non-dispersive traveltime. The total phase angle θ(f) including the non-dispersive part and the dispersive part of the total phase angle θ(f) is

θ(f) = 2πft + B(f) (1.29)

where t is the non-dispersive (pure) traveltime — the linear-with-frequency part from the "time-shifting theorem" — as the frequency component f travels the distance x in the time t.

Thus equation (1.29) can be written on the form

θ(f) = 2πft = 2πft (x/c(f) (1.30)

where t is the total time equal to the linear-with-frequency time plus the dispersive time. From (1.29) and (1.30) we have

2πft (x/c(f) = 2πft + B(f) (1.31)

Solving (1.31) for c(f) we get:

(1.32)

In the denominator of (1.32) we substitute

B(f) = bH(-f) (1.33)

We make one other substitution. The quantity τ in (1.32) is a non-dispersive traveltime. If we choose an upper frequency limit f_N above which there is assumed to be no dispersion then we can write

τ = x/(c/f_N ) (1.34)

The upper frequency limit can be conveniently taken to be the Nyquist frequency. It may be, as with Wuenchel's [200] experimental data for Plexiglas, that this will be a real data point on an actual dispersion curve. Note that we have switched from continuous notation for c(f) to discrete notation c(f_N ) because we will be taking the discrete Hilbert transform. Substituting (1.33) and (1.34) in (1.32) we get

= (1.35)

The quantity b in (1.35) is

b = the slope of α(f) versus f = πx/(Qc/f_N ) (1.36)

Making this substitution in (1.35) we get

C(f) =

And after some calculations:

(1.37)

Clearly, the phase velocity c[f] cannot depend on a sampling interval, which is implied by the value of the velocity c(f_N) at the Nyquist frequency f_N. What the velocity c_N does is to scale the dispersion curve without changing its shape; that is, c(f) scales the quantity

as defined by (1.37). As already noted, c(f_N) at f_N might be an actual observed coordinate on a dispersion curve as is the case for Wuenchel's Plexiglas data.

Equation (1.37) is Ecevitoglu's equation for the computation of the exact shape of the body wave phase velocity p[f] versus frequency f in a dispersive medium for any value of Q. Every value of dispersive velocity can be computed using (1.37) up to f = f_N; there are no arbitrary constant(s) to be chosen. Although the notation H(-f) looks like one is taking the Hilbert transform of the variable f it can be interpreted as taking the discrete Hilbert transform of a straight line defined by mf + b where the slope m =- 1 and the intercept b = 0.

If Q and one other phase velocity c(f) are known then c(f_N) can be computed from Equation (1.37). For example, Wuenschel's [200] experiment with the Pierre Shale yielded a coordinate {150 Hertz, 6835 ft/sec}. For Q = 30 using this pair in (1.37) and computing the value H(-f) =H(-150) gives

C(f_N ) = c(150)[1+(1/2x30) (165.008/150)]

= 6960.31

This value can also be obtained by trial-and-error, guessing at any value of c(f_N ) above 150 Hertz, and changing it until c[150] = 6835. The result will be c(f_N )= 6960. This is so because c(f_N ) is just a scale factor and does not influence the shape of (1.37).

If Q is unknown but c(f_N ) is known as well as any one coordinate {f, c(f)} then Q can be computed from (1.37). Wuenschel's experiment with the Pierre Shale yielded observed coordinates {150 Hertz, 6835 ft/sec} and {500, 6960}. Using these in (1.37) and computing the value H(-f) = gives

Q = (6835) (165.008) / 2(150)(6835 - 6960) = 30

It will be useful to display graphically the intermediate and final steps in the derivation of (1.37). A plot of —f versus f is shown in Figure 1.4 a, then H(—f) (Figure 1.4 b), and finally c(f) (Figure 1.4.c). First in the context of body wave dispersion, we have the equations:

(1.38)

(1.39)

The Hilbert transform of -f is obtained by first taking the inverse Fourier transform of -f, which is the outer integration in (1.39):

will be an even function in the time domain

Calculations

This outer integration will be an even function in the time domain because -f (or the dummy variable "-u") is even in the frequency domain. We do the inverse transform defined by the outer integral of (1.39). That is, in order to evaluate the outer integral of (1.39) an even function of —f was simply constructed as illustrated at the top of Figure 1.4.a.

Figure.1.4. (a) Plot of —f before Hilbert transformation, which is an even function by construction because we wish to use Mathematica's InverseFourier subroutine for computations,

We are operating on an even frequency-domain function so the inverse operation is really just an inverse cosine transform, as required by (1.39). Now we are in the time domain. We know that to obtain the causal function associated with — f (the real part of the Fourier transform of the unknown causal function) we must suppress the negative times of this even time-domain function (whatever it is, we don't know or care yet) and multiply the result by "2". This is exactly what |t| and multiplication by "2" do in (1.39). Now we have a causal function (that we have not yet inspected). Take the Fourier transform of this causal function. Equation (1.39) says to take the sine transform but using Fourier gives us that if we just look at the odd part of the transform of the causal function obtained by use of |t| and "2". The odd part of the resulting Fourier transform is shown in Figure 1.4.b. From (1.) the odd part of the Fourier transform is the Hilbert transform of the even part.

Figure.1.4.b: (b) H(-f), the Hilbert transform of —f, which gives an odd function,

Figures 1.4 a and b constitute a Hilbert transform pair. The transformation is the discrete version of Equation (1.26).

How do the above results for the computation of body wave dispersive velocities compare with those of others who approached the problem differently? In a classic paper Futterman [72] gave an approximation to the dispersive phase velocity c(f) versus f as

(1.40)

where c_q is the velocity at f₀ , an arbitrary low-frequency cutoff. His Equation (1.40) requires the choice of two constants c₀ and f₀ whereas (1.37) requires none. Following Futterman, c₀ is chosen below the frequency range of interest. Because the choice is arbitrary we can choose f₀ = 1 Hertz. Having chosen f₀, the constant c₀ can be chosen such that the theoretical dispersion curve defined by (1.37) passes through some known coordinate {f,c(f)} that has been determined from observed data.

As an example, Wuenschel chose a value of c₀ such that the dispersion curve defined by (1.37) passed through the point whose coordinates are {150 Hertz, 6835 ft/sec}, the latter values determined from observed data from the Pierre Shale. Equation (1.37) can then be used to determine the entire dispersion curve as shown by the upper thin solid line in Figure 1.4.c. Equation (1.37) on the other hand requires no arbitrary constants. It describes an exact shape. The value of c(f_N) is actually the dispersive velocity value that would be obtained experimentally for a rock of a given value of Q at the frequency f_N. The fact that f_N happens to be designated as the "Nyquist frequency" here is immaterial and is simply computationally convenient.

Figure.1.4. (c) Dispersive velocity c(f) as determined from Equation (1.37) for Q = 250 and c(f_N )= 5000 m/sec. Only velocity values for positive f need to be displayed.

The physical process of absorption between a source and a receiver is complete when the wavelet is recorded at the receiver. The absorption process is independent of any sampling interval and therefore of any value of f_N. In order to implement (1.37), however, we must assign a value of f_N and c(f_N)_.

Equation (1.37) requires no arbitrary constants (actually, not even a value for c(f_N)) to give the correct shape of the dispersion curve. It describes an exact shape. A value of a coordinate {f₀, c(f₀)} is required to compute Futterman's dispersion curve but the shape of the curve is approximate. Equation (1.37) can therefore be used to approximate velocities over a limited frequency range but not the shape of the entire curve unless many values of c₀ are used and it is determined piecewise. Futterman's approach acknowledges the observation by Aki and Richards [3] that it is not possible to take the Hilbert transform of a straight-line (the linear-with-frequency assumption for attenuation) that goes off to infinity and therefore does not converge. Futterman approached this problem by deriving an approximation that does converge but that predicts velocity c(f) versus frequency f only over a limited frequency interval.

Ecevitoglu and Costain [68] approached the problem by assuming that the linear-with-frequency variation was periodic over the interval f = 0 to f = f_N and that velocity values c(f) could be computed over the entire range of frequencies to give the correct shape of the dispersion curve. (Review the Dirichlet conditions on page 64 for the conditions under which a Fourier series—read "Hilbert transform"—converges to the function wherever it is continuous.) The assumption of periodicity results in a phase that is periodic and therefore bounded versus Aki and Richard's phase that is unbounded and increases without limit. The theoretical result of (1.37) is nicely supported by the Plexiglas dispersion experiment of Wuenschel [200] . Although c(f_N) in (1.37) can be thought of as a scale factor that plays the same role that c₀ does in Equation (1.37), it is more than that; in the case of (1.37) it is a real, physical velocity c(f_N) at the frequency f_N associated with a particular value of Q. We compare numerical values from each method below.

Fig.1.5 Summary of equations for absorption-dispersion pairs.Frequency is denoted ν in the figure (f in the text).

See Figure 1.5 for a summary of differences between the two approaches to body wave dispersion as well as those of others.

Comparison of dispersive phase velocity values

The theoretical dispersion curves of Futterman [72] and Ecevitoglu and Costain [68] both require prior knowledge of a single {velocity, frequency} coordinate through which, in the case of (1.37), the exact shape of the dispersion curve will pass. As an example, c(f_N) can be chosen by trial and error such that p[150] passes through Wuenschel's observed coordinate of {150 Hertz, 6835 ft/sec}, as shown in Figure 1.6. This value of c[f_N], once chosen, is the actual dispersive rock velocity that would be obtained if we were to make a velocity measurement at that frequency. Once c(f_N) has been determined the entire exact body wave dispersion curve can be computed from (1.37) instead of from (1.39).

Using Wuenchel’s results for the Pierre Shale and a value of Q = 30 and selecting a value of f_N = 500 Hertz (clearly f_N = 500 from his figure) then we can compute the entire theoretical dispersion curve. The value of c(f_N ) is changed by trial and error until it is found to be equal to 6960 ft/sec when f = 150 Hertz and c[150] = 6835 ft/sec. Now, with c[f_N] set equal to 6960 ft/sec in (1.37) we have the dispersion curve for all frequencies from f = 0 to f_N = 500. Wuenschel set the constants {f₀, c₀} in (1.39) equal to {1 Hertz, 6522 ft/sec} so that the theoretical curve from (1.39) can approximate the observed data. From (1.37) we compute the pair {1, 6477} in agreement with Wuenschel to within about 0.6 %.

As a comparison with (1.39) take several values for Q = 30 obtained from (1.37) and compare them with those obtained from (1.39). The values of {frequency, velocity} obtained using (1,37) with c(f_N) = 6960 are:

{10., 6645. 54}, { 20. ,6692.48}, { 30. ,6720.3 }, {40. , 6740.21},

{50., 6755. 77}, { 60. ,6768.56}, { 70. ,6779. 44}, {80. , 6788.93},

{90., 6797. 35}, {100. ,6804.93}, {110. ,6811 .83}, {120. ,6818. 18},

{150 ., 6834 . 69} , {200 . , 6856 . 74} , {300 . , 6890 . 78} , {400 . , 6920 . 47}

These values can be compared with Futterman's equation after choosing f₀ = 1 Hertz and c₀ = 6522 ft/sec, the latter value having been chosen as described above.

Additional comparisons with Futterman are shown graphically in Figure 1.6 from Ecevitoglu and Costain [68] where Futterman's [72] velocity dispersion expression is superimposed in a piecewise manner upon the curve denned by Equation (1.40). The excellent agreement using constants c_q and f₀ for Futterman obtained from 1.37) is apparent in the figure. Note that it is not possible to superimpose Futterman's entire results with a single selection of his constants. On the other hand, there are no arbitrary constants that govern the shape of the velocity dispersion curves computed from discrete Hilbert transformation. This leads us to Figure 1.6 and generalized dispersion curves.

As noted earlier several attempts have been made to formulate analytical expressions for absorption and the body wave dispersion that must accompany absorption. It is beyond the scope of this volume to delve into the details of this history. An excellent review can be found in [66]. Aki and Richards [3, pages 172-177] summarized difficulties with analytic expressions that involve the total phase (linear-with-frequency phase plus dispersive phase) and frequency limits of integration that extend to infinity. Ecevitoglu and Costain [68] suggested that the problems of a divergent integral as described by them [3, Equation 5.74] can be overcome by invoking discrete numerical Hilbert transformation and by recognizing that the total phase spectrum can be split up into a linear-with-frequency nondispersive phase defined by the traveltime of a reflected event plus a pure dispersive phase spectrum that is associated with the body wave dispersion required by causal absorption. A summary of the phase definitions of Futterman [72], Strick [172], Kjartansson [94] and the later results of Ecevitoglu and Costain [68] that have been discussed here is shown in fig.1.5. For additional details the reader is referred to Ecevitoglu [66].

This completes our discussion of body wave dispersion, which is one member of the "absorption-dispersion pair". We now turn our attention to absorption, the other member of the pair, and its effect on the shape of the wavelet.

Figure 1.6 Upper light line is the dispersion curve in Pierre Shale (dashed line) from Wuenschel [200] computed using Equation (1.40) and Q = 30. Lower heavy line is the exact dispersion curve from Equation (1.40) using the discrete Hilbert transform. The curves are not identical in shape because Futterman's Equation (1.40) is an approximation. The solid line could just as well have been scaled by the value p[vn] so that it passed through the dot that represents the coordinate of the observed data, as shown in the bottom diagram. But then the solid line would not have agreed with (1.40) elsewhere.

Absorption

Wuenschel [200] noted that for purposes of seismogram synthesis, a few percent variation in phase velocity can produce a significant change in pulse waveform for long propagation distances. The pulse waveform can be thought of as a vernier to magnify small differences in phase velocity.

If only absorption is taken into account without the B[f] causal phase term then the initial wavelet spreads into negative as well as positive time as shown earlier in Figure 1.6, and therefore violates causality. With the addition of the B[f] term the lower frequency sinusoids are required to travel more slowly in just the right amount to keep the wavelet causal.

The plots in Figure 1.4 a and b constitute a Hilbert transform pair. It is important to note that the shape of Figure 1.4 b is actually the shape found by Ecevitoglu [66, Figure 27, p. 95] from real conventional CDP reflection data analyzed by singular value decomposition. His results might be the first to demonstrate causal absorption from conventional CDP reflection data. It is also apparent that if absorption is not linear with frequency then this will show up in the analysis because the shapes of Figures 1.4 a and b are Hilbert transforms of each other. For example there might be a "hump" or a local depression on the linear-with-frequency curve that will also affect the causal phase.

In our absorption Earth model (absorption only, no spherical spreading for example) the source emits an impulse δ(t) that contains all frequencies f and all of equal amplitude. The amplitudes of frequency components from this impulsive source are gradually attenuated by absorption along the travel path a. This results in a lengthening of the impulse δ(t) until after a distance x = 1500 m, for example, in an absorptive medium with a quality factor Q = 20. The Dirac delta function has lengthened and now looks like that shown in Figure 1.7.

Figure 1.7. (a) The absorption impulse response of the Earth for Q = 20 and a = 1500 m. This is a plot of Equation (1.41).

As discussed above, this absorption impulse response h(t) of the Earth to an impulsive plane wave is

(1.41)

Here B(f) is the Hilbert transform of –bf and from (1.8.b)

B= π/(cQ)

We could convolve h(t) with a seismic wavelet to get the new causal shape at the distance x. Alternatively we can multiply the Fourier transform of h(t) by the spectrum of the source wavelet at x = 0. This is done on fig.1.8.b

X=1000 Q=30

X=100 Q=10

X=250 Q=10

X=250 Q=30

X=500 Q=10

X=800 Q=10

Figure 1.8: Wavelet spreading for Q = 30 for various distances x=0 and x=. For this example the "wavelet" is a rectangular function of width 2 ms. Wavelet broadening increases for increasing values of the offset distance a but all wavelets remain causal.

Gladwin and Stacey [78] used the theories of Futterman [72], Strick [171, 172] and Azimi et al. [7] to calculate the shapes of impulse responses after propagation in a medium of Q = 50 to a distance of x = 5 m. Their results for Futterman's theory are shown in Figure 1.9 along with the shape computed from Equation (1.40). The comparison is excellent. The different arrival times shown in their Figure 4 for Futterman [72], Strick [171, 172], and Azimi et al. [7] are due to the "hidden" phases summarized in Fig.1.5 and discussed in detail by Ecevitoglu [66]. As noted by Gladwin and Stacey [78] the use of Futterman's theory by them resulted in a waveform that arrives too early. Ecevitoglu [66, p. 7 and Figure 2] explained this phenomenon by showing that a second term (a "hidden" phase term resulting from the approximation in Futterman's theory) in the phase angle spectrum subtracts from the linear-with-frequency term and causes the early arrival.

Normalized dispersion D

The normalized amount D of body wave dispersion, hereafter referred to as just "dispersion", over the bandwidth from f = 0 to f = f_N can be defined as

D = (c(f_N ) – c(0)) / c(f_N ) (1.42)

The Nyquist frequency component f_N and the frequency component at f = 0 both travel the same distance x at their respective phase velocities. So

c[f_N] τ = c[0] [τ + ∆W] = x

where "c" denotes phase velocity and ∆W is a dispersive time delay due to body wave dispersion. Solving for ∆W

∆W = [(c(f_N ) – c(0)) / c(f_N )] (1.43)

From (1.37) and taking into account limiting values of H(-f)/f for f = 0 and f = f_N the phase velocity at f = 0 is

C(0) = c(f_N)/(1+2/πQ) (1.44)

Repeating (1.42) and substituting (1.44) gives

D =( c[f_N] – c[0])/c[f_N] = 1/(1+π/2 Q) (1.45)

which says that the maximum amount of velocity dispersion for a given value of Q is independent of any frequency f, even f_N. For example, if Q = 20 then

D = 1/ (1+π/2 20 ) = 0.031

which says that the maximum amount of body wave dispersion between f= 0 and f =∞ for Q = 20 is 3.1 %. Values of D over a more restricted range of frequencies are easily determined by measurements in the frequency domain as discussed below.

Comparison of dispersion D values with real data

It is possible to estimate the amount of f₂ using

D = (c(f₂ ) – c(f₁)) / c(f₂ ) (1.46)

This formula can be applied to several excellent data sets such as the pulse-propagation experiment in the Pierre Shale by McDonal et al. [118]. Wuenschel [200] reported an observed variation in phase velocity of 2.6 % between 25 and 400 Hertz for the Pierre shale, to which he assigned a value of Q ≈ 30. According to (1.46) this value of Q corresponds to D = 2.1%. Waters [187, p. 33] assigned an average value of Q = 17 to samples of the Pierre shale for phase velocity determinations over the frequency range 75 — 555 Hertz. According to (1.45) a value of Q = 17 corresponds to D = 3.6 %. Considering the apparent range of values of Q for the Pierre shale the values of total dispersion from (1.45) appear to be in excellent agreement with published values over more restricted frequency ranges.

Wuenschel [200] published a composite dispersion curve for Plexiglas. From his experimental data he observed that for Plexiglas the phase velocities vary by 2.3 % over a five-octave range from 4 to 128 kilocycles per sec. Using (1.40) it is possible to estimate a value for D over this same frequency range by setting c[f_N] equal to 7720 ft/sec, which is the observed value of velocity that can be visually estimated from his Figure 12. Wuenschel gives an estimate of 55 for Q.

A time-domain method for determination of Q

Body wave dispersion and absorption result in a change of wavelet shape. Glad-win and Stacey [78] and Kjartansson [94] showed that wavelet broadening due to absorption and body wave dispersion is proportional to the wavelet traveltime and is related to Q by the empirical equation

W₂ = (C/Q)(τ₂ –τ₁) + W₁ (1.47)

where τ₁ is the wavelet traveltime on one trace and τ₂ is the traveltime to, say, an adjacent trace. W₁ is the pulse width at traveltime τ₁, W₁ is the pulse width at time τ₂ , and C is an arbitrary constant to be determined empirically from some other means of measuring Q such as spectral ratio techniques. Ecevitoglu [66, p. 23] pointed out that there is no consensus in the literature about the exact value of C. Kjartansson [94] used C = 0.485 and Gladwin and Stacey [78] and Badri and Mooney [9] used C = 0.5. Ecevitoglu [66] derived an exact value of C as

C=2/π

Ecevitoglu's [66] exact equation corresponds to the empirical Equation (1.47) of Gladwin and Stacey [78] and Kjartansson [94]. From (1.43) and (1.44) we have

∆W = (( c[f_N] – c[0])/c[0]) τ =

= D ( c[f_N])/c[0]) τ

=(2/πQ) τ

Or

Q= 2τ/π∆W (1.48)

This is Ecevitoglu and Costain's [66, 68, Equation 28, p. 23] result for the determination of Q in the time domain. See Figures 1.9 and 1.10 for an application. The quantity τ is the total traveltime during which body wave dispersion in the amount of ∆W = W₂ — W₁ has taken place. For our model example W₁ is the pulse width of the wavelet at the source at x = 0. W₂ is the pulse width after the wavelet has traveled for the time τ. For values of W₁ and W₂ one measures the same peak to peak, trough to trough, or peak to trough time duration even though the wavelet is now distorted by causal absorption. Values of τ and Wi will normally have to be in fractions of a millisecond.

Fig.1.9

The delay of the absorbed wavelet of about 10 ms shown in Figure 1.9 is an artifact of the processing and does not affect the results. The delay is due entirely to the choice of the model sampling interval of 0.1 ms in order to facilitate accurate picking of the values of Wi and Wz. The small sampling interval of the modeled data implies high frequencies that are not present in the wavelet so these higher frequencies do arrive before the lower frequencies, as they should, but they have essentially zero amplitude. The result is the 10 ms "delay". If actual field data were to be recorded at a sample interval of 0.1 ms (same as the modeled data) then the attenuated wavelet would be shifted to the left and the "delay" would disappear. If the wavelets modeled in Figure 1.9 were instead generated using a sampling interval of 4 ms then the "delay" would disappear but it would not be possible to pick accurate times for the peaks or troughs. On the other hand, the resolution could be increased to 0.1 ms.

Figure 1.10: Ricker wavelet shape change for Q = 20, x = 300 m, p[c _n] = 2000 m/sec. Sampling interval dt = 0.1 ms. Times shown are in ms. Solid line is initial wavelet. Dashed line is attenuated and dispersed Ricker wavelet at a source-receiver distance of 300 m. Value of Q recovered from the analysis using Equation (1.48) is 21.2 versus model value of 20.

If a field study is designed to examine body wave dispersion in low-Q rocks then a small sampling interval, say 0.1 ms, should be used. If the data are already recorded at, say 4 ms, then (8.139) can be used to resample at 0.1 ms with equivalent results to search for body wave dispersion.

It is clear that Ecevitoglu's time-domain approach to the measurement of Q has several advantages. For example, it is not necessary to see the entire wavelet. Under ideal conditions only a single peak or trough and another adjacent to it are necessary. So overlapping wavelets and windowing effects become less important. The redundancy of reflection seismic data means that the occasional noisy trace can be tolerated. Q can be determined using single CDP gathers and changes monitored horizontally and vertically. Ecevitoglu used singular value decomposition along with a great multiplicity of data to estimate Q in a crustal volume using (1.48) as well as to produce what might be the first measured observation of the phase angle spectrum associated with body wave dispersion. From this, the Hilbert transform yielded the attenuation coefficient as a function of frequency.

Figure 1.11: Klauder wavelet shape change for Q = 20, c(f_N)= 1500 m/sec at a distance of x = 100 meters. ∆t = 0.1 ms. Solid line is initial Klauder wavelet. Dotted line is attenuated and dispersed Klauder wavelet. The value of Q recovered from the analysis using Equation (8.138) and using a search algorithm for troughs is 25 versus the model value of 20. Decreasing W\ by 0.1 ms and increasing Wiz by 0.2 ms results in a value of Q = 20.2

Figure 1.12: Klauder wavelet shape change for Q = 10, c(f_N)= 1500 m/sec at a distance of x = 100

Summary

The results described in this section for the computation of body wave dispersion are in agreement with those of earlier workers but the approach is more general because no arbitrary constants are required to determine the shape of the body wave dispersion curve at any frequency, thus defining "the" pair of absorption-dispersion equations as discussed on page 519 and summarized in Table 8.1. Athough a linear-with-frequency attenuation is used here as an example the numerical approach described by Ecevitoglu [66] is appropriate for any behavior of o(f) versus v. This has not been emphasized in the above discussion except for Figure 8.83 but it means that the inverse problem, that of determining a(v} versus v from Hilbert transformation of the phase spectrum, the phase spectrum being derived from real data, will reveal the nonlinear dependence of a(t') versus v if it is present in the data. This should have important implications for detecting frequency ranges over which values of Q are different but constant over limited segments.

Figure 8.81: Generalized body wave dispersion curves p(f) for Q = 5, 10, 15, 20. p[f] /p[f _n] determined from Equation (1.40). Frequency v normalized by division by l/(2∆t) =f_N where ∆t is the sampling interval. The quantity c[f_N] is the rock velocity at the Nyquist frequency, not at f = ∞. If Q is known it is not necessary to know c[f_N]. In general a control point {f,c(f)] will be used to determine c[f_N ] as per the example given for the Pierre Shale.

Figure 8.82: Solid, thinner continuous line is the result of using Equation (8.127) to compute phase velocity p[i/] versus i/. Short and thicker discontinuous line segments superimposed on the continuous curve are computed from Putterman's [72] Equation (8.129) for phase velocity, which is shown on the figure. The constants cq and vq used in Futterman's equation are taken from Equation (8.127). The agreement is excellent but the use of Futterman's equation requires a choice of arbitrary constants Co and v$ whereas the solid line resulting from the use of discrete Hilbert transformation "H[—v\ is exact and requires no arbitrary constants. Comparison made for Q ~ 250. Figure from Ecevitoglu and Costain [68].

Estimates of fractional values can be obtained from the Mathematica graphics output of program Absorption. nb by holding down the "Ctrl" key while moving the cursor over the graphics output to the desired locations shown in Figure 8.87. A better approach might be to use

(8.139)

n=l

for an interpolation to a finer sampling interval and then use (8.138).

In (8.138) the time r is to be compared with the difference between travel-times ti and T2 in Gladwin and Stacey's Equation (8.137). That is, for Gladwin and Stacey, dispersion has taken place during the time ti — r\ whereas in (8.138) dispersion has taken place during the entire time r. That is, in Equation (8.137) t! is zero and T₂ in (8.137) is then the same as t in (8.138). So (8.137) and (8.138) are the same except for the constant C in their equation, which was found by Ecevitoglu to be

c=²-

7T