• Nie Znaleziono Wyników

Methods for source-independent Q estimation from microseismic and crosswell perforation shot data in a layered, isotropic viscoelastic medium


Academic year: 2021

Share "Methods for source-independent Q estimation from microseismic and crosswell perforation shot data in a layered, isotropic viscoelastic medium"


Pełen tekst


Methods for source-independent Q estimation from microseismic and crosswell perforation shot data

in a layered, isotropic viscoelastic medium

Renat Shigapov∗and Boris Kashtan, Saint Petersburg State University, Alexander Droujinine, Shell Global Solutions International BV, Wim Mulder, Shell Global Solutions International BV and Delft University of Technology


We consider the problem of anelastic full waveform inversion in a multi-layered, isotropic viscoelastic medium from micro-seismic and cross-hole perforation shots. Usually, the source wavelet is not well known, so we focus on Q estimation tech-niques that can handle this problem. We revisit the spectral ra-tio method (SRM) and discuss its domain of applicability. We conclude that SRM can be used for interval Q estimation if the sources and receivers are contained in the same homogeneous layer. To overcome the restrictions of SRM, we consider the more sophisticated source-independent full waveform inver-sion (SIFWI). We compare three different source-independent misfit functions on synthetic data and show that they can be used for Q estimation in a viscoelastic medium. We use 3-D complex ray theory to implement SIFWI in a layered vis-coelastic medium.


Accurate Q estimates are used for seismic image correction and can be related to the physical properties of the subsurface, such as saturation, porosity and permeability. Tonn (1991) made a comprehensive comparison of ten methods for Q esti-mation from vertical seismic profile (VSP) data and concluded that there is no single method that is better than the others. Our work is motivated by the problem of attenuation estima-tion from microseismic and perforaestima-tion shot data. The prob-lem has some characteristic features: 1) the source spectrum is difficult to estimate in practice, and 2) the data are usually so noisy that we can work with direct arrivals only.

The simplest method of attenuation estimation that handles the unknown source wavelet is SRM. Romero et al. (1997) applied SRM to wavefields excited by microearthquakes in order to get Q’s estimates. Using SRM, Eaton (2011) estimated Qp and

Qsvalues from real microseismic data in order to improve the

calculation of the microseismic source parameters. Shekar and Tsvankin (2012) applied SRM for the P-wave angle-dependent attenuation estimation from cross-hole data generated by per-foration shots in a horizontal borehole and recorded in a verti-cal borehole. For an arbitrary acquisition geometry, even in a layered medium SRM cannot be applied for interval Q estima-tion, at least without applying layer stripping algorithms. To overcome this restriction, FWI might be preferable.

For traveltime tomography without picking, Van Leeuwen and Mulder (2010) proposed a weighted correlation-based misfit function that reduces sensitivity to phase shifts in the incor-rect source wavelet. As a natural generalization, Luo and Sava (2011) proposed a deconvolution-based misfit function. For attenuation tomography, Plessix (2006) used a misfit function

based on the centroid frequency shift method (Quan and Har-ris, 1997), which requires a source wavelet. In contrast, SIFWI algorithms (Choi and Alkhalifah, 2011; Choi and Min, 2012) completely eliminate the source spectrum during the inver-sion process. These algorithms are driven by convolution- or deconvolution-based misfit functions using reference synthetic and observed data. Although Xu et al. (2006) have shown for some cases in an elastic medium that SIFWI yields results infe-rior to the results obtained by FWI with estimation of a source function, the capabilities of SIFWI for Q estimation are still unclear.

To implement FWI in a laterally homogeneous, multi-layered medium, we need a relatively accurate and fast algorithm for the computation of the direct arrival’s traveltime and ampli-tude. To this end, we focus on asymptotic ray theory. Hanyga and Seredyˇnska (1999) and Thomson (1997) used complex ray theory. ˇCerven`y (2005) and ˇCerven`y and Pˇsenˇcik (2009) de-scribed elastic ray tracing based on the perturbation approach. Vavry˘cuk (2008) proposed real viscoelastic ray tracing. Hearn and Krebes (1990a,b) showed that saddle point evaluation of a wavefield integral provides complex rays. Droujinine (2006) proposed generalized anelastic asymptotic ray theory.

