Original paper

 

QVOA techniques for fracture characterization

 

Vladimir Sabinin

 

Instituto Mexicano del Petróleo, México D.F., México. Corresponding author: vsabinin@yahoo.com

 

Received: September 30, 2011.
Accepted: April 02, 2013.
Published on line: September 30, 2013.

 

Resumen

Se desarrollan nuevas técnicas de cálculo en el análisis del factor de calidad frente a la distancia fuente receptor y el acimut (QVOA) para caracterización de fracturas. Estas técnicas se aplican a datos sintéticos superficiales de reflexión con ruido.

Palabras clave: QVOA, factor de calidad, medio HTI, anisotropía sísmica, caracterización de yacimientos fracturados.

 

Abstract

New computational techniques of QVOA analysis (Quality factor Versus Offset and Azimuth) for fracture characterization are developed. The techniques are applied to synthetic surface data of reflection with noise.

Key words: QVOA, quality factor, HTI medium, seismic anisotropy, fracture-reservoir characterization.

 

Introduction

Quality factor Q is a seismic characteristic of attenuation property of a medium. By the value of Q-factor, one may estimate what kind of liquid fills pores of the medium (a reservoir characterization), because the attenuation of reservoir is due to presence of a liquid in the pores. For estimating factor Q, many different computational methods have been suggested (e.g., Tonn, 1991; Quan & Harris, 1997; Dasios et al., 2001; Zhang & Ulrych, 2002), as well as methodologies of its application to real data (Dasgupta & Clark, 1998).

QVOA analysis (Q-factor Versus Offset, and Azimuth) was first intended by Chichinina et al. (2004, 2005). Then it was further developed by Chichinina et al. (2006a, 2006b, 2006c), and Zhu and Tsvankin (2006). In fractured anisotropic media, anisotropy of Q is present (Clark et al., 2001; Chichinina et al., 2004, 2005). Namely, the value of quality factor Q depends on offset and azimuth. An approximate formula of this dependence was first suggested by Chichinina et al. (2006a, 2006b) for transversely isotropic media with a horizontal axis of symmetry (HTI), as well as an idea to apply it to estimation of principal fracture directions. Formulation of the QVOA analysis by Chichinina et al. (2006b) is similar to formulation of AVOA analysis (Amplitude Versus Offset, and Azimuth) by Rüger (1998) which gives good results in fracture characterization for media without attenuation (Sabinin & Chichinina, 2008).

Here computational techniques of QVOA analysis for estimation of fracture directions are suggested. Some of them are developed by analogy with known AVOA techniques (Mallik et al., 1998; Jenner, 2002). The others are original.

The QVOA techniques deal with estimating an attenuation factor (or quality factor Q) in surface data of reflection, what is not a trivial problem. Methodologies of such estimation of Q are proposed by Behura & Tsvankin (2009), Shekar & Tsvankin (2011), Reine et al. (2012), and Sabinin (2012).

 

Estimation of a symmetry axis angle by QVOA

Attenuation of wave amplitudes due to dispersion obeys a law (e.g., Tonn, 1991):

where r is a travel distance, and β is an absorption coefficient.

Quality factor Q is defined (e.g., Sheriff, 2002) as the ratio of 2π times the peak energy to the energy dissipated in a cycle:

It relates to β as follows (Futterman, 1962):

where , ƒ is frequency, and τ is travel time. For large values of Q, the approximation is valid.

Let us define the attenuation factor α as

Value α varies between 0 and 1.

For a model of viscoelastic media with complex stiffness tensor, Carcione (1995) defined the formula where V is the complex velocity, which was used by Chichinina et al. (2006a, 2006b). This definition is valid only for large Q values. In the general case, following the definition of Q (Equation 2), it must be

For anisotropic media with attenuation, from (5), by following Chichinina et al. (2006b), one can describe a dependence of the attenuation factor on a wave-ray travel (incidence) angle θ, and on a source-receiver azimuth Φ by the approximate formula:

where (Chichinina et al., 2006c)

A, B₀, C₀, are constants, provided B₀ > 0 for HTI media (Chichinina et al., 2004, 2006c), and Φ₀ is the symmetry axis angle, usually unknown, which has to be estimated.

The incidence angle θ has a sense of the angle between axis z and the wave ray in the anisotropic medium with attenuation. The approximation was made under the assumption of a small value of sin⁶θ, which practically means: sin²θ<0.36, or θ <37º.

 

Numerical methods for estimating the symmetry axis angle

One can compute the value α in the target anisotropic layer with attenuation from seismograms of P-waves reflected from the top and the bottom of this layer (see next section). Usually, the data used are 3D seismic data from spaced receivers and sources within nodes of a rectangle grid at the surface. The symmetry axis angle is usually obtained for a square surrounding a node of the grid (for a bin), by using seismic traces which Common Middle Point (CMP) is within this bin. If such traces are few, then neighbor bins are combined into a superbin, and calculations are made for it. Therefore, a preliminary stage of the estimation is a selection and an extraction of seismic traces of the superbin from the 3D data.

All methods described below are applied to the individual superbin. Let the superbin have n traces (i = 1,..., n) with incidence angles θ_i in the target layer and azimuthal angles Φ_i.

The first two methods use also sectoring; i.e., the traces of the superbin are sorted into m azimuthally equal sectors (j = 1,..., m), with k_j traces in each sector, with incidence angles θ_k (k = 1,..., k_j) in the target layer. It is adopted that all traces in the individual sector j have the same value of azimuth Φ_j, equal to the middle azimuth of the sector.

 

Approximate method of sectors (ASM)

It is the most computationally simple method based on the idea by Mallik et al. (1998) for the AVOA. From (6) one can write for the sector j:

where α_jk is the value of α calculated from the trace k in the sector j.

Having α_jk and θ_k for all k in the sector j (k_j ≥ 3), one can calculate A_j, B_j, and C_j from (8) by the least-squares method (e.g., Sabinin & Chichinina, 2008). These calculations have to be made for all sectors.

The unknown value Φ₀ may be obtained from the first of equations (7) written for each sector j:

As a rule, B_j is calculated with errors incorporated by the least-squares method from the data. Equation (9) is valid only for precise, theoretical values of B_j, and reasonably can be replaced by another. Let B_j = B_j^(p) δ_j where B_j^(p) is the precise value, and δ_j is an error. Then equation (9) can be replaced by the following approximate equation, which is proved to be more convenient for calculations:

If the error δ_j is not large, then a and b must be close to each other for an HTI medium (a ≈ b), although generally saying a ≠ b.

System (10) has three unknowns (a, b, and Φ₀), therefore m should be at least 3 for obtaining a solution. It is the same system as in the AVOA technique by Sabinin & Chichinina (2008). System (10) is solved by the least-squares method, and has two solutions (two Φ₀ differing by π/2 with the same b but of opposite sign). The condition b > 0 may be used to distinguish symmetry-axis and fracture-strike directions (see Equation 7).

 

Method of sectors (SM)

It differs from ASM by the method of solving the azimuthal problem (9). Instead of the approximate equation (10), equation (9) is written in its precise form:

where

This simplification leads to a more complicated solution. Let us consider the functional of error

Functional F must be minimized over parameters b and Φ₀. For this, it is necessary to solve the system of two equations:

Transformation of (13) gives the following non-linear equation for obtaining Φ₀:

Additionally, b = C / A

Equation (14) is trigonometric non-linear and is solved by the method of half-dividing. It usually has more than one solution for Φ₀. From these solutions, one chooses the one that gives the minimum for the functional (12). This choosing can be erroneous for highly noised data.

 

Approximate truncate method (ATM)

This method is similar to Jenner (2002) for the AVOA. The sectoring is not needed. Equation (6) is truncated for simplicity, and it is combined with (10), as in ASM. It gives after transformation:

where α_i is the value α calculated from the trace i of the superbin.

If we define s_i = sin² θ_i, S = sin(2Φ₀), C = cos(2Φ₀), g_i = cos(2Φ₁), and h_i = sin(2Φ_i), then equation (15) can be rewritten into a more convenient form as:

The values s_i, g_i, and hi are known because they can be calculated from the headers of the seismograms and the parameters of the medium.

Let us consider the functional of error

Functional F must be minimized over parameters a, b, d, and Φ₀. For this, it is necessary to solve the system of four equations:

The solution of system (17) gives the equation for obtaining the parameter Φ₀:

The other parameters are:

From (18), one can see that the solution Φ₀ has a period of . This value of the period means that we must use an additional condition for understanding which value of Φ₀ is the symmetry-axis azimuth. Application of the method has shown that the condition b > 0, as by analogy with ASM, fails sometimes.

 

Truncate method (TM)

This method differs from ATM by using equation (11) instead of (10). It leads to the following equation for superbin instead of (15):

where ξ_i = cos² (Φ₁ - Φ₀).

In order to solve (19) by the least-squares method, let us consider the functional of error

Functional F must be minimized over parameters a, b, and Φ₀. For this purpose, it is necessary to solve the system of three equations:

Transformation of (21) gives the non-linear equation for obtaining Φ₀:

Additionally, b = B₁ / C₁, and a = A₁ / C₁

Equation (22) is trigonometric non-linear and it is solved by the method of half-dividing. It usually has more than one solution for Φ₀. From these solutions, one chooses the one that yields the minimum for the functional (20).

 

General method (GM)

This method differs from TM by using the full form of equation (6). One can write instead of (19):

Let us consider the functional of error

Functional F must be minimized over parameters a, b, c, and Φ₀. For this, it is necessary to solve the system of four equations:

Three first equations of the system (25) give expressions for the parameters a, b, and c:

The fourth equation of (25) can be transformed into a non-linear equation for obtaining Φ₀:

The system (26) – (27) is solved by the method of half-dividing. Similar to SM and TM, the system usually has more than one solution for Φ₀. One must choose the one that gives the minimum for the functional (24). The condition b > 0 is used to choose the symmetry-axis direction.

 

Estimation of the attenuation factor

A correct estimation of the attenuation factor α from seismograms is very important for the proposed QVOA techniques. There are some good methods of estimating α (or quality factor Q), see Introduction. Here, it is suggested a methodology adapted to surface data of reflection.

The Frequency Shift Method (FS) by Quan & Harris (1997) was chosen here because it operates with integral values and hence is less sensitive to noise and gives values of α close to the classical Spectral Ratio Method (e.g., Tonn, 1991).

The ratio of spectral amplitudes of waves reflected from the bottom and from the top of the target layer with attenuation (see Figure 1) can be expressed from equations (1) and (3) as (e.g., Jannsen et al., 1985)

where A_t and A_b are amplitudes of spectra of reflected waves from top and bottom of the target layer, respectively, R(ƒ) is a coefficient that combines reflectivity set and geometrical spreading, and τ is the travel time of the ray inside the target layer.

With the assumption of weak dependency of R with frequency ƒ, equation (28) can be rewritten in the form:

where η = dτ.

Following Quan & Harris (1997), the coefficient h of equation (29) is calculated by the algorithm:

where sums are taken over an interval of frequencies where η is nearly constant (see Figure 1).

Value a is calculated via equations (3) and (4), where d = η/τ.

Data with noise can introduce significant error in estimating α values. For more ability of QVOA techniques, only values with , where is the arithmetic mean value of a over traces, are taken into consideration by the techniques.

For applying the FS method, it is necessary to choose a window for the impulse in the time domain, and a window (interval) in the spectral domain. If the seismic source generates a Ricker wavelet, I would choose the time window including three central phases of the impulse. For real data, if it is difficult to determine the three phases of the impulse because of noise, then one can chose an impulse of less number of phases. Then I would compose a wave with much more samples than the impulse truncated by this window, by adding zeros after the impulse. Typically for seismic problems, the window has near 300 samples, and the composed wave has 16384 samples. Then I would make the Fast Fourier transform of this wave to obtain the spectrum. The more samples the better spectrum.

Also, as it was noticed by Dasgupta & Clark (1998) for real data whose amplitudes (and spectra) are deformed by the normal moveout (NMO) correction, it is necessary to restore the data to its pre-NMO form.

The spectra of the impulse usually have the shape of a bell. The spectral interval should be chosen there where a plot of the logarithm of spectral ratio 1nS (see Figure 1) has a part with linear behavior (η ≈ const), and near the peak of the spectra in order to decrease errors in calculations. By observing spectra of impulses from many synthetic seismograms, it was found that the spectral interval between the peak frequency of the top spectrum (t) and 0.8 of the peak frequency of the bottom spectrum (b) is quite satisfactory. For making the position of interval more precise, one can use an adaptive procedure of the best fitting to the linear part of 1nS with the least-square method.

An important value in the calculation of α is the travel time τ of the ray inside the target layer. It is greater than the difference in time between reflected impulses at the same trace, Δt.

A system of non-linear equations can be derived by the tracing method for calculating τ in a multilayered reservoir (Sabinin, 2012). It allows to obtain the difference between Δt and τ, and also the incidence angle θ, and the refraction point x shown in Figure 2.

Alternative methods were developed by Behura & Tsvankin (2009), Shekar & Tsvankin (2011), Reine et al. (2012).

The difference between Δt and τ can be seen in Figure 3. It increases significantly with incidence angle or offset.

As can be seen in Figure 2, one should use the wave reflected from the point x at the top of the target layer to calculate the spectral amplitude At (Reine et al., 2012; Sabinin, 2012). This reflected wave can arrive at the surface between receivers. Knowing the value of x, one may obtain this wave (or its spectrum) by interpolation from waves (or spectra) of adjacent traces, taking into account different values of geometrical spreading.

 

Example of aplication of the techniques

The techniques are compared in ability to give the most precise value of symmetry axis angle Φ₀ for HTI medium. At present, reliable field methods of obtaining Φ₀ do not exist. Therefore, I generated synthetic seismograms by applying the technique by Sabinin (2012) of 2D viscoelastic modeling. A three-layer area was constructed with an anisotropic viscoelastic layer in the middle. I set Φ₀ = 60º, and derived models for Φ_j = 0º, 30º, 45º, 75º, 90º by rotating the stiffness tensor for the anisotropic layer.

The stiffness tensor for HTI medium (e.g., Chichinina et al., 2006b) was rotated by 60º to obtain the model for Φ_j = 0º, by 30º for Φ_j = 30º, 90º and by 15º for Φ_j = 45º, 75º. Anisotropic parameters Δ_N = 0.35, and Δ_T = 0.2 were used in the stiffness tensor.

The attenuation was assumed isotropic in the anisotropic layer, and values of relaxation times were chosen to obtain the factor Q near 60. Host rock velocities V_P in three layers from surface were given the values 3200, 4000, and 4800 m/s, V_S were twice less than V_P ; densities were the same for the three layers, and thicknesses of the first two layers from surface were 1590 and 410 m.

A source of explosion type generated one Ricker impulse of frequency 45 Hz. Receivers were spaced over every 200 m beginning from the position of the source, and measured the z-component of velocity.

For data being quasi-real, I added a random Gauss normal (white) noise to the synthetic seismograms generated. Maximum amplitude of the noise was chosen as 10% of the maximum amplitude of the impulse reflected from the top boundary of the target middle layer in the first trace.

In Figure 4, the seismogram for azimuth 0º is presented. As one may see, the amplitude of noise reached up to 50% of the wave amplitude in the middle traces, and up to 100% in far traces.

In Table 1, the symmetry axis angles calculated with the proposed numerical methods are compared. The methods were applied both, to traces smoothed by splines of third order and to non-smoothed (natural) traces.

For noisy data, the result of estimating Φ₀ depends on the choice of time windows. The results in Table 1 are obtained with an automatic choice of time windows by a hyperbolic dependence on offset with the initial window set manually. This manual setting was varied up to 10% of the time window.

The results presented in Table 1 are for the best choice of some smart attempts of defining the time windows at the traces. The maximum difference between the estimated Φ₀ and the correct value of 60º in these attempts, characterizes a sensitivity of the method to the window choice. The approximate values of the sensitivity for smoothed traces are also presented in Table 1.

 

Discussion and conclusion

The results in Table 1 show superior quality of the General method (GM). It gives stable non-sensitive values with small error. This may be explained by the synthetic nature of the data.

The other methods give good results too, but with a large sensitivity to the choice of time window. For example, 20% for ATM and SM means that one may obtain the value of Φ₀ with error up to ± 12º (0.2x60º).

Smoothing impulses seems to be useful for reducing sensitivity. Smoothing spectra also improves the results.

In synthetic data without noise, all techniques give nearly precise results with low sensitivity.

One can strongly conclude from results for GM that the third term of the approximation in Equation (6) is important when considering synthetic data.

In applying to real field data, the methods can give other results because of other structure of noise.

In real data, the role of the third (squared) term of approximation (6) can become small or wrong. Therefore, the precision of non-truncated methods may be reduced.

The least-squares method used in the techniques gives the better solution, the operating with more traces. Therefore, the methods of sectors (ASM, SM) should give worse results than the others because they operate with much less traces: coarsely k_j ≈ n/m.

It should also be noted that SM and ASM are constructed for sectored data and must really on having an additional error because of the need to set a unique value of azimuth for all traces in a sector. This error must be less by decreasing the width of the sector. It must also depend on disposition of the traces inside the sectors.

The maximum value of this error can be easily estimated. For example, let the width of the sector be equal to 10º, and all traces in the sector in the example of previous section be disposed not at the middle azimuth, as above, but at the border of the sector; i.e. at Φ_j = 5º, 35º, 50º, 80º, 95º. This case will consequently correspond to Φ₀=65º, and to an error of 5º in Φ₀.

All suggested methods are not ideal and can be non-reliable in real data. The methods GM, TM, and SM demand the solution of non-linear equations and further choosing one solution from many. This choosing may fail in real data. Contrary, the methods ASM and ATM give only one pair of solutions but use approximate equation (10) instead of (9).

The criterion B₀ > 0 for distinguishing the pair of symmetry axis and fracture strike directions, see (7), may unexpectedly fail in real data, too.

The logarithm of spectral ratio may not have the straight linear part near the peaks of the spectra, even in smoothed real data, and therefore the proposed algorithm of setting spectral windows may lead to wrong values of α.

The smoothing by splines is also not the best solution for improving real data.

In such conditions of problems with real data, the best strategy in fracture-reservoir characterization is to apply all these methods together and check if the results obtained with different methods are close to each other. The more methods, the better forecasting.

Here, we discussed some common problems of QVOA techniques. A more careful investigation of the techniques in synthetic data needs to be done in more complicated experiments. The investigation of different methods of smoothing is also interesting. This is a deal of future publications.

Numerical experiments with real data must not use Φ₀, but other parameters for a comparison of the techniques, because of the correct value Φ₀ is not known in field data. Such work is beyond the scope of the present publication, and might be hold yet.

Here, the numerical techniques for estimation of fracture directions by applying theory of QVOA analysis were suggested. One of them (GM) is proved good in synthetic data with noise.

 

Acknowledgement

This work was partially supported by SENER and CONACYT of the Mexican government in the framework of project Y.00108 "Atributos sismicos azimutales de atenuacion y amplitud en datos multicomponentes en el mapeo de fracturas".

 

Bibliography

Behura J., Tsvankin I., 2009, Estimation of interval anisotropic attenuation from reflection data, Geophysics, 74, 6, A69-A74.         [ Links ]

Carcione J.M., 1995, Constitutive model and wave equations for linear, viscoelastic, anisotropic media, Geophysics, 60, 2, 537-548.         [ Links ]

Chichinina T., Sabinin V., Ronquillo-Jarillo G., 2004, P-wave attenuation anisotropy in fracture characterization: numerical modeling for reflection data: 74th Annual International Meeting, SEG, Expanded Abstracts, 143–146.         [ Links ]

Chichinina T., Sabinin V. Ronquillo-Jarillo G., 2005, QVOA analysis as an instrument for fracture characterization: SEG 75th Annual International Meeting, paper ANI2.2, p.127-130.         [ Links ]

Chichinina T.I., Sabinin V.I., Ronquillo-Jarillo G., Obolentseva I.R., 2006a, Аzimuthal QVO analysis applied to oil and gas exploration, Russian Geology and Geophysics, 47, 2, p. 265-283.         [ Links ]

Chichinina T., Sabinin V., Ronquillo-Jarillo G., 2006b, QVOA analysis: P-wave attenuation anisotropy for fracture characterization, Geophysics, 71, 3, C37-C48.         [ Links ]

Chichinina T., Sabinin V., Ronquillo-Jarillo G., 2006c, QVOA method or Q-anisotropy for fracture-reservoir characterization. EAGE 68th Conference & Exhibition – Vienna, Austria, 12-15 June 2006. P147.         [ Links ]

Clark R.A., Carter A.J., Nevill P.C., Benson P.M., 2001, Attenuation measurements from surface seismic data: azimuthal variation and time–lapse case studies: 63rd Meeting, European Association Exploration Geophysicists, Extended Abstracts, paper L28.         [ Links ]

Dasgupta R., Clark R.A., 1998, Estimation of Q from surface seismic reflection data. Geophysics, 63, 2120-2128.         [ Links ]

Dasios A., Astin T.R., McCann C., 2001, Compressional-wave Q estimation from full-waveform sonic data, Geophysical Prospecting, 49, 353-373.         [ Links ]

Futterman W.I., 1962, Dispersive body waves: J. Geophys. Res., 67, 5279-5291.         [ Links ]

Jannsen D., Voss J., Theilen F., 1985, Comparison of methods to determine Q in shallow marine sediments from vertical reflection seismograms, Geophysical Prospecting, 33, 479-497.         [ Links ]

Jenner E., 2002, Azimuthal AVO: Methodology and data examples, The Leading Edge, 21, 8, 782-786.         [ Links ]

Mallik S., Craft K.L., Meister L.J., Chambers R.E., 1998, Determination of the principal directions of azimuthal anisotropy from P-wave seismic data: Geophysics, 63, 692-706.         [ Links ]

Quan Y., Harris J.M., 1997, Seismic attenuation tomography using the frequency shift method, Geophysics, 62, 3, 895-905.         [ Links ]

Reine C., Clark R., van der Baan M., 2012, Robust prestack Q-determination using surface seismic data: Part 1 – Method and synthetic examples, Geophysics, 77, 1, R45-R56.         [ Links ]

Rüger A., 1998, Variation of P-wave reflectivity with offset and azimuth in anisotropic media: Geophysics, 63, 935-947.         [ Links ]

Sabinin V., 2012, Viscoelastic modeling and factor Q for reflection data. Geofisica Internacional, 51, 4, 377-391.         [ Links ]

Sabinin V., Chichinina T., 2008, AVOA technique for fracture characterization: resolving ambiguity, Geofisica Internacional, 47, 1, p.3-11.         [ Links ]

Shekar B., Tsvankin I., 2011, Estimation of shear-wave interval attenuation from mode-converted data, Geophysics, 76, 6, D11-D19.         [ Links ]

Sheriff R.E., 2002, Encyclopedic dictionary of exploration geophysics, SEG, 4^th ed., 429 p.         [ Links ]

Tonn R., 1991, The determination of the seismic quality factor Q from VSP data: a comparison of different computational methods, Geophysical Prospecting, 39, 1-37.         [ Links ]

Zhang C., Ulrych T.J., 2002, Estimation of quality factors from CMP records, Geophysics, 67, 5, 1542-1547.         [ Links ]

Zhu Y., Tsvankin I., 2006, Plane-wave propagation in attenuative transversely isotropic media: Geophysics, 71, 2, T17 - T30.         [ Links ]