1. technology and computing

Quality factors estimation using wavelet’s envelope peak instantaneous frequency
Jinghuai Gao*, Senlin Yang, Inst. Wave & Information, Xi'an Jiaotong University, Xi'an, China
Daxing Wang, Research Inst. of E & D, Chang-Qing Oil-Field Company of CNPC, Xi’an, China
Rushan Wu, Modeling and Imaging Laboratory, IGPP, University of California, Santa Cruz, California, USA
Summary
We derive an approximate equation combining the quality
factor Q, the traveltime of a wave, and the variation of the
instantaneous frequency (IF) at the envelope peaks of two
successive seismic wavelets, along the wave propagating
direction, based on the theory of one-way wave
propagation in a 1D viscoelastic medium. We then propose
a method (called the WEPIF method) to estimate Q by
measuring the variations of the wavelet envelope peak IF
(WEPIF) with the traveltime of seismic wavelet. A test on
synthetic VSP data shows that the WEPIF method is less
sensitive to interference from the reflector than the
logarithm spectral ratio and the centroid frequency shift
methods. Applied to field VSP data, the WEPIF method
gives a Q-curve with nearly the same distribution as the
results from a known well. Applied to poststack seismic
data, it produces a Q-profile that indicates an intense
absorption zone corresponding to the excellent gas-bearing
reservoir. This allows us to predict a potential highproductivity gas well. Drilling confirmed this prediction.
The WEPIF method can be applied to poststack seismic
data and zero-offset VSP data.
Introduction
Seismic waves propagating through the earth suffer
attenuation and dispersion due to the viscosity of the media
(Ricker, 1953; Futterman, 1962). Q is a diagnostic tool for
reservoir characterization and hydrocarbon detection
(Toksöz et al., 1979). Moreover, Q is very important when
interpreting the effects of AVO, improving the resolution
of seismic imaging, and advancing the study of material
properties.
In general, amplitude variations of seismic signals are used
to estimate Q values (Tonn, 1991). Time-domain methods
involve pulse amplitude decay (Brzostowski and
McMechan, 1992), pulse rising time (Gladwin and Stacey,
1974; Kjartansson, 1979), and pulse broadening (Wright
and Hoy, 1981), which all use pulse amplitude information.
Nevertheless, amplitude information of seismic pulses is
often influenced by scattering, geometric spreading, and
other factors. In the Fourier-frequency domain, approaches
for estimating Q include the logarithm spectral ratio (LSR)
(Hauge, 1981; Stainsby and Worthington, 1985), centroid
frequency shift (CFS) (Quan and Harris, 1997), and peak
frequency shift (Zhang and Ulrych, 2002; Gao and Yang,
2007) methods, all of which require Fourier transforms to
calculate the frequency spectra of seismic records sampled
within a time window. In practice, it is difficult to properly
select the window function and window length. Li et al.
(2006) suggested using peak scale variations in the wavelet
domain to estimate Q by assuming an idealized pulse as the
seismic source wavelet. Its application is, however,
restricted due to the fact that the difference between the real
source wavelet and an idealized pulse may be substantial.
Q can also be estimated by the variation of the
instantaneous frequency (IF) of a seismic signal. Tonn
(1991), Barnes (1991), and Engelhard (1996) obtained the
relationship between the measured instantaneous spectra
and seismic attenuation. Assuming that the source wavelet
is an idealized band-pass wavelet, Barnes (1991) derived
the relation between Q and IF variations of seismic waves,
establishing a new approach for Q estimation.
In this work, we develop Barnes’s work (Barnes, 1991),
and propose a method for estimating Q. we derive an
approximate equation combining Q, the traveltime of the
wave, and the variation of the IF at the envelope peaks of
two successive seismic wavelets, along the wave
propagating direction, based on the theory of one-way
wave propagation in a 1D viscoelastic medium. We then
propose a method (called the WEPIF method) to estimate Q
by measuring the variations of the wavelet envelope peak
IF (EPIF) with its traveltime. We finally test our method
using both synthetic and field seismic data.
Theories of four Q estimation methods
Considering a one-way plane wave propagating in a
horizontally layered anelastic medium with a frequencyindependent Q, we can determine the wavefield of a source
wavelet traveling through a distance z by (Aki and
Richards, 1980; Stainsby and Worthington, 1985)
⎡ iω z
ωz ⎤,
−
Uˆ ( z , ω ) = Uˆ ( 0 , ω ) e x p ⎢ −
⎥
2 c (ω ) Q ⎦
⎣ c (ω )
(1)
where i= -1 , ω is the angular frequency, z and c(ω )
denote the travel-distance and phase velocity, respectively,
Q is the quality factor of the medium, and Uˆ (0, ω ) is the
wavefield of the source signature. For VSP and reflection
seismograms, we assume that the source wavelet can be
approximated by
SEG Houston 2009 International Exposition and Annual Meeting
2457
Quality factors estimation using wavelet’s envelope peak instantaneous frequency
⎛δ2 ⎞
u (0, t ) = A ⎜
⎟
⎝ π ⎠
1/ 4
⎡
δ 2t 2 ⎤
exp ⎢ i(σ t + ϕ ) −
⎥,
2 ⎦
⎣
(2)
where σ is the modulated frequency (i.e., dominant
frequency), A and ϕ are the amplitude and phase factor,
respectively, and δ is the energy decay factor of the
wavelet. The wavelet, defined by equation 2, has four
parameters ( A,ϕ ,σ , δ ) to be determined properly so that it
can match the actual source signature. Applying the Fourier
transform to equation 2, the source wavefield can be
expressed as
Uˆ (0, ω ) = A
( )
1/ 4
4π
2
δ
2
⎡ (ω − σ )
⎤
exp ⎢ −
+ iϕ ⎥ .
2δ 2
⎣
⎦
(3)
According to the definition of Barnes (1991), equation 3 is
a constant-phase wavelet. A constant-phase wavelet, being
non-causal, is not physically realistic; however, a real
causal source wavelet can be converted to a constant-phase
wavelet by an appropriate phase rotation (Toksöz and
Johnston, 1981; Barnes, 1991). Also, wavelets derived by
seismic-to-well correspondence often have a near-constant
phase (Longbottom et al., 1988). Neglecting velocity
dispersion (i.e., c(ω ) = c ) and substituting equation 3 into
equation 1, we achieve
1/ 4
⎛ 4π ⎞
Uˆ ( z , ω ) = A ⎜ 2 ⎟
⎝δ ⎠
⎡ (ω − σ )2 ω z
ω z ⎞⎤
⎛
−
+ i ⎜ϕ −
exp ⎢ −
⎟⎥ .
c ⎠ ⎥⎦
2δ 2
2Qc ⎝
⎢⎣
(4)
In order to find the relation between Q and the EPIF, we
first discuss the EPIF of a seismic wavelet. For a constantphase wavelet propagated for time τ through a
homogeneous anelastic medium, its EPIF is exactly equal
to the average Fourier frequency weighted by its amplitude
spectrum (Sheriff, 1984; Robertson and Nogami, 1984;
Bodine, 1986; Barnes, 1991), i.e.,
∞
f p (τ ) = ∫ fA(τ , f )df
0
∫
∞
0
A(τ , f )df ,
the amplitude spectrum and the EPIF of the wavelet after
traveltime τ , respectively, A(τ , f ) =| Uˆ ( z , f ) | , and
τ = z / c . From equations 3 and 5, we get
⎡
exp ⎢ −
⎣
∞
⎡ 2π 2
exp ⎢ −
2
0
⎣ δ
δ
4π
∫
2
2
2
⎛ σ ⎞ ⎤
⎜
⎟ ⎥
⎝ 2π ⎠ ⎦
,
2
σ ⎞ ⎤
⎛
⎜ f −
⎟ ⎥ df
2π ⎠ ⎦
⎝
2π
δ
2
2
σ
τδ
f p (τ ) =
−
+
2π 4π Q
2
(6)
⎡ 2π 2 ⎛ σ
δ2
τδ 2 ⎞
exp ⎢ − 2 ⎜
−
⎟
2
4π
⎢⎣ δ ⎝ 2π 4π Q ⎠
2
⎤
⎥
⎥⎦ . (7)
2
2
⎡ 2π 2 ⎛
⎤
∞
⎞
σ
τδ
∫ 0 exp ⎢⎢ − δ 2 ⎜⎝ f − 2π + 4π Q ⎟⎠ ⎥⎥ df
⎣
⎦
After analyzing characteristics of the seismic wavelets
which are extracted from zero-offset VSP data and
poststack seismic data, using the first-order Taylor
approximation for equations 6 and 7, we obtain an
approximated relationship as follows (Gao et al., 2008)
Q ≈ τδ 2κ (η ) /(4 π Δ f p )
(8)
where σ and δ are the modulated angular-frequency and
the standard deviation of the wavelet, Δf p is the EPIF
variation between the reference pulse and received pulse,.,
2 2
η = σ (2πδ )-1 , κ (η ) = 1 − 2πη e −2π η φ −1 (2πη ) , and φ ( x)
is
the
probability
integration
of
x
standard
normal
distribution, i.e., φ ( x ) = (2π) −0.5 ∫ e −0.5t dt . We call the
2
−∞
method of Q estimation by equation (8) as the WEPIF
method for short.
Estimation of Q by equation 8 requires the parameters such
as EPIF, traveltime and parameters of the reference pulse to
be known.
For a given depth z , wavelet u ( z , t ) is a function of time
determined by parameters ( A, δ ,σ ,ϕ ) . Letting
ψ (t , A,σ , δ ,ϕ ) be the real part of u ( z, t ) and w(t ) be the
reference pulse, the wavelet parameters can be solve a
optimal problem
(5)
where f = ω 2π is the frequency, A(τ , f ) and f p (τ ) are
σ
f p (0) =
+
2π
where f p (0) is the EPIF of the source wavelet. Similarly,
by inserting equation 4 into equation 5, we obtain
min
( A ,δ ,σ ,ϕ )
∫ | u ( z, t ) − w(t ) | dt .
(9)
For seismic data with a high signal-to-noise ratio (SNR),
the IF is usually computed using the derivative of the
instantaneous phase. However, field seismic data is usually
contaminated by random noise. Large spikes in the IF occur
when the denominator of equation 12 approaches zero more
rapidly than the numerator suggested that the IF should be
replaced by the damped IF when the large spikes are
uninterested for our purpose. For stabilization in IF, the
damped IF or weighted IF may be employed (Matheney
and Nowack 1995). For noisy signals, IFs can also be
extracted in the wavelet domain (Gao et al., 1999).
SEG Houston 2009 International Exposition and Annual Meeting
2458
Quality factors estimation using wavelet’s envelope peak instantaneous frequency
a)
b)
a)
LSR
source
density: 2.0 g/cc
c=1.0 km/s
Q=30
3420
3440
3460
3480
T rue Q
LSR
c)
Depth (m)
280
300
320
340
360
480
500
560
580
600
-50 -25
1200
3480
-10
0
10
20
30
40
Q
3380
3400
3400
3420
3420
3440
3440
3460
3460
0
50
100
GR (API)
520
540
geophones
1100
3380
3480
420
440
460
layer 3
600
Depth (m)
Depth (m)
240
260
380
400
density: 3.0 g/cc
c=2.0 km/s
Q=100
1000
WEPIF
200
220
400
900
CFS
140
160
180
layer 2
3440
3460
100
120
density: 2.5 g/cc
c=1.25 km/s
Q=60
3420
Time (ms)
60
80
200
T arget
reservoirs
3400
800
0
20
40
layer 1
CFS
WEPIF
3400
b)
0
3380
3380
Depth (m)
We examine the validity of our method using synthetic
zero-offset VSP data. Figure 1a shows a 3-layer depth
model. Based on this model, we calculate the synthetic
seismogram using the method proposed by Ganley (1981).
The source signature is a 50 Hz constant-phase wavelet
with δ = 75 located at the surface with zero-offset. There
are 51 geophones; the interval between two successive
geophones between 0 and 400 m is 10 m, while between
400 and 600 m, it is 20 m. The sampling interval in time is
2 ms. For the sake of comparison, we use the LSR, CFS,
and WEPIF methods to estimate Q values as shown in
Figure 1b. Each of these methods has a few depth ranges in
which it can not work. The ranges are on the order of 160210 and 350-420 m for the LSR method, 150-210 and 320420 m for the CFS method, and 170-210 and 360-420 m for
the WEPIF method. These methods all fail near interfaces
due to serious overlapping of reflected waves, such as at
170-210 and 360-420 m, but the WEPIF method works
well where direct waves are partially overlapped by
reflection waves, such as at 150-170 and 320-360 m.
Apparently, the Q values estimated by the WEPIF method
are better than those of the LSR and CFS methods. This test
shows that the WEPIF method is less sensitive to
interference reflections and has a higher resolution.
and the distance between two adjacent geophones is 20 m.
The interpretation log shows that the well is drilled through
23.7-m gas sandstone at the target depth of 3395-3470 m.
Figure 5b is the Q curve estimated by the LSR, CFS, and
WEPIF methods. The depth ranges of target reservoirs are
indicated by blue crosses in Figure 2b. Figure 2c shows the
natural gamma ray (GR) and acoustic (AC) logs. In general,
gas-bearing sandstone has a lower Q and more intensive
absorption than other rocks. The lower Q values estimated
by the LSR method are at depths of 3420-3480 m, and that
estimated by the CFS method are at depths of 3400-3460 m.
The intensive absorption zone estimated by our method is
at the depth ranges of 3380-3400 and 3420-3440 m, which
agrees primarily with the distribution of the gas reservoirs.
This test shows that the intensive absorption zone relates to
the productive gas reservoir.
Depth (m)
Results
0
25
50
75 100 125 150
Q
Figure 1 The VSP model (a) and (b) Q-factors estimated by LSR
method, CFS method, and WEPIF method.
We then test our method using one-shot zero-offset VSP
records of 6-levels as shown in Figure 2a in which the
depth range is 3380-3480 m, the sample interval is 1 ms,
150
200
3480
170
230
290
350
AC (us/m)
Figure 2 (a) One-shot zero-offset VSP data of six levels, (b) Qfactors estimated by the LSR (black solid line), CFS (green solid
line), and WEPIF methods (red line: IF calculated with Hilbert
transform method). The GR and AC logs are given in (c). The
depth ranges of target reservoirs are indicated by blue crosses in
(b).
We also use reflection seismic data to test the validity of
our method. Figure 6a is a poststack seismic profile through
wells in a gas field in which the sample rate is 1 ms. We get
a Q −1 -profile as shown in Figure 3b by the WEPIF method.
For visualization, we also smooth the result in Figure 3b
with a 2-D Gaussian function in the f-k domain and show it
SEG Houston 2009 International Exposition and Annual Meeting
2459
Quality factors estimation using wavelet’s envelope peak instantaneous frequency
in Figure 3c. Two wells shown in Figure 3b-6c, well 1 (the
solid line in white) and well 2 (the dotted line in black),
were drilled in this survey. Well 1 is drilled through 16.4 m
of gas-bearing sandstone at the target reservoir (Shan 23
layer) and obtained 12.9815 m3/d industrial gas streams, yet
well 2 is drilled through only 8 m of sandstone without a
developed reservoir and proved to be a dry well. The GR
logs and AC logs of wells 1 and 2 are given in Figures 3d
and 3e, respectively. At the target depth (around 1.3 s)
shown in Figures 3b-3c, well 1 corresponds to an intensive
absorption zone, whereas well 2 corresponds to a weak
absorption zone. The result also shows a relation between
the intensive absorption zone and the well-developed gas
reservoir. Based on this observation, we predict a gas well
shown as well 3 denoted by a dashed line in purple.
Fortunately, well 3 is drilled through 9.2 m of gas-bearing
sandstone and obtained an open flow potential of 4.1481
m3/d at the gas test site. The GR logs and AC logs of well 3
are shown in Figure 3f. Although some unavoidable
unphysical negative Q values appeared due to factors such
as reflection/transmission coefficients and scattering of thin
beds, this test shows the validity of our prediction of a gas
reservoir based on the Q estimation.
Conclusions
We propose the WEPIF method for estimating Q. A test
using synthetic VSP data shows that the WEPIF method is
more stable, more convenient when selecting parameters,
and less sensitive to interface reflections than the LSR and
CFS methods. The results of zero-offset VSP data and
poststack seismic data show that this method is valid for
seismic attenuation estimation and gas reservoir
characterization.
Acknowledgements
This work was supported by the NSFC (No. 40730424, and No.
40674074), and the HTRDPC (No. 2006AA09A102).
1200
1225
1225
1250
1250
1250
1275
1275
1275
1275
1300
1300
1300
1300
1325
1325
1325
1325
1350
1350
1350
1375
10
1375
140
1375
10
1200
1225
90
170
GR (API)
250
220
300
AC (us/m)
380
90
170
GR (API)
250
f)
1200
1200
1225
1225
1225
1250
1250
1250
1275
1275
1300
1300
1325
1325
1350
1350
1350
1375
140
1375
10
Time (ms)
e)
1200
1200
Time (ms)
Time (ms)
d)
220
300
AC (us/m)
380
90
170
GR (API)
250
1375
140
220
300
380
AC (us/m)
Figure 3 Q-value estimations of reflection seismic data by WEPIF method, (a) poststack seismic data of a certain gas-field, (b) the Q -1
profile, (c) smoothed Q -1 profile, and (d)-(f) are the GR logs and AC logs for wells 1-3 respectively. In (b) and (c), wells 1-3 are denoted by
solid line in white, dotted line in black, and dashed line in purple, respectively.
SEG Houston 2009 International Exposition and Annual Meeting
2460
EDITED REFERENCES
Note: This reference list is a copy-edited version of the reference list submitted by the author. Reference lists for the 2009
SEG Technical Program Expanded Abstracts have been copy edited so that references provided with the online metadata for
each paper will achieve a high degree of linking to cited sources that appear on the Web.
REFERENCES
Aki, K., and P. G. Richards, 1980, Quantitative seismology: Theory and methods: W. H. Freeman & Co.
Barnes, A. E., 1991, Instantaneous frequency and amplitude at the envelope peak of a constant-phase wavelet: Geophysics, 56,
1058–1060.
Boashash, B., 1992, Estimating and interpreting the instantaneous frequency of a signal, part I: Fundamentals, part II: Algorithms
and applications: Proceedings of the IEEE, 80, 520–568.
Bodine, J. H., 1984, Waveform analysis with seismic attributes: 54th Annual International meeting, SEG, Expanded Abstracts,
505–509.
Brzostowski, M., and G. McMechan, 1992, 3-D tomographic imaging of near-surface seismic velocity and attenuation:
Geophysics, 57, 396–403.
Engelhard, L., 1996, Determination of seismic-wave attenuation by complex trace analysis: Geophysical Journal International,
125, 608–622.
Futterman, W. I., 1962, Dispersive body waves: Journal of Geophysical Research, 67, 5279–5291.
Gabor, A., 1946, Theory of communication: Journal of the Institution of Electrical Engineers, 93, 429–457.
Gao, J. H., and S. L. Yang, 2007, On the method of quality factors estimation from zero-offset VSP data: Chinese Journal of
Geophysics, 50, 1026–1040.
Gao, J. H., X. L. Dong, W. B. Wang, Y. M. Li., and C. H. Pan, 1999, Instantaneous parameters extraction via wavelet transform:
IEEE Transactions on Geoscience and Remote Sensing, 37, 867–870.
Gao J. H., S. L. Yang, D. X. Wang, and R. S. Wu, 2008, Estimation of quality factor Q from the instantaneous frequency at the
envelope peak of a seismic signal: submitted to Geophysics for publication.
Ganley, D. C., 1981, A method for calculating synthetic seismograms which include the effects of absorption and dispersion:
Geophysics, 46, 1100–1107.
Hauge, P. S., 1981, Measurements of attenuation from vertical seismic profiles: Geophysics, 46, 1548–1558.
Kjartansson, E., 1979, Constant Q-wave propagation and attenuation: Journal of Geophysical Research, 84, 4737–4748.
Li, H. B., W. Z. Zhao, H. Cao, F. C. Yao, and L. Y. Shao, 2006, Measures of scale based on the wavelet scalogram with
applications to seismic attenuation: Geophysics, 71, no. 5, V111–V118.
Matheney, M. P., and R. L. Nowack, 1995, Seismic attenuation values obtained from instantaneous frequency matching and
spectral ratios: Geophysical Journal International, 123, 1–15.
Quan, Y., and J. M. Harris, 1997, Seismic attenuation tomography using the frequency shift method: Geophysics, 62, 895–905.
Ricker, N., 1953, The form and laws of propagation of seismic wavelets: Geophysics, 18, 10–40.
Robertson, J. D., and H. H. Nogami, 1984, Complex trace analysis of thin beds: Geophysics, 49, 344–352.
Sheriff, R. E., 1984, Encyclopedia dictionary of exploration geophysics: SEG, 194.
Sheriff, R. E., and L. P. Geldart, 1995, Exploration seismology: Cambridge University Press.
Taner, M. T., F. Koehler, and R. E. Sheriff, 1979, Complex seismic trace analysis: Geophysics, 44, 1041–1063.
Toksöz, M. N., and D. H. Johnston, 1981, Seismic waves attenuation: SEG.
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.
Waters, K. H., 1978, Reflection seismology: John Wiely.
Wright, C., and D. Hoy, 1981, A note on pulse broadening and anelastic attenuation in near-surface rocks: Physics of the Earth
and Planetary Interiors, 25, P1–P8.
Zhang, C. J., and T. J. Ulrych, 2002, Estimation of quality factors from CMP records: Geophysics, 67, 1542–1547.
SEG Houston 2009 International Exposition and Annual Meeting
2461