Here, we have implemented complex ray tracing for a lay-ered viscoelastic medium. We consider Q estimation tech-niques that eliminate the source function. We begin with SRM and discuss its capabilities and limitations. We also consider SIFWI and analyze different misfit functions. We finish with an example of attenuation estimation in a layered medium.


The analytical solution for the wavefield excited by a point source in a layered viscoelastic medium can be constructed with the Fourier-Bessel transform. We replace the related Han-kel functions by there approximations for large arguments. Af-ter that, the inner integral can be evaluated by the saddle point method. This leads to a formula for the complex ray fields.


The k-th (k = 1, 2, 3) displacement component of a direct ar-rival, uk= uk(xs, xr,t), at the receiver point xrexcited by a

ver-tical point force at xswith a frequency-domain wavelet F3(ω)

in a 3-D layered viscoelastic medium can be described as

uk= F3∗ Gk3= 1 (2π)2 ρ1 Re ∞ Z 0 X n=P,S (−1)δnS γ n kγ3n− δnSδk3  v2n1γ3n q X ζnτn00 ·TnF3(ω) · e ω (τn−it)−i 2arg  ζnτn00  +iθ dω, (1)

where δi j and δnSare Kronecker symbols, with n = P for the

P-wave and n = S for the S-wave; τn= τn(xs, xr, ω) is the


complex traveltime function (n = P, S), defined as τn= iX ζn+ m P j=1 hj q

ζn2− v−2n j; m is the number of ray segments; i =

√ −1; X= s 2 P k=1

(xsk− xrk)2; ζn is the stationary n-wave complex

slowness defined as the solution of the equation τn0= 0, where

the prime denotes the derivative with respect to the complex slowness; vn jare the complex n-wave (n = P, S) velocities and

ρjis the real-valued density in the j-th layer, which is

asso-ciated with the j-th ray segment, with a thickness hj; the

in-terfaces of layers are assumed to be orthogonal to the x3-axis;

Re(·) takes the real part; Gk3stands for Green’s tensor; Tnis

the product of the anelastic single-interface transmission coef-ficients for nn-waves (n = P, S); γkn= vn1ζn·xkr for k = 1, 2;

γ3n= ivn1 q ζn2− v−2n1; arg  ζnτn00 

denotes the phase angle of ζnτn00; θ = −π4 for k = 1, 2, and θ =3π4 for k = 3. For the sign

of Fourier transform used here we define rootsqζn2− v−2n j = q ζn− v−1n j· q ζn+ v−1n j as follows: q ζn2− v−2n j = ζn− v −1 n j · ζn+ v −1 n j · e ϕ1+ϕ2 2 , where ϕ1= arg  ζn− v−1n j  ∈ [−3π 2,π2] and ϕ2= arg  ζn+ v−1n j  ∈ [−π 2,3π2].

The stationary slownesses ζP and ζSare the saddle points of

the P- and S-wave phase functions or, in other words, the cor-responding solutions of the following equations

iX+ m X j=1 ζnhj q ζn2− v−2n j = 0, n= P, S. (2) a) b)

Figure 1: (a) SRM for receivers, and (b) SRM for sources. The red star denotes a source. The green triangle stands for a receiver.


SRM for interval Q estimation from microseismic and cross-hole perforation shot data (SRM for receivers)

We abbreviate Unk(xsi, xr j, ω) with Unki j, where n = P or S and

k= 1, 2, 3. Note that in our notation, the real part of τn

cor-responds to attenuation. The logarithmic spectral ratio LSRnk

has the form

LSRnk= ln Unki j Unkil = ln Ξi, jlnk + ω Reτni j− τnil  , n= P, S, (3)

where the ratios Ξi, jlnk = ϒi, jlnkTnki, jl/Γi, jlnk of the directivity func-tions ϒi, jlnk = ϒi jnk/ϒilnk, transmission coefficients Tnki, jl= Tnki j/Tnkil and geometrical spreading Γi, jlnk = Γi jnk/Γil

nkat the j-th and l-th

