Measuring Velocity Dispersion And Attenuation

GEOPHYSICS, VOL. 74, NO. 2 共MARCH-APRIL 2009兲; P. WA113–WA122, 20 FIGS. 10.1190/1.3068426

Measuring velocity dispersion and attenuation in the exploration seismic frequency band

Langqiu F. Sun1, Bernd Milkereit1, and Douglas R. Schmitt2

alter perfect elasticity. For example, Zener 共1948兲 presented the standard linear solid model, in which the internal friction Qⳮ1 共where Q is quality factor; Knopoff and MacDonald, 1958兲 shows a peak at a certain frequency. This type of attenuation model is referred to as the Debye peak model. By investigating some attenuation models, Knopoff and MacDonald 共1958兲 concluded that Q and wave velocity are frequency dependent in a linear medium, i.e., one in which stress is related linearly to strain. The strain induced by passing seismic waves is small, and the earth materials can be reasonably assumed to be linear. However, attenuation observations 共e.g., Knopoff, 1964兲 support a constant-Q model; that is, Q is constant in a broad frequency band 共f L through f H兲. Liu et al. 共1976兲 explained the contradiction between observation and theory by superposing a series of Debye peaks to produce a nearly constant-Q model 共Figure 1兲. The causality of seismic waves requires seismic velocity to be frequency dependent 共i.e., velocity dispersion exists兲 in a medium with attenuation. Furthermore, the linkage between velocity dispersion and attenuation is the Kramers-Krönig relation 共Futterman, 1962兲. Velocity dispersion and Q can be written explicitly as 共Bourbié et al., 1987兲

ABSTRACT No perfectly elastic medium exists in the earth. In an anelastic medium, seismic waves are distorted by attenuation and velocity dispersion. Velocity dispersion depends on the petrophysical properties of reservoir rocks, such as porosity, fractures, fluid mobility, and the scale of heterogeneities. However, velocity dispersion usually is neglected in seismic data processing partly because of the insufficiency of observations in the exploration seismic frequency band 共⬃5 through 200 Hz兲. The feasibility of determining velocity dispersion in this band is investigated. Four methods are used in measuring velocity dispersion from uncorrelated vibrator vertical seismic profile 共VSP兲 data: the moving window crosscorrelation 共MWCC兲 method, instantaneous phase method, time-frequency spectral decomposition method, and cross-spectrum method. The MWCC method is a new method that is satisfactorily robust, accurate, and efficient in measuring the frequency-dependent traveltime in uncorrelated vibrator records. The MWCC method is applied to the uncorrelated vibrator VSP data acquired in the Mallik gas hydrate research well. For the first time, continuous velocity dispersion is observed in the exploration seismic frequency band using uncorrelated vibrator VSP data. The observed velocity dispersion is fitted to a straight line with respect to log frequency to calculate Q. This provides an alternative method for Q measurement.

V共␻ 兲 ⳱

冑

2兩M共␻ 兲兩2 , ␳ 共兩M共␻ 兲兩 Ⳮ M R共␻ 兲兲

and

Q共␻ 兲 ⳱

INTRODUCTION The energy of seismic waves propagating in an anelastic medium is attenuated by various dissipation mechanisms. This phenomenon has been observed in experiments on different solids, liquids, and seismograms 共e.g., Gutenberg, 1958兲. Early attempts to account for this phenomenon include various modifications of Hooke’s law to

M R共 ␻ 兲 , M I共 ␻ 兲

where ␳ is density, ␻ ⳱ 2␲ f is angular frequency, and M is the complex elastic modulus defined by the ratio of stress to strain. The real part of M, M R, and the imaginary part of M, M I, constitute a Kramers-Krönig pair:

M R共␻ 兲 ⳱ M R共0兲 Ⳮ

2␻ 2 ␲

冕

⬁

0

M I共 ␣ 兲 d ␣ ␣ ␻2 ⳮ ␣ 2

Manuscript received by the Editor 20 January 2008; revised manuscript received 15 September 2008; published online 11 March 2009. 1 University of Toronto, Department of Physics, Toronto, Ontario, Canada. E-mail: [email protected]; [email protected]. 2 University of Alberta, Department of Physics, Institute for Geophysical Research, Edmonton, Alberta, Canada. E-mail: [email protected]. © 2009 Society of Exploration Geophysicists. All rights reserved.

WA113

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WA114

Sun et al.

M I共 ␻ 兲 ⳱

2 ␲

冕

⬁

共M R共␣ 兲 ⳮ M R共0兲兲

0

␻ d␣ ␣ ⳮ ␻2 2

Thus, in theory, Q can be calculated from velocity dispersion. This alternative method of determining Q was deemed impractical because velocity dispersion is very difficult to measure in reflection seismograms, especially with noise 共Jannsen et al., 1985兲. In a constant-Q model, seismic velocity increases with frequency 共Kjartansson, 1979兲:

冉冊 冉冊

V共f 2兲 f2 ⳱ V共f 1兲 f1 where

␥⳱

␥

,

1 1 tanⳮ1 , ␲ Q

where f 1 and f 2 are in the f L through f H frequency band. For example, if Q ⳱ 20, and the seismic velocity is 2.6 km/s at 20 Hz, then from equation 1, the velocity is 2.7 km/s at 200 Hz, increased by 4%. When 兩␥ · ln共f 2 /f 1兲兩
1 and tanⳮ1共1/Q兲 ⬇ 1/Q, using Taylor’s expansion, the above equation becomes

1 V共f 2兲 f2 ⬇1Ⳮ ln . V共f 1兲 ␲Q f1

共1兲

In the exploration frequency band, 兩log共f 2 /f 1兲兩
4, hence equation 1 is valid when Q ⬎ 10. This type of velocity dispersion is referred to as linear velocity dispersion 共see Figure 1兲. Equation 1 is identical to the velocity dispersion in Liu et al. 共1976兲 for a nearly constant-Q model. The later research on attenuation and velocity dispersion includes modeling studies on various petrophysical conditions, e.g., partial gas saturation 共White, 1975兲, squirt flow 共Mavko et al., 1998兲, patchy-saturated porous media 共Johnson, 2001兲, double-porosity dual-permeability materials 共Pride and Berryman, 2003兲, and fractured porous media with fluid flow 共Brajanovski et al., 2006兲. These modeling studies indicate that attenuation and velocity dispersion in saturated porous rocks depend on petrophysical properties, e.g., porosity, fractures, and fluid fill. Velocity dispersion in the sonic and ultrasonic frequency band is observed in rock experiments 共e.g., Spencer, 1981; Winkler, 1986; Jones, 1986; Adam et al., 2006; Batzle et al., 2006兲. Velocity dispersion in the seismic and sonic frequency band is observed in field surveys 共e.g., Brown and Seifert, 1997; Sams et al., 1997兲. These obser-

Figure 1. The constant-Q model resulted from superposed Debye peaks and the corresponding velocity dispersion. Here f L and f H denote the frequency range in which Q is constant. After Liu et al. 共1976兲.

vations demonstrate that velocity dispersion could be a strong indicator of petrophysical properties. However, most velocity dispersion data are in the sonic and ultrasonic frequency band, and are sparse in the seismic frequency band 共⬃10 through 200 Hz兲. Therefore, the existing velocity dispersion observations are insufficient because they do not provide enough information in the seismic frequency band, which is of greater interest in exploration seismology because the range of seismic wavelength is a few meters through a few hundred meters, making it possible to assess the bulk rock properties in a rock volume. Velocity dispersion usually is neglected in conventional seismic data processing because the effect is small and difficult to measure in a medium with Q ⱖ 30 共Futterman, 1962兲. However, in high-attenuation media with Q ⬍ 30, velocity dispersion is not negligible 共Molyneux and Schmitt, 1999兲. A frequency-dependent phase shift will be introduced into the seismic data because the traveltime varies with frequency. Figure 2 shows how attenuation and velocity dispersion distort the conventional crosscorrelation process for vibrator data. In the absence of velocity dispersion, the correlation wavelet is zero phase. Assuming a constant Q and linear velocity dispersion in the frequency band of the pilot vibrator sweep, as the sweep propagates, the high-frequency component is attenuated quickly as a result of the constant-Q attenuation, and the correlation wavelet becomes mixed phased because of velocity dispersion. In summary, attenuation and velocity dispersion can provide new insight into physical rock properties. However, present observations in the exploration seismic frequency band are insufficient to make any valuable deductions. In addition, seismic data processing can be

a)

b)

Figure 2. Distortion of correlation wavelets of vibrator sweeps in a medium with attenuation and velocity dispersion. The input vibrator sweep is linear, 12 s long, 8–180 Hz, with 0.15 s tapering. 共a兲 The Q and corresponding velocity dispersion model in the 8–180-Hz band. 共b兲 Correlation wavelets of the vibrator sweeps changing with the source-receiver distance.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

Velocity dispersion in the seismic band affected severely if attenuation and velocity dispersion are intense. Hence, it is necessary to develop a robust method to measure velocity dispersion in field seismic data.

Methodology Suitable data type: Uncorrelated vibrator data Attenuation alters the shape of a signal’s amplitude spectrum, whereas phase velocity dispersion changes the phase spectrum. Therefore, uncorrelated vibrator data are appropriate to measure the small velocity dispersion in the seismic frequency band. The advantages of using uncorrelated vibrator sweeps include 共1兲 the ability to control or to measure both the amplitude and phase spectra of the pilot sweep, and 共2兲 the retention of both the amplitude and phase spectra of the received sweep, which is no longer possible once the signal has been correlated. In a vibrator survey, seismic energy of varying frequencies is launched from the source and arrives at the receiver at different times, forming a time-frequency 共t-f兲 relation. This t-f relation of a pilot sweep 共pilot t-f relation兲 is predetermined in a vibrator survey. Figure 3 shows the predetermined t-f relation of a pilot sweep, the pilot sweep, received sweep, and their t-f spectra. Different events, such as the direct wave, reflections, harmonics, and noise, are separated satisfactorily in the t-f spectra, making it possible to analyze the t-f relation of a particular event. The most visible event in a received t-f spectrum usually is the direct wave. Only the direct wave is discussed in this paper; “received t-f relation” refers to the t-f relation of the direct wave in a received sweep. When velocity dispersion is negligible, the received t-f relation will be parallel to the pilot t-f relation. Otherwise, the received t-f relation will deviate. In a constant-Q attenuation model, seismic velocity is larger at higher frequencies; hence the higher-frequency component always arrives earlier than expected. Although the physics of the problems differs, in many ways this frequency dispersion could be compared with the Doppler effect in marine vibrator surveys using a moving boat 共Dragoset, 1988兲. Thus, if a vibrator sweep is distorted by velocity dispersion, the distortion can be corrected with the similar method presented by Dragoset. The difference between the traveltimes at a high frequency and at a low frequency ⌬t can be used to determine velocity dispersion. The ⌬t usually is very small in the seismic band. For example, if f 1 ⳱ 20 Hz, f 2 ⳱ 200 Hz, V共f 1兲 ⳱ 2.6 km/s, V共f 2兲 ⳱ 2.7 km/s, and the source-to-receiver distance is 1 km, then ⌬t is only 14 ms in the f 1 through f 2 band. From equation 1, ⌬t increases as f 1 /f 2 increases. In addition, ⌬t increases if the source-receiver distance increases. Therefore, to optimize observability of the small velocity dispersion in uncorrelated vibrator sweeps, broadband, long-baseline data are desirable. As for data acquisition geometry, vertical seismic profile 共VSP兲 data are preferred because transmission seismograms are easier to analyze than reflection seismograms. Moreover, the method to measure velocity dispersion in the seismic frequency band needs to be accurate and robust for field data.

Methods of measuring velocity dispersion in uncorrelated vibrator data There are various methods of measuring velocity dispersion in uncorrelated vibrator sweeps. A desirable method must be 共1兲 accu-

WA115

rate because velocity dispersion is small in the exploration frequency band, 共2兲 robust in the presence of noise, and 共3兲 capable of separating different events because a field vibrator record usually is contaminated by noise and contains multiple events 共e.g., reflections, tube waves, and harmonics兲. The cross-spectrum 共CS兲 method 共Donald and Butt, 2004兲, instantaneous phase 共IPH兲 method 共Taner and Sheriff, 1977兲, t-f spectral decomposition 共TFSD兲 method 共Castagna and Sun, 2006兲, and moving window crosscorrelation 共MWCC兲 method have been tested using synthetic data. The MWCC method is presented below. The other methods and synthetic tests are discussed in the appendix. The MWCC method measures the t-f relation of an uncorrelated vibrator sweep in the time domain. The difference between the pilot and received t-f relations provides the frequency-dependent traveltime, T共f兲. Dispersion of this raypath-average velocity can be calculated by dividing the source-receiver distance by T共f兲. The procedure consists of six steps: Step 1兲 Take a part of the pilot sweep with a tapering window. Denote f i as the central frequency of this part of the pilot sweep. The window must be longer than a few wavelengths so that the correlation in the next step is valid.

b)

a)

d)

e)

c)

f)

Figure 3. An example of uncorrelated vibrator sweeps. 共a兲 The designed t-f relation of a pilot sweep; 共b兲 the pilot sweep; 共c兲 t-f spectrum of the pilot sweep; 共d兲 a received sweep; 共e兲 t-f spectrum of the received sweep; 共f兲 a schematic illustration of the events in 共e兲. The grayscale in 共c兲 and 共e兲 is amplitude. The spectral decomposition method for 共c兲 and 共e兲 is discrete Fourier transform using a 1-s-long cosine window. The data are from the 3L-38 Mallik gas hydrate research well, MacKenzie Delta, Northwest Territories, Canada. The receiver is in the borehole at a depth of 1085 m. The vibrator-borehole offset is 22 m. The sampling interval is 2 ms.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WA116

Sun et al.

Step 2兲 Crosscorrelate this part of the pilot sweep with the entire received sweep to produce a wavelet of correlation coefficient. Step 3兲 Calculate the envelope of the correlation wavelet so that a peak is produced at the arrival time of an event. This is the arrival time of the frequency component centered at f i. The most significant peak usually is the direct wave. Step 4兲 Move the tapering window along the pilot sweep for another f i, and repeat steps 1 through 4 so that a pseudo t-f spectrum is produced for an uncorrelated vibrator sweep. Similarly to a t-f spectrum from spectral decomposition, different events can be separated satisfactorily in a pseudo t-f spectrum, and the most significant event usually is the direct wave. Step 5兲 Pick the arrival time of the direct wave in the pseudo t-f spectrum for different f i to obtain the received t-f relation. For high-quality data, the arrival time can be picked continuously at each frequency component using an automatic arrival picking algorithm. Otherwise, the arrival time must be picked manually at a few control frequencies. Step 6兲 Calculate the difference between the pilot and received t-f relations to obtain the T共f兲.

provides satisfactory results in the presence of different types of noise and multiple events, but it requires the pilot t-f relation to be obtained from another method. Combining the IPH and MWCC methods provides the best resolution because the IPH method can measure the pilot t-f relation, and the MWCC can be used to measure the received t-f relation. The CS and TFSD methods can be applied for comparison and quality control.

Field data examples The Mallik gas hydrate field, located in the MacKenzie Delta, Northwest Territories, Canada, is of historic interest to gas hydrate investigators and a site of ongoing research since 1971 共Dallimore et al., 2005兲. There has been abundant research on the petrophysical and geophysical properties of this area. For example, Lee 共2002兲 models the P- and S-wave velocities in gas hydrate-bearing sediments using the Biot-Gassmann theory. Winters et al. 共2005兲 presents profiles of different petrophysical properties. In general, the zone of interest within the upper kilometer of sediment consists of permafrost, water-saturated sediments, and gas hydrate-saturated sediments.

Steps 1 through 3 are shown in Figure 4, and steps 4 through 6 are shown in Figure 5. The uncorrelated vibrator data in Figures 3–5 were acquired in the Mallik gas hydrate research well 共Dallimore et al., 2005兲 at a depth of 1085 m and offset 22 m. It must be noted that the MWCC method by itself cannot determine the central frequency f i of any part of a pilot sweep. This method determines only the arrival time of a wave pattern from the pilot sweep. The f i must be determined by other methods, for example, the IPH method. Although a pseudo t-f spectrum from the MWCC method resembles a t-f spectrum from spectral decomposition, they are different. The MWCC method is a time-domain method, and the values in a pseudo t-f spectrum are envelopes of correlation coefficients. A peak in a pseudo t-f spectrum implies the maximum likelihood of the arrival of an event, whereas a peak in a t-f spectrum represents high amplitude at a certain frequency. According to the results of synthetic tests 共see Appendix A兲, the CS method is fast, but it fails when noise or multiple events exist in received vibrator sweeps. The IPH method is fast and suitable for measuring the t-f relation of pilot sweeps, but it is not appropriate for received sweeps. The TFSD method is the most expensive and cannot separate multiple events as satisfactorily as the MWCC method. The MWCC method is less expensive than the TFSD method and

Figure 4. Steps to calculate the arrival time of part of a pilot sweep with central frequency f i using the MWCC method. The dot in step 1 stands for multiplication of the tapering window and pilot sweep. The asterisk in step 2 stands for crosscorrelation of this part of the pilot sweep and the entire received sweep.

Figure 5. The pseudo t-f spectrum of the uncorrelated received sweep shown in Figure 3d from the MWCC method using a 2-s-long cosine window. The pseudo t-f spectrum is obtained by repeating the steps shown in Figure 4 for different f i. 共a兲 Displayed in grayscale, which represents the envelope of the correlation coefficients. 共b兲 Zoomed in on part of the pseudo t-f spectrum and displayed in wiggles. The dashed line shows the t-f relation of the pilot sweep. The crosses show the arrival time picks of the direct wave at different frequencies.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

Velocity dispersion in the seismic band A three-component multioffset vibrator VSP survey was conducted in the 3L-38 Mallik gas hydrate research well in 2002. The uncorrelated vibrator data from this survey are analyzed to determine the seismic attenuation and velocity dispersion in this area. For the vibrator data set discussed below, the source-borehole offset is 22 m. The receivers are at 15-m intervals in a 560–1145-m depth range. The sampling interval is 2 ms. The pilot sweep is linear, 14 s long, 8–120 Hz above 935 m, and 8–180 Hz at, and deeper than, 935 m. Figure 6 shows the vertical component of the correlated VSP data, and the raypath-average P-wave velocity derived from well logging. The uncorrelated pilot and received sweeps are analyzed with the IPH and MWCC methods, respectively. Arrival time of the direct wave is picked continuously at each frequency in the pseudo t-f spectrum and compared with the pilot t-f relation to obtain T共f兲. As an example, T共f兲 measured in the vibrator sweep recorded at the depth of 1085 m is shown in Figure 7. The seismograms and pseudo t-f spectrum for Figure 7 are shown in Figures 3 and 5, respectively. Figure 8 shows the dispersion of the raypath-average velocity, indicating that, in general, seismic velocity increases with frequency. This trend can be explained with the linear velocity dispersion in a constant-Q model. By obtaining a best-fit straight line with respect to log frequency, the linear velocity dispersion is determined to be 6% in the 8–180-Hz band. Using equation 1, a raypath average Q of 15 is obtained.

a)

WA117

The above procedure is applied to all uncorrelated sweeps from the 3L-38 well at the 22-m source-borehole offset. Figure 9 shows the T共f兲 profile for all depth levels. For convenience, only the 8– 120-Hz band is shown. The trend of the high-frequency component being faster is clear. Superposed on this trend are narrow bands with earlier or later arrivals, a phenomenon that is systematic in this data set and inexplicable with constant-Q attenuation. Figure 10 shows the linear velocity dispersion calculated from the T共f兲 profile, com-

Figure 8. Solid line: dispersion of the raypath-average P-wave velocity in the vibrator sweep shown in Figure 3d, calculated from the T共f兲 shown in Figure 7. Dashed line: linear fitting of the solid line. The linear velocity dispersion is 6% in the 8–180-Hz band, which gives a Q estimate of 15. Note that frequency is in log scale.

b)

Figure 6. 共a兲 The 0.1–0.5 s of the correlated vibrator VSP section from the 3L-38 Mallik gas hydrate research well, vertical component. The source-borehole offset is 22 m. 共b兲 The raypath-average Pwave velocity obtained from well logging, compared with the geologic setting.

Figure 7. The T共f兲 of the direct wave in the vibrator sweep shown in Figure 3d, obtained by comparing the pilot and received t-f relations. The t-f relation of the pilot sweep, which is shown in Figure 3b, is calculated with the IPH method. The received t-f relation is measured with the MWCC method, with the pseudo t-f spectrum shown in Figure 5.

Figure 9. The T共f兲 profile for the direct wave in the vibrator sweeps received at different depth levels from the 3L-38 Mallik gas hydrate research well. The source-borehole offset is 22 m.

Figure 10. The profile of linear velocity dispersion for the Mallik data, calculated from the T共f兲 profile in Figure 9, compared with the geologic setting. The contours are raypath-average P-wave velocity. Note that frequency is in log scale.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WA118

Sun et al. samples will be studied to determine the relationship between the physical rock properties and observed attenuation and velocity dispersion. Various petrophysical models will be tested to interpret the observed attenuation and velocity dispersion in terms of fracture distribution, fluid fill, porosity, and other rock properties.

ACKNOWLEDGMENTS The Mallik 2002 Gas Hydrate Production Well Research Program was supported by the International Continental Scientific Drilling Program. Funding for the VSP data acquisition was provided by a University of Toronto start-up grant to Bernd Milkereit. We acknowledge the excellent field operations conducted by Schlumberger Ltd. for the offset VSP survey.

APPENDIX A Figure 11. The Q profiles for the Mallik data compared with the geologic setting. The QKK: calculated from the linear velocity dispersion shown in Figure 10. The QSR: calculated from the spectral ratio of first arrivals in the correlated vibrator data, which are shown in Figure 6. pared with the geologic setting. Although the velocity values are the raypath-average velocities, the relationship between the velocity dispersion and geologic setting can be observed. A raypath-average Q profile is calculated from linear velocity dispersion using equation 1, and is shown in Figure 11. This Q profile is compared with the Q profile obtained from the spectral ratio 共Tonn, 1991兲 of the first arrivals in the correlated vibrator traces. The two Q profiles are comparable and consistent with the results of Guerin et al. 共2005兲, Gauer et al. 共2005兲, and Pratt et al. 共2005兲. The Q values of a P-wave can be as small as 15 in gas hydrate-bearing sediments and are higher in gas hydrate-free sediments. Uncorrelated vibrator VSP data from other areas with different geologic settings have been analyzed 共Sun et al., 2007a, 2007b兲. Velocity dispersion as large as 10% in the 30–150-Hz band and extremely low Q have been observed in the data from McArthur River Mines, Athabasca Basin, Canada, where the metamorphic sediment is highly fractured and saturated with fluids. Velocity dispersion less than 1% in the 30–160-Hz band and moderate Q have been observed in the uncorrelated vibrator data acquired in fractured crystalline rocks in Outokumpu, Finland.

CONCLUSIONS A new method, the moving window crosscorrelation method, is presented for determining velocity dispersion with satisfactory accuracy and affordable cost in the exploration seismic frequency band, using uncorrelated vibrator VSP data. This determination of velocity dispersion provides an alternative method for measuring Q. Using the MWCC method, continuous velocity dispersion has been observed in the uncorrelated vibrator VSP data acquired in the Mallik gas hydrate research well. By fitting a straight line of the velocity dispersion data with respect to log frequency, a Q profile is obtained for the Mallik data. This Q profile is consistent with the one calculated using the spectral ratio method. Attenuation and velocity dispersion can provide better insight into physical rock properties. Therefore, future research will focus on the analysis of uncorrelated vibrator data from other regions with different geologic and petrophysical settings. Well logs and core

METHODS TO MEASURE VELOCITY DISPERSION IN UNCORRELATED VIBRATOR SWEEPS Cross-spectrum (CS) method Donald and Butt 共2004兲 measured velocity dispersion in their rock experiments from the cross-spectrum of the pilot and received signals

Ssr ⳱ Sr · Ss* , where Sr is Fourier transform of the received signal, and Ss* is the conjugate of Fourier transform of the pilot signal. The frequency-dependent traveltime T共f兲 then is

T共f兲 ⳱

1 d⌰ 共f兲 , 2␲ df

共A-1兲

where ⌰ 共f兲 is the unwrapped phase of Ssr. The CS method directly calculates the difference between the pilot and received t-f relations, other than individually. This method is applicable to seismic data of any source type. For example, the source signal in rock experiments by Donald and Butt 共2004兲 was an impulse-like wavelet, whereas in this case it consists of uncorrelated vibrator sweeps. This method relies on the assumption that the wave group propagates as a single mode and therefore is not applicable if multiple events exist in a seismic record, as shown in Figure A-1.

Instantaneous phase (IPH) method The instantaneous phase was introduced by Taner and Sheriff 共1977兲. The instantaneous frequency at a time t is

f共t兲 ⳱

1 d⌽ 共t兲 , 2␲ dt

共A-2兲

where ⌽ 共t兲 is the unwrapped phase of the Hilbert transform of an uncorrelated vibrator sweep. Similarly to the CS method, this method relies on the assumption that the wave group propagates as a single mode. This method measures the frequency at a certain time but not the arrival time as a function of frequency. To obtain a received t-f rela-

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

Velocity dispersion in the seismic band tion, f共t兲 must be converted to obtain the frequency-dependent arrival time. Thus, this method is convenient for measuring the t-f relation of a pilot sweep, but not for a received sweep. The derivatives in equations A-1 and A-2 can be calculated with difference quotient formulas.

Time-frequency spectral decomposition (TFSD) method The t-f relation of an uncorrelated vibrator sweep can be detected by picking the arrival time of an event at different frequencies in the t-f spectrum. A t-f spectrum is calculated using a time-frequency spectral decomposition algorithm. The clearest event in a received t-f spectrum usually is the direct wave. The T共f兲 is the difference between the pilot and received t-f relations. The TFSD method can be used to calculate the t-f relation of both pilot sweeps and received sweeps. The key procedure of this method is an appropriate TFSD to generate a t-f spectrum in which different events are separated satisfactorily. There are different approaches to conduct TFSD, as discussed by Castagna and Sun 共2006兲. The spectral decomposition algorithm used in this appendix is the discrete Fourier transform 共DFT兲.

Computing cost of the four methods The CS, IPH, TFSD, and MWCC methods are implemented in MATLAB 7.0 on a PC with a dual-core CPU and 2 GB memory. The computing time to obtain a t-f spectrum using the TFSD method, or a pseudo t-f spectrum using the MWCC method, increases with window length and decreases with the time step at which the window moves along a vibrator sweep:

Computing time ⬀

共Window length兲␣ , 共Time step兲␤

共A-3兲

where ␣ ⬇ 1, ␤ ⬇ 2 for the TFSD method, and ␣ ⬇ 0.7, ␤ ⬇ 0.8 for the MWCC method. Comparing the computing time of the four methods, and setting the CS method to be 1, the relative computing time of the other methods are as follows: the IPH method is 2; the TFSD method is 7 ⫻ 104 when the window length is 1000 sampling intervals and the time step is 1 sampling interval; and the MWCC method is 102 when the window length and time step are 2000 and 50 sampling intervals, respectively.

Measuring the t-f relation of a pilot sweep In a vibrator survey, the pilot sweep often is not exactly the same as predetermined because of site response, harmonics, and so forth. Hence the recorded pilot sweep should be used. The pilot t-f relation is obtained by measuring the instantaneous frequency of the pilot sweep using the IPH or TFSD method. Synthetic vibrator data were used in the tests discussed in this appendix. The sampling interval is 1 ms. The pilot sweep is linear, 12 s long, 8–180 Hz, and with 0.15 s taper, as shown in Figure A-1. Figure A-2 shows the errors in the pilot t-f relation measured using the IPH method without noise. The errors, small in the central part but large at both ends, are caused mainly by improper tapering of the pilot sweep. This is because the Hilbert transform is implemented using Fourier transform, in which a short taper results in a wavy spectrum. These errors can be reduced using a mean filter, as shown in Figure A-2.

WA119

Figure A-3 shows errors in the pilot t-f relation measured using the TFSD method without noise. The t-f spectrum is calculated with DFT using a 1-s cosine window. The errors in the central part are caused by the reduced frequency resolution caused by the short window. The errors always are greater than 0.1 Hz and are large at both ends in the time range marked by Tunst. The large errors at both ends result because a partial sweep being extracted from an end of a pilot sweep contains less effective data points. Magnitude of errors in the central part and the Tunst depend on the window length:

兩Error in f共t兲兩 ⱕ 0.5/Window length, and

Tunst ⬇ 0.5 ⫻ Window length. In Figure A-3, the window length is 1 s, so that the error amplitude is 0.5 Hz, and the Tunst is about 0.5 s. Therefore, increasing the window length reduces errors in the central part, but increases the Tunst. The TFSD method is remarkably more expensive but not necessarily more accurate than the IPH method. Therefore, the IPH method is preferable when determining pilot t-f relations. On the pilot sweep in Figure A-1, white noise is added to test the IPH method. The errors can be controlled by applying a mean filter on ⌽ 共t兲. Increasing filter length reduces errors in the central part but expands the unstable range at both ends.

a)

b)

c)

Figure A-1. The pilot sweep used in the synthetic tests, without noise. This sweep is linear, 12 s long, 8–180 Hz, sampling interval 1 ms. 共a兲 The first 2 s of the pilot sweep; 共b兲 the designed t-f relation of the pilot sweep; 共c兲 amplitude spectrum in the 0–200-Hz band of the pilot sweep.

Figure A-2. Errors in the measured f共t兲 of the pilot sweep 共Figure A-1兲 using the IPH method. No noise is added. Black line: no smoothing applied. Gray line: a 10-point mean filter has been applied on ⌽ 共t兲.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WA120

Sun et al.

If stationary noise exists at a frequency, the result of the IPH method is not stable in the time range in which the instantaneous frequency of the signal is close to the noise frequency. The large errors can be reduced by applying a median filter on ⌽ 共t兲. In theory, the IPH method cannot handle harmonics because of the single mode assumption. Fortunately, synthetic tests show that when the harmonics are weak, the result of the IPH method still is acceptable, as long as ⌽ 共t兲 is smoothed properly. Figure A-4 is an example. Amplitude of the harmonic is about 20% of the signal. Errors in the 0.74–11.92-s range are less than 0.1 Hz. However, if the harmonic is strong, e.g., 70% of the signal, the IPH method fails. A solution is to apply a pure phase-shift filter 共Li et al., 1995兲 on the pilot sweep before measuring f共t兲 using the IPH method. To conclude, the IPH method is efficient and applicable for pilot sweeps of ordinary data quality. The TFSD method is expensive and does not provide better results than the IPH method. Therefore, the IPH method is preferable.

Figure A-6 compares results from the CS, TFSD, and MWCC methods when the received sweep is free of noise. No smoothing is applied for the CS method. Errors in the CS method are negligible. A 1-s cosine window is used for the TFSD and MWCC method. Results of the TFSD and MWCC methods are not stable at both ends. For the TFSD method 共Figure A-6b兲, when there is no noise, the magnitude of the errors in the central part depends on the time step at which the window moves along the received sweep:

兩Error in T共f兲兩DFT ⱕ 0.5 ⫻ Time step. The errors are 0.5 ms in the central part because the window moves sample by sample, i.e., 1 ms. Therefore, to obtain T共f兲 measurements of satisfactory accuracy, the time step must be small enough, which makes the TFSD method extremely time-consuming 共equation A-3兲. For the MWCC method 共Figure A-6c兲, when there is no noise, magnitude error in the central part is determined by the sampling interval:

Measuring the frequency-dependent traveltime The CS method measures T共f兲 directly, whereas the other methods compare the pilot and received t-f relations to calculate T共f兲. Measuring the pilot t-f relation was discussed in the previous section. This section is a discussion of measuring the received t-f relation, i.e., frequency-dependent arrival time. In the tests below, the noise-free pilot sweep 共Figure A-1兲 is used to better assess the influence of noise in the received sweep. The received sweep, shown in Figure A-5, is calculated by assuming that the pilot sweep has propagated for 1 km in a medium with Q ⳱ 20 and linear velocity dispersion 共Figure 2b兲. The four methods are tested in the presence of noise and multiple events.

a)

b)

c)

Figure A-3. Errors in the measured f共t兲 of the pilot sweep, using the DFT method with a 1-s cosine window. 共a兲 The full length; 共b兲 zooming in for the central 1 s; 共c兲 zooming for the last 1 s, with the arrow and Tunst marking the time range in which the DFT method could not provide stable results.

兩Error in T共f兲兩MWCC ⱕ 0.5 ⫻ Sampling interval. The errors are within 0.5 ms in the central part because the sampling interval is 1 ms.Accuracy of the MWCC method does not depend on the time step at which the window moves along the pilot sweep; hence, the time step can be increased to save computing time 共equation A-3兲. This results in fewer f i readings. If white noise is added to the received sweep, none of the methods are stable when the signal-to-noise ratio 共S/N兲 is smaller than 1. When S/N⬎ 1, mean filters need to be applied to ⌰ 共f兲 for the CS method. For the MWCC method, a longer window length should be used to better tackle the noise. However, increasing the window length in the MWCC method smears the T共f兲. After proper smoothing, the results of the CS, TFSD, and MWCC methods are comparable. For stationary noise at a frequency f noise, results of the TFSD and MWCC methods are not affected significantly because stationary noise can be separated satisfactorily in a t-f spectrum or a pseudo t-f

a)

b)

c)

Figure A-4. Errors in the measured f共t兲 of the pilot sweep using the IPH method in the presence of a harmonic with an amplitude of 20% of the signal. Two 100-point mean filters were applied on the ⌽ 共t兲.

Figure A-5. The received sweep used in the synthetic tests without noise. This sweep was obtained by propagating the pilot sweep in Figure A-1 for 1 km in a medium with the Q and velocity dispersion shown in Figure 2a. 共a兲 The first 4 s of the received sweep. 共b兲 The T共f兲 calculated from the velocity model. 共c兲 Amplitude spectrum of the received sweep compared with that of the pilot sweep in the 8–180-Hz band.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

Velocity dispersion in the seismic band spectrum. The CS method failed around f noise, but the large errors can be eliminated using a medium filter. In field data, there always are multiple events in one received sweep, whereas the T共f兲 of the direct event is desired. The CS method cannot separate different events. Notches appear in the measured T共f兲, resembling the “ghosts” in marine seismic surveys as a result of interference. The TFSD and MWCC methods can provide satisfactory results if the events are well separated in time; otherwise, interference occurs. As an example, a secondary event is added with a time delay ⌬T to the received sweep in Figure A-2. As shown in Figure A-7, the CS method cannot separate multiple events even if the ⌬T is large and the amplitude of the second event is significantly smaller than the main event. Capability of the TFSD method in separating multiple events depends on the window length; a shorter window better separates different events because the events cannot be distinguished in an amplitude spectrum. Figure A-8 is an example. In Figure A-8, T共f兲 measurements of the main event are satisfactory when the window length is 0.8 s but have large errors when the window length is 2 s. For the MWCC method, the ability to separate multiple events increases with window length. Figure A-9 shows that the errors when the window length is 4 s are smaller than those when the window length is 3 s. However, a long window will smear the result. Conclusions of synthetic tests on determining T共f兲 of a received sweep are as follows: 1兲

The CS method is efficient but not suitable for field data because it cannot provide accurate results when multiple events exist.

WA121

a)

b)

Figure A-7. Measuring T共f兲 of a main event using the CS method when another event exists in the received sweep. 共a兲 The correlated seismogram of the received sweep. 共b兲 The measured T共f兲 compared with the model T共f兲.

a)

b)

c)

a)

Figure A-8. Measuring the T共f兲 of a main event in the presence of a second event using the TFSD method. 共a兲 The correlated seismogram of the received sweep. 共b兲 Errors in the T共f兲 measurements when a 0.8-s cosine window is used. 共c兲 Errors in the T共f兲 measurements when a 2-s cosine window is used.

b)

a)

c) b)

Figure A-6. Errors in the measured T共f兲 from 共a兲 the CS method, 共b兲 the DFT method, and 共c兲 the MWCC method. The pilot sweep is shown in Figure A-1, and the received sweep is shown in Figure A-5. No noise is added. No smoothing is applied for the CS method. A 1-s cosine window is used for the DFT and MWCC methods.

Figure A-9. Measuring the T共f兲 of a main event in the presence of a second event using the MWCC method. 共a兲 The correlated waveform of the received sweep. 共b兲 Errors in the measured T共f兲 using a 3 -s 共black line兲 and a 4-s 共gray line兲 cosine window.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WA122 2兲

3兲

4兲

Sun et al.

The IPH method is suitable for measuring the pilot t-f relation because it is efficient and provides satisfactory results for the pilot sweeps of ordinary data quality, e.g., when S/N of white noise is larger than 1 and harmonics are weak. The TFSD method is stable with white noise and stationary noise, but the capability of separating different events is not as satisfactory as that of the MWCC method. This method is the most expensive. The MWCC method is appropriate for measuring T共f兲 in a received vibrator sweep because it is stable with noise and can satisfactorily separate different events in the pseudo t-f spectrum. However, the window length must be selected carefully, as a longer window better tackles the noise and better separates different events, but more severely smears the variation in the T共f兲.

REFERENCES Adam, L., M. Batzle, and L. Brevik, 2006, Gassmann’s fluid substitution and shear modulus variability in carbonates at laboratory seismic and ultrasonic frequencies: Geophysics, 71, no 6, F173–F183. Batzle, M. L., D. H. Han, and R. Hofmann, 2006, Fluid mobility and frequency-dependent seismic velocity — Direct measurements: Geophysics, 71, no 1, N1–N9. Bourbié, T., O. Coussy, and B. Zinszner, 1987, Acoustics of porous media, translated by N. Marshall: Gulf Publication Co. Brajanovski, M., T. Müller, and B. Gurevich, 2006, Characteristic frequencies of seismic attenuation due to wave-induced fluid flow in fractured porous media: Geophysics Journal International, 166, 574–578. Brown, R. L., and D. Seifert, 1997, Velocity dispersion: A tool for characterizing reservoir rocks: Geophysics, 62, 477–486. Castagna, J. P., and S. Sun, 2006, Comparison of spectral decomposition methods: First Break, 24, April, 75–79. Dallimore, S. R., T. S. Collett, T. Uchida, and M. Weber, 2005, Overview of the science program for the Mallik 2002 Gas Hydrate Production Research Well Program, in S. R. Dallimore, and T. S. Collett, eds., Scientific results from Mallik 2002 Gas Hydrate Production Research Well Program, MacKenzie Delta, Northwest Territories, Canada: Geological Survey of Canada Bulletin 585. Donald, J. A., and S. D. Butt, 2004, Experimental technique for measuring phase velocities during triaxial compression tests: Rock Mechanics and Mining Sciences, 42, 307–314. Dragoset, W. H., 1988, Marine vibrators and the Doppler effect: Geophysics, 53, 1388–1398. Futterman, W. I., 1962, Dispersive body waves: Journal of Geophysical Research, 67, 5279–5291. Gauer, K., C. Haberland, R. G. Pratt, F. Hou, B. E. Mediole, and M. H. Weber, 2005, Ray-based cross-well tomography for P-wave velocity, anisotropy, and attenuation structure around the JAPEX/JNOC/GSC et al., Mallik 5L38 gas hydrate production research well, in S. R. Dallimore, and T. S. Collett, eds., Scientific results from Mallik 2002 Gas Hydrate Production Research Well Program, MacKenzie Delta, Northwest Territories, Canada: Geological Survey of Canada Bulletin 585. Guerin, G., D. Goldberg, and T. S. Collett, 2005, Sonic attenuation in the JAPEX/JNOC/ GSC et al., Mallik 5L-38 gas hydrate production research well, in S. R. Dallimore, and T. S. Collett, eds., Scientific results from Mallik 2002 Gas Hydrate Production Research Well Program, MacKenzie Delta, Northwest Territories, Canada: Geological Survey of Canada Bulletin 585. Gutenberg, B., 1958, Attenuation of seismic waves in the earth’s mantle:

Bulletin of the Seismological Society of America, 48, 269–282. Jannsen, D., J. Voss, and F. Theilen, 1985, Comparison of methods to determine Q in shallow marine sediments from vertical reflection seismograms: Geophysical Prospecting, 33, 479–497. Johnson, D. L., 2001, Theory of frequency dependent acoustics in patchysaturated porous media: Journal of the Acoustical Society of America, 110, 682–694. Jones, T. D., 1986, Pore fluids and frequency-dependent wave propagation in rocks: Geophysics, 51, 1939–1953. Kjartansson, E., 1979, Constant Q — Wave propagation and attenuation: Journal of Geophysical Research B, 84, 4737–4748. Knopoff, L., 1964, Q: Reviews of Geophysics, 2, 625–660. Knopoff, L., and G. J. F. MacDonald, 1958, Attenuation of small amplitude stress waves in solids: Review of Modern Physics, 30, 1178–1192. Lee, M. W., 2002, Biot-Gassmann theory for velocities of gas hydrate-bearing sediments: Geophysics, 67, 1711–1719. Li, X.-P., W. Söllner, and P. Hubral, 1995, Elimination of harmonic distortion in Vibroseis data: Geophysics, 60, 503–516. Liu, H.-P., D. L. anderson, and H. Kanamori, 1976, Velocity dispersion due to anelasticity; implications for seismology and mantle composition: Geophysical Journal of the Royal Astronomical Society, 47, 41–58. Mavko, G., T. Mukerji, and J. Dvorkin, 1998, The rock physics handbook: Cambridge University Press. Molyneux, J. B., and D. R. Schmitt, 1999, First break timing: Arrival onset times by direct correlation: Geophysics, 64, 1492–1501. Pratt, R. G., F. Hou, K. Gauer, and M. Weber, 2005, Waveform tomography images of velocity and inelastic attenuation from the Mallik 2002 crosshole seismic surveys, in S. R. Dallimore., and T. S. Collett, eds., Scientific results from Mallik 2002 Gas Hydrate Production Research Well Program, MacKenzie Delta, Northwest Territories, Canada: Geological Survey of Canada Bulletin 585. Pride, S. R., and J. G. Berryman, 2003, Linear dynamics of double-porosity dual-permeability materials: 1 — Governing equations and acoustic attenuation: Physical Review, E68, 036603. Sams, M. S., J. P. Neep, M. H. Worthington, and M. S. King, 1997, The measurement of velocity dispersion and frequency-dependent intrinsic attenuation in sedimentary rocks: Geophysics, 62, 1456–1464. Spencer, J. W., 1981, Stress relaxations at low frequencies in fluid-saturated rocks — Attenuation and modulus dispersion: Journal of Geophysical Research, 86, 1803–1812. Sun, L. F., B. Milkereit, and D. Schmitt, 2007a, Detecting heterogeneity near a borehole using Vibroseis VSP Data, in B. Milkereit, eds., Proceedings of Exploration 07: 5th Decennial International Conference on Mineral Exploration, 1091–1094. ——–, 2007b, Measuring attenuation and velocity dispersion using vibrator sweeps: 77th Annual International Meeting, SEG, Expanded Abstracts, 26, 3115–3119. Taner, M. T., and R. E. Sheriff, 1977, Application of amplitude, frequency, and other attributes to stratigraphic and hydrocarbon determination, in C. E. Payton, ed., Applications to hydrocarbon exploration: AAPG Memoir 26, 301–327. Tonn, R., 1991, The determination of the seismic quality factor Q from VSP data: A comparison of different computational methods: Geophysical Prospecting, 39, 1–27. White, J. E., 1975, Computed seismic speeds and attenuation in rocks with partial gas saturation: Geophysics, 40, 224–232. Winkler, K. W., 1986, Estimates of velocity dispersion between seismic and ultrasonic frequencies: Geophysics, 51, 183. Winters, W. J., S. R. Dallimore, T. S. Collett, B. E. Medioli, R. Matsumoto, T. J. Katsube, and P. Brennan-Alpert, 2005, Relationships of sediment physical properties from the JAPEX/ JNOC/GSC et al., 5L-38 gashydrate production research well, in S. R. Dallimore, and T. S. Collett, eds., Scientific results from Mallik 2002 Gas Hydrate Production Research Well Program, MacKenzie Delta, Northwest Territories, Canada: Geological Survey of Canada Bulletin 585. Zener, C. M., 1948, Elasticity and anelasticity of metals: University of Chicago Press.

Downloaded 19 Apr 2009 to 117.98.183.40. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/