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
SUMMARY
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.
INTRODUCTION
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.
COMPLEX RAY THEORY
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.
Formula
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.
Q ESTIMATION TECHNIQUES
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.
CONCLUSIONS
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.
ACKNOWLEDGEMENTS
Shell Global Solutions International BV provided financial sup-port to this project through CRDF grant RUG1-30019-ST-11.
http://dx.doi.org/10.1190/segam2013-0948.1
EDITED REFERENCES
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.
REFERENCES