receiver are assumed to be independent of frequency. Strictly speaking, this assumption is valid only for a homogeneous vis-coelastic medium, in which the rays are straight lines (Fig-ure 1a), which can be verified by substituting the analytical solution for the stationary slowness into equation 2. We are therefore able to estimate interval Q factors only by assuming that a couple of receivers and a microseismic event lie in the same homogeneous m-th layer. Using a robust least squares method, for example, with bisquare weights, we determine the slope Reτni j− τnil

in equation 3 on a LSRnk(ω)-plot and

ob-tain estimates of Qmn via

Qmn = π ri j− ril  vm nRe  τni j− τnil  , n= P, S, (4)

where ri jand rilare distances from the source i to the receivers

jand l correspondingly.

If Re 

τni j− τnil

is not a linear function of frequency in the source bandwidth, then the area where rays propagate is not homogeneous or/and Q depends on frequency. Figures 2a and 2b show true (SLS mechanism with Q = 11 at the peak fre-quency) and estimated Q values in a homogeneous viscoelastic medium calculated by SRM, assuming a) Q is constant, and b) Qdepends on frequency via the SLS mechanism. Thus, SRM can be used for frequency-dependent Q estimation if we know the attenuation mechanism and rays propagate in a homoge-neous region. Figure 1c shows the case of modeling with a constant Q = 11. The SRM estimates are QP= QS≈ 11.5,

be-cause we performed modeling using vn= vn(ω) (Figure 2d),

whereas for the estimate 4, we used a constant vn(ωp),

corre-sponding to the velocity at the peak frequency. Nevertheless, in the case of constant Q and small velocity dispersion (Fig-ure 2d), SRM provides quite accurate estimates of the Q’s. To obtain the true values of Q using SRM, we should use the true velocities vn= vn(ω) in equation 4.

SRM for interval Q estimation using microseismic events with the same source functions (SRM for sources) Here we propose a method for Q estimation from microseismic events with identical source functions. Consider two events, i and l, and one receiver, j. Assuming that both the sources and the receiver lie in the same laterally homogeneous layer (Figure 1b), m, and considering spectral ratio of wavefields recorded at j from two events, i and l, we have the estimate

Qmn = π ri j− rl j  vm nRe  τni j− τnl j  , n= P, S, (5)

which is the same as equation 4 with source and receiver posi-tions swapped.

Based on the same idea, we could separate sets of events with similar source functions if the nonlinearity of Reτnl j− τnli

 in the source bandwidth is induced by different source func-tions rather than by the nonlinear behavior of Q and the


a) b)

c) d)

Figure 2: Panels (a) and (b) show the true Q’s and the SRM es-timates when the modeling was performed for a standard linear solid with Q = 11 at the peak frequency. (a) SRM estimates, assuming LSRnk(ω) is a linear function with two unknowns, so

Qis frequency-independent. The estimated Qp= Qs= 13.6.

(b) SRM estimates, assuming LSRnk(ω) is a nonlinear

func-tion of three unknowns. Panel (c) displays the true and esti-mated Qp= Qs≈ 11.5, using SRM, when the modeling was

performed for a constant Q = 11. Panel (d) shows P- and S-wave velocities, density, and source S-wavelet used in each ex-ample.

mogeneity of a layer. Separation of events into such subsets might be helpful in the analysis of microseismic data.

Source-independent full waveform inversion

As opposed to inversion for Q with a well-known wavelet (Ri-bodetti and Virieux, 1998; El Yadari et al., 2008), SIFWI elim-inates the source function but requires reference traces. The choice of reference traces can strongly affect the quality of the attenuation estimates. It is therefore important to clarify the SIFWI capabilities for Q estimation on a simple example. We consider three source-independent (SI) misfit functions and a source-dependent (SD) one. Let Di j(ω) be the Fourier

trans-form of the observed data di j(t) and Ui j(ω, Q) of the modeled data ui j(t, Q), where Q is the vector with N

Qelements of Qp

and Qsvalues in each layer. A typical SD misfit function is just

the well-known `2-norm of the residuals between the modeled

and observed data,

ESD= Ns X i=1 Nr X j=1 3 X k=1 1 2 U i j nk− D i j nk 2 `2, n= P, S, P ± S. (6)

A convolution-based misfit function in the frequency domain, proposed by Choi et al. (2005), is

E1SI= Ns X i=1 Nr X j=1 3 X k=1 1 2 D i nk,re fU i j nk−U i nk,re fD i j nk 2 `2 , (7)

where Dink,re f is the k-th component of the observed refer-ence wavefield averaged over all the Nrreceivers from the

i-th source, and Unk,re fi stands for the corresponding modeled reference wavefield.

A deconvolution-based misfit function, proposed by Zhou and Greenhalgh (2003) and Lee and Kim (2003), looks like

E2SI= Ns X i=1 Nr X j=1 3 X k=1 1 2 Unki j Unk,re fi ! − D i j nk Dink,re f ! 2 `2 . (8)

Finally, the logarithmic SI misfit function (Choi and Min, 2012) can be viewed as a natural generalization of SRM:

E3SI= Ns X i=1 Nr X j=1 3 X k=1 1 2 ln U i j nk Unk,re fi ! − ln D i j nk Dink,re f ! 2 `2 . (9)

Unfortunately, the last two misfit functions in any case require additional regularization of the wavefields because of the divi-sions and the domain of definition of the logarithm. We first analyzed the behavior of these functions in a homogeneous viscoelastic medium. We assume that all parameters except Qpand Qsare known. Hak and Mulder (2010, 2011) showed

that causality and the acquisition geometry have a strong in-fluence on the sensitivity function in a viscoacoustic medium. Our numerical tests confirm that this is also the case for a vis-coelastic medium. We therefore assume from now on that the frequency-dependent term in the velocity definition is known. The acquisition geometry should at least satisfy the condition that each receiver does not lie in a shadow zone of a source. Here, we use only four widely spaced receivers (Figure 5b). We have also found that for E1SI the best choice for a

refer-ence trace is not just averaging over all the receivers, but we should also average over all three component. Then, the ob-tained reference wavefield is a scalar field for a given source.

Figure 3 shows the normalized values of the various misfit functions as a function of Qpand Qsin a homogeneous

vis-coelastic medium. We use the same a 50-Hz Ricker wavelet for both ‘real’ and ‘synthetic’ wavefields. In the case of the same source wavelets for ‘real’ and ‘synthetic’ data and with the assumption that we know all the medium parameters except the Q’s, all misfit functions have their minumum in the correct place. The behavior of deconvolution-based misfit function is somewhat inferior. Figure 4 shows the normalized values of the various misfit functions with different source wavelets in ‘real’ and ‘synthetic’ wavefields. We used a 50-Hz Ricker wavelet for ‘synthetic’ and a 60-Hz derivative of a Gaussian for the ‘real’ data. Figure 3a shows that the SD misfit function provides incorrect values for the attenuation, whereas the SI ones have the correct minimum. We also see that for the cho-sen parameters, the convolution-based misfit function has the best behavior. Therefore, SIFWI has a potential for Q estima-tion under idealized circumstances.


a) b)

c) d)

Figure 3: Normalized misfit functions: (a) ESD, (b) E1SI, (c)

E2SI, (d) E3SIwith the true values Qp= Qs= 60.1. We used

the same 50-Hz Ricker wavelet in both ‘real’ and ‘synthetic’ simulations. For (a) and (b), we did not apply any regulariza-tion. In case (c), we regularized the reference fields by adding 10−3of its maximum absolute value. In case (d), we regular-ized also the ‘real’ and ‘synthetic’ wavefields to avoid prob-lems with the logarithm.

Iterative minimization in a layered medium

We consider iterative SIFWI based on E1SI for Q estimation

in a layered medium. We minimized E1SI by three different

iterative methods: a) steepest descent, b) conjugate-gradient, and c) full Newton. We calculated the gradient using Fr´echet derivatives. In cases b) and c), the minimization procedure often breaks down (negative Q values arise) if no constraint on the gradient is applied. We therefore clip the gradient of the misfit function in order to eliminate unrealistic negative Q values. Figure 5 shows results obtained by a steepest-descent method in a three-layered medium using only P-waves. We used a 50-Hz Ricker wavelet for the synthetic and a 50-Hz derivative of a Gaussian for simulating the real data. The small difference between the estimated and true Q values is caused by the constant step length for each iteration. The estimation of the optimal step length can reduce this difference.


We have considered the problem of Q estimation from syn-thetic microseismic and cross-hole perforation shot data. We revisited the capabilities of SRM and proposed simple scenar-ios of Q estimation from microseismic data. We also checked

a) b)

c) d)

Figure 4: Normalized misfit functions: (a) ESD, (b) E1SI, (c)

E2SI, (d) E3SI. We used a 50-Hz Ricker wavelet for the

‘syn-thetic’ and a 60-Hz derivative of a Gaussian for simulating the ‘real’ data. All the conditions except for the ‘real’ source wavelet are the same as in Figure 3.

Qtrue = 35 Qtrue = 55 Qtrue = 60 Qini = 50 Qini = 70 Qini = 45 Qest = 36 Qest = 54 Qest = 60 a) b)

Figure 5: (a) Three-layered true, initial and estimated models. (b) Acquisition geometry.

the capabilities of SIFWI for attenuation estimation. We im-plemented a complex ray tracing algorithm for fast calcula-tion of the direct arrivals in a layered viscoelastic medium. With this forward-modeling algorithm, we have shown that convolution-based SIFWI with a poorly known source wavelet can be used for interval Q estimation in a laterally homoge-neous multi-layered medium.


Shell Global Solutions International BV provided financial sup-port to this project through CRDF grant RUG1-30019-ST-11.




Note: This reference list is a copy-edited version of the reference list submitted by the author. Reference lists for the 2013 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.


Aki, K., and P. Richards, 2002, Quantitative seismology: University Science Books.

Cerveny, V., 2005, Seismic ray theory: Cambridge University Press.

Cerveny, V., and I. Psencik, 2009, Perturbation Hamiltonians in heterogeneous anisotropic weakly

dissipative media : Geophysical Journal International, 178, 939–949,


Choi, Y., and T. Alkhalifah, 2011, Source-independent time-domain waveform inversion using convolved

wavefields: Application to the encoded multisource waveform inversion: Geophysics, 76, no. 5,

R125–R134, http://dx.doi.org/10.1190/geo2010-0210.1.

Choi, Y., and D.-J. Min, 2012, Source-independent elastic waveform inversion using a logarithmic

wavefield : Journal of Applied Geophysics, 76, 13–22,


Choi, Y., C. Shin, D.-J. Min, and T. Ha, 2005, Efficient calculation of the steepest descent direction for

source-independent seismic waveform inversion: An amplitude approach: Journal of Computational

Physics, 208, 455–468, http://dx.doi.org/10.1016/j.jcp.2004.09.019.

Droujinine, A., 2006, Generalized anelastic asymptotic ray theory: Wave Motion, 43, 357–367,


Eaton, D. W., 2011, Q determination, corner frequency and spectral characteristics of microseismicity

induced by hydraulic fracturing: 81st Annual International Meeting, SEG, Expanded Abstracts,

1555–1559, http://dx.doi.org/10.1190/1.3627499.

El Yadari, N., F. Ernst, and W. Mulder, 2008, Near-surface attenuation estimation using

wave-propagation modeling: Geophysics, 73, no. 6, U27–U37, http://dx.doi.org/10.1190/1.2976548.

Hak, B., and W. A. Mulder, 2010, Migration for velocity and attenuation perturbations : Geophysical

Prospecting, 58, 939–952, http://dx.doi.org/10.1111/j.1365-2478.2010.00866.x.

Hak, B., and W. A. Mulder, 2011, Seismic attenuation imaging with causality: Geophysical Journal

International, 184, 439–451, http://dx.doi.org/10.1111/j.1365-246X.2010.04848.x.

Hanyga, A., and M. Seredynska, 1999, Asymptotic ray theory in poro- and viscoelastic media : Wave

Motion, 30, 175–195, http://dx.doi.org/10.1016/S0165-2125(98)00053-5.

Hearn, D., and E. Krebes, 1990a, Complex rays applied to wave propagation in a viscoelastic medium:

Pure and Applied Geophysics, 132, no. 1-2, 401–415, http://dx.doi.org/10.1007/BF00874371.

Hearn, D., and E. Krebes, 1990b, On computing ray-synthetic seismograms for anelastic media using

complex rays: Geophysics, 55, 422–432, http://dx.doi.org/10.1190/1.1442851.

Lee, K., and H. Kim, 2003, Source-independent full-waveform inversion of seismic data: Geophysics, 68,

2010–2015, http://dx.doi.org/10.1190/1.1635054.


Luo, S., and P. Sava, 2011, A deconvolution-based objective function for wave-equation inversion: 81st

Annual International Meeting, SEG, Expanded Abstracts, 2788–2792,


Plessix, R.-E., 2006, Estimation of velocity and attenuation coefficient maps from crosswell seismic data:

Geophysics, 71, no. 6, S235–S240, http://dx.doi.org/10.1190/1.2345051.

Quan, Y., and J. M. Harris, 1997, Seismic attenuation tomography using the frequency shift method:

Geophysics, 62, 895–905, http://dx.doi.org/10.1190/1.1444197.

Ribodetti, A., and J. Virieux, 1998, Asymptotic theory for imaging the attenuation factor Q: Geophysics,

63, 1767–1778, http://dx.doi.org/10.1190/1.1444471.

Romero, A., T. McEvilly , and E. Majer, 1997, 3-D microearthquake attenuation tomography at the

Northwest Geysers geothermal region, California : Geophysics, 62, 149–167,


Shekar, B., and I. Tsvankin, 2012, Anisotropic attenuation analysis of crosshole data generated during

hydraulic fracturing: The Leading Edge, 31, 588–593, http://dx.doi.org/10.1190/tle31050588.1.

Thomson, C., 1997, Complex rays and wave packets for decaying signals in inhomogeneous, anisotropic

and anelastic media : Studia Geophysica et Geodaetica, 41, no. 4, 345–381.

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,


Van Leeuwen, T., and W. A. Mulder, 2010, A correlation-based mifit criterion for wave-equation

traveltime tomography: Geophysical Journal International, 182, 1383–1394,


Vavrycuk, V., 2008, Real ray tracing in anisotropic viscoelastic media: Geophysical Journal International,

175, 617–626, http://dx.doi.org/10.1111/j.1365-246X.2008.03898.x.

Xu, K., S. Greenhalgh, and M. Wang, 2006, Comparison of source-independent methods of elastic

waveform inversion: Geophysics, 71, no. 6, R91–R100, http://dx.doi.org/10.1190/1.2356256.

Zhou, B., and S. Greenhalgh, 2003, Crosshole seismic inversion with normalized full-waveform

amplitude data: Geophysics, 68, 1320–1330, http://dx.doi.org/10.1190/1.1598125.


Powiązane dokumenty

jest przedm iotem sprzedaży, działu, dziedziczenia, a także obciążenia hipoteką, służebnościam i

Miłoszowe Teraz nieustannie odnosi się do Zawsze, a pojawiający się w polu świadomości byt oglądany jest przez odwrotną stronę gobelinu. Kiedy więc powiada Miłosz: „od

W analizę owych dodatkowych pieśni (od trzeciej zwrotki począwszy) nie wdajemy się narazie: zlepek ich z »Bogurodzicą« jest bowiem całkiem mechaniczny, tylko

Warto zaznaczyć, że kategoria ariergardy jest przez Perloff stosowana w sposób anachroniczny: aby opisać pewien nurt, który nie mieści się ani w ramach awangardy, ani w

W części czwartej Winkler zanalizował wydarzenia roku 1923 stanowiące wyraźną cezurę w dziejach ruchu ro ­ botniczego i niejako zamykające zapoczątkowany pod

Wydaje się również, że we w nikliw ym i obszernym wstępie zabrakło miejsca dla dwóch kwestii: zasygnalizowano jedynie zmasowaną nagonkę na Kota jako historyka

De vraag is: hoe brengen we de ondergrond terug in de stedenbouwkundige plannen?Als de afstemming van de verschillende belangen in stedelijke ontwikkeling wordt meegenomen in

Z obszaru polskiej k u ltu ry znane są jego badania twórczości Jan a Kochanowskiego (w perspektyw ie kom paratystycznej) i Ignacego Krasickiego (studia nad Panem