Print version ISSN 0016-7169
Geofís. Intl vol.47 no.1 México Jan./Mar. 2008
AVOA technique for fracture characterization: resolving ambiguity
V. Sabinin and T. Chichinina1
1 Instituto Mexicano del Petróleo. Eje Central Lázaro Cárdenas, 152, Mexico, City, Mexico. Email: email@example.com
Received: July 14, 2006
Accepted: October 25, 2007
Se desarrolla una nueva técnica de cálculo en el análisis de amplitud contra distancia fuentereceptor y contra acimut (AVOA) que resuelve la ambigüedad en la estimación de las direcciones principales de fracturas. Se analiza el tercer término de la aproximación de Rüger para el coeficiente de reflexión de ondas P, para distinguir la dirección de fracturas y la dirección de simetría. La técnica se prueba con datos reales del campo y datos sintéticos.
Palabras clave: AVOA, el medio HTI, anisotropía sísmica, caracterización de yacimientos fracturados.
A new computational technique of AVOA analysis (Amplitude Versus Offset and Azimuth) is developed which resolves an ambiguity in fracturedirection estimation. It analyzes a third term of Rüger's approximation for Pwave reflection coefficient to distinguish symmetryaxis and fracturestrike directions. The technique is tested on synthetic and field data.
Key words: AVOA, HTI media, seismic anisotropy, fracturereservoir characterization.
The analysis of azimuthal variation in reflection coefficients or AVOA analysis (Amplitude Versus Offset and Azimuth) is widely applied for detecting and mapping highly fractured zones with azimuthallyoriented vertical cracks. Based on the Rüger (1998) approximation for the reflection coefficients in HTI medium, the AVOA technique gives, in general, two principal symmetry directions of HTI medium without any formal indication of which of these two directions points to the symmetry axis, and which to the fracture strike. As shown by Zheng et al. (2004), some additional information is required besides the reflection coefficients to resolve the ambiguity in fracturedirection estimation.
We develop a computational technique for PPwave AVOA analysis improved in comparison with existing techniques, which resolves the ambiguity using only information obtained from amplitudes (or reflection coefficients). The technique analyzes the third term of Rüger's approximation for determining symmetryaxis and fracturestrike azimuths uniquely. We test the algorithm on synthetic and field data.
The methodology of AVOA analysis is based on the concept of azimuthal anisotropy caused for the most part by parallel vertical fractures. It leads to the azimuthal anisotropy of amplitudes, in particular, to azimuthal variation in reflection coefficients. The fractured reservoir is represented by a model of a transversely isotropic medium with horizontal symmetry axis (HTI medium). The FFwave reflection coefficient R at the interface (or reflecting boundary) between weakly anisotropic HTI media (or between isotropic and HTI media) is defined by the approximate formula (Rüger, 1998):
where θ is the incidence angle, and is the sourcereceiverline azimuth with respect to the coordinate axis x. The term A is the normalincidence reflection coefficient
where is the vertical Pwave impedance, Vp11 is the vertical velocity (or velocity in the isotropy plane) of the P wave, Δ denotes the difference between the values of a parameter below and above the reflecting boundary, and the bar ... indicates the mean of these values.
The coefficient B () is a socalled AVO gradient, which can be written (Rüger, 1998)
where 0 is the angle of the symmetry axis with the xaxis. The term BISO is the AVOgradient isotropic part (equal to the AVO gradient for isotropic media), and BANI is the anisotropic part of the AVO gradient.
The last term in equation (1), that is, the term with the coefficient C(), gives a marked contribution to the reflection coefficient value only for distant offsets (large incidence angles). The coefficient C() can be written (Rüger, 1998)
Following the terminology of Rüger (1998) and Tsvankin (1997), the Thomsenstyle anisotropy parameters ε(v), δ(v) and g are used to describe the anisotropy of HTI media, which are effective models of vertically fractured rocks. The parameters ε(v), and δ(v) are negative, γ is positive for HTI medium (Bakulin et al., 2000), and they are equal to zero for an isotropic medium. The difference in anisotropy parameters across the boundary is written with the help of Δ in (4) (5) as above, and the subscripts |1 and |2 refer to the upper or lower medium respect to the boundary.
The main problem is to estimate the symmetryaxis angle 0. The technique which we propose is based on equations (1) (5). Note that equation (1) is intended for calculation of reflection coefficients, while in real data, one deals with amplitudes of reflected waves, not with reflection coefficients. This brings some problems. We may not use instantaneous amplitude due to its high variability. Instead, we must take into consideration some integral characteristic of a signal for imitation of amplitude. A proper definition of this integral amplitude from wave samples of seismic records is required because amplitude estimation is very sensitive to noise and errors. Depending on this definition, the sign of amplitude may sometimes differ from the sign of the reflection coefficient; therefore it is necessary to know the sign of the difference of impedances across the boundary (see equation (2)). Also, it is necessary to correct the amplitudes for geometrical spreading, which procedure is not simple or exact for anisotropic media. Our solution of these problems is presented in Appendix A.
The algorithm for estimating the symmetry axis angle 0 may be divided into three steps. The first step consists in estimating the coefficients A, B and C in equation (1). The second step is defining two principal symmetry directions of HTI medium using the formula (3) for the AVO gradient B. The third step is analyzing the coefficient C, equation (4), to distinguish the symmetry axis direction from the fracture strike.
The first step of the algorithm is the estimation of A(), B(), and C() for each individual azimuthsectored CMP gather. The coefficients A, B and C in equation (1) are determined by the leastsquares method. For each sourcetoreceiver line j ( j = 1,...,n ), the functional
is to be minimized. Here i is the number of the offset (corresponding to the incidence angle θi) in the jth azimuthsectored CMP gather. The coefficients A, B and C are determined from the system of three linear equations
This procedure is repeated for each j (j = 1,...n ). Due to azimuthal anisotropy of the medium, the AVO gradient B acquires different values B1 B2, ..., Bn for n azimuthsectored CMP gathers.
The second step of the algorithm consists in the estimation of the angle 0 from equation (3). For this purpose, at least three AVOgradient values are required, each of them estimated for a different sourcetoreceiver line azimuth (or sector). For n=3 the problem was solved by Mallik et al. (1998). For n>3 the problem is solved by the leastsquares method as follows.
As the A value should not depend on azimuth, but the A values may be calculated different for different azimuths, then we should normalize the AVO gradient B as Bj / Aj = Bj , where Aj and Bj are calculated above by the leastsquares method for each azimuth sector j. Formula (3) for the normalized AVO gradient can be rewritten by trigonometric transformation in the form
where j is the mean azimuth of the jth azimuth sector for real data, or the azimuth value for synthetic data. The values B = B/ A are already known as estimated in the first step. The parameters a, b, and 0 are unknown, and can be obtained by minimizing the functional
This yields the system of equations
From (10), we derive an analytical trigonometric equation for 0 (given in Appendix B) and we find that the solution 0 has a period of 0. This value of the period means that we cannot be sure whether Φ0 is the symmetryaxis azimuth or the orthogonal direction, i.e. the fracturestrike. Also from equations (8), one can infer two different solutions, namely (a, b, 0) or (a, b, 0 ± , see also Zheng et al. (2004). These solutions are the local minima of the functional (9), as may be confirmed by evaluating the second partial derivatives of the functional f.
How to distinguish fracture strike from symmetry axis
For some cases, the ambiguity in fractureorientation detection may be resolved by analyzing the AVOgradient dependence on , equation (3). As shown by Hall & Kendall (2000), and by Zheng et al. (2004), when the sign of BANI is known beforehand it is possible to solve the problem. However, in general, the sign of BANI may be arbitrary. For many anisotropic media BANI > 0. But as shown by Chichinina et al. (2003), anisotropic media with BANI < 0 are not exotic, and one can be sure that BANI > 0 for a given HTI medium as long as, for hostrock velocities, the ratio VS/VP > 0.56 is valid.
Let us rewrite equation (4) as
where j = 1,...,n, Δ η = Δ δ Δ ε, and Δ ε is given by equation (5).
When substituting the value 0 ± instead of 0 into equation (11), the sign of the second term in the sum, Δ ε cos2(j 0), switches to the opposite sign, because equation (11) takes the form
Note that Δ ε should be negative for the upper reflecting boundary of the HTI layer (the case of HTI layer overlain by an isotropic overburden), and positive for the lower reflecting boundary (i.e., for the interface HTI medium isotropic medium). This follows from equation (5), taking into account ε(v)= 0 in the isotropic medium and ε(v)< 0 in the HTI medium.
Thus, the sign of Δε is predefined for the model. If, for a given 0, the sign of Δε from equation (11) does not satisfy the condition for the boundary of HTI layer mentioned above, it means that 0 is the fracturestrike direction, because the sign of the term Δεcos2(j 0) is changed. Thus, if the sign of Δε could be estimated, the ambiguity problem would be solved.
At first, we will normalize both sides of equation (11), dividing by A, as was done for equation (8), and assume from now on that C = C/ A, Δα = Δα/Aj. Δε Δε/Aj, and Δη Δη /Aj.
Then the value of the new normalized Δ ε can be estimated by minimizing the following functional f1, using the leastsquares method:
where 0 is one of the known solutions of system (10).
The system of equations (13) is solved by analogy with the system of equations (6) (7).
Keeping in mind the sign requirements for nonnormalized Δε, we suggest the following criterion for identification of the solution 0 as the symmetryaxis azimuth:
for upper reflecting boundary of the HTI layer
for lower reflecting boundary of the HTI layer
where the value of normalized Δ ε is assumed to be calculated from equations (13).
In (14) and (15), the term sign (A) is the sign of the normalincidence reflection coefficient A, which can be determined from the impedance change over the boundary, equation (2), or from equations (7). If one deals with data of amplitudes (not with data of reflection coefficients), then one can not be certain of using sign (A) from the system of equations (7), because the sign of the amplitude might not match the sign of the reflection coefficient. For example, this might arise from distortion of wave form that could lead to a changed sign of integral amplitude (see Background, and Appendix A).
Examples of algorithm application
Tests with synthetic data. The synthetic data are represented by three CSP seismograms corresponding to three sourcetoreceiverline azimuths 0°, 45° and 90°. These synthetic seismograms were generated by the ray method (Obolentseva and Grechka, 1989). Figure 1 shows one seismogram (for the azimuth of 45°) as an illustration. The model consists of 11 layers and describes the subsurface of the UrubchenTokhomo Area (Eastern Siberia), in which the HTI layer is the tenth (460 mthick Riphean carbonates). The symmetryaxis azimuth of the HTI layer is 0 = 60°. The Thomsen parameters and density for this layer as well as vertical velocities (V||S, V||P) and density for isotropic layers above and below the HTI layer are given in Table 1.
The proposed AVOA technique can be applied to the upper as well as the lower boundary of HTI layer. Here we performed AVOA analysis for the lower reflecting boundary of the HTI layer at a depth of 2665 m, and for the upper reflecting boundary of the HTI layer at a depth of 2205 m, both marked by arrows in Figure 1. Note from Table 1 a negative change of the impedance over the lower boundary, that is A < 0, and a positive change of the impedance over the upper boundary (A > 0).
For calculating the incidence angles from offsets, we used raytracing in the multilayered model, or a onelayer approximation which considers the layered medium above the reflecting boundary as a homogeneous layer. The raytracing was performed by Snell's law using the formula (Tsvankin, 1997) VP = VP||[1+δ(v)sin2θcos2θ+ε(v) cos4θ] for the velocity in the anisotropic layer, if needed.
When using the raytracing for calculating q, we estimated from (6) (10) f0 = 60.53° for the upper boundary, and 0=60.9° for the lower boundary. When using the onelayer approximation, we estimated 0 =61.26° for the upper boundary, and 0 = 61.26° for the lower boundary. Thus, the accuracy was better than 1.5% for the raytracing, and 2% for the onelayer approximation.
With a knowledge of these values of 0, we now use the algorithm for distinguishing the symmetryaxis direction by calculating the parameter Δ ε found from (13). We obtained Δ ε < 0 for both boundaries in all considered cases, which yields sign (A)Δ ε < 0 for the upper boundary and sign (A)Δ ε > 0 for the lower boundary, as it should be for the symmetry axis, equations (14) and (15).
To simulate real data, we added a random Gaussian white noise to the synthetic seismograms. Maximum amplitude of the noise was chosen as 10%, and as 20% of the maximum amplitude of the signal reflected from the upper boundary of the HTI layer at the first trace. As the amplitudes of reflected waves in the traces decrease with offset up to 2.3 times, the actual error caused by the 20% noise was up to 45% for distant offsets.
For the upper boundary, in the case of ray tracing, we obtained 0 = 63.64° for 10% noise, and 0 = 64.12° for 20% noise. For the lower boundary the respective values were 0 = 59.39° for 10% noise, and 0 = 58.03° for 20% noise. Such a difference between the results for the boundaries can be explained by interference of the waves reflected from the upper boundary with the waves reflected from the preceding boundary (see Fig. 1), the superposition interval is equal to 10 ms. To reduce the influence of interference, an additional smoothing can be applied (see Appendix A) for obtaining better results for the upper boundary.
In the case of calculation with the onelayer approximation, the results were as follows: for the upper boundary (with additional smoothing) 0 = 62.98° for 10% noise, and 0 = 63.4° for 20% noise; for the lower boundary 0 = 61.52° for 10% noise, and 0 = 57.91° for 20% noise. The onelayer approximation seems to be a good alternative to raytracing for the data used.
Application to field data. The field data contain 12 seismograms from 12 azimuth sectors (Fig. 2) for a superbin of 13x13 bins acquired at the Vancor Area (Eastern Siberia, Russia) in 2004. The target reflection is from the top of the HTI layer (overlain by an isotropic overburden) at the depth of 3.4 km marked by an arrow in Fig. 2.
The analysis of azimuthal variation in the AVO gradient provided two principal directions of symmetry: 4.3° and 85.7° (Fig. 3). From these two directions we distinguished the azimuth of 85.7° as the symmetryaxis direction. We used our criterion taking into account that the coefficient A was negative due to negative impedance change across the reflection boundary, ΔZ < 0 (estimated from the sonic logs). We substituted the calculated values of Cj and 0 = 4.3° into equations (13), and obtained the parameter Δ ε = 19.4. Following equation (14), the condition for the symmetry axis at the boundary "isotropic medium HTI medium" requires sign (A)Δ ε < 0, but in this case sign(A)=1, and therefore sign (A)Δ ε > 0. We concluded that 0 = 4.3° was the fracturestrike direction, and the symmetryaxis azimuth was 85.7°.
Discussion and conclusion
We solve the problem of ambiguity by analyzing the sign of the parameter Δ ε which is a coefficient in the second term of (11) for the coefficient C at the third, highangle, term in Ruger's equation (1). We realize that it is hard to obtain a reliable estimation of Ruger's highangle term from real field data. However, actually only the sign and not the value of Δ ε is required for the criterion (14) (15).
For determining the sign of Δ ε, we use equations (13), in which we take into account only sectors j with reliable values of coefficients Cj. For instance, in the above example of field data, the data at the azimuth sector of 7° (j=1) which give an unaccountably large value of Cj was omitted as shown in Figure 3, where the omitted points (from the polynomial fit as well as from the linear fit) are marked by crosses. The point was eliminated from the total of 12 points in the plot and only 11 points remained for the fit. Note that the corresponding AVOgradient value B for the first sector, j=1, exhibits another unreliable value like the Cj value, and should be excluded from the cos2( 0 ) fitting for the variation B (0 ) .
Following equations (1), and (6) (7), it is obvious that the estimate of Cj is linked to the Bj estimate. By selecting a reliable Bj estimate, we can provide a reliable Cjestimate. To select a reliable Bjestimate, if necessary, we perform the linear fitting y=A+Bx (where xsin2θ), additionally to the polynomial fitting y=A + Bx + Cx2, which is carried out in the first step, equations (6) (7). Figure 3 shows the trend of the polynomial fit similar to the trend of the linear fit, that gives a difference of 11.8° in the 0estimate (0 =7.5°). It seems sufficient to consider the used data as reliable.
Also, Ajestimate demonstrates close values for reliable sectors, and extremely unlike values for unreliable sectors.
We believe that the technique is applicable to field data by using controlled reliable values of C, that are consistent with reliable values of Bj and Aj.
We thank NPO EniseyGeofizika, Russia, for the field data provided, and the Institutes of Petroleum Geology and Geophysics of SB RAS, and VNIIGEOFIZIKA, Russia, for access to the synthetic seismograms.
Appendix A Using amplitudes instead of reflection coefficients
While the background of AVOA analysis is based on Ruger's formula for the reflection coefficient, equation (1), in real data AVOA analysis should use signal amplitudes. It is true that the signal amplitude is not equal to the reflection coefficient. Moreover, no picked instantaneous amplitude (sample) in the signal can be used because the signal changes form during propagation for many reasons. One should use an integral amplitude characteristic of the signal which adequately corresponds to the reflection coefficient. Let's call this characteristic simply by amplitude and denote it as P.
From our experiments, a simple definition of the amplitude as an extremum of samples in the time window of the signal , where Pn is a sample, n is a number of the sample, n1 and n2 are limits of the time window) was found nonsatisfactory.
We found that the estimated value of 0 is very sensitive to the definition of P, especially for data with noise. We suggest the following procedure for definition of P which gives good and stable results. The procedure calculates a maximum value of the signal envelope with a sign of the central peak of signal. In calculating the envelope, the Fourier transform of this signal is used: F = F++F-, where F+ is the part of spectrum corresponding to positive frequencies (ω>0), and F- is the part of negative frequencies. The envelope of the signal is given by the absolute value of inverse Fourier transform of 2F+, with F- 0 (Sheriff & Geldart, 1983).
Limits of the signal (n1 and n2) are very important parameters in this definition. In the examples above, for a Ricker wavelet which has a central peak and two adjacent peaks of opposite sign, we define the time limits of the signal by including the central peak and 0.9 of the adjacent peaks. For finding the limits for data with noise, we smooth the signal by a sequence of cubic polynomials, each defined on 12 points, using the leastsquares method.
For accuracy of Fast Fourier Transform one should use as many points as possible. For this, we calculate the envelope of a wave composed from the signal repeated five times, and we find the maximum from the middle part of this envelope.
Equation (1) should be rewritten for using the amplitudes. The amplitude of reflected PPwave can be expressed as
where cg is the coefficient of geometrical spreading (divergence) of this wave, cg = cg (θ ,), Pini is the amplitude of the source (the initial amplitude), and R is the reflection coefficient R = R (θ,) in the equation (1).
The amplitude for the normalincidence wave can be written as
where cg0 is the normalincidence coefficient of geometrical spreading, which does not depend on (θ,), and A is the normalincidence reflection coefficient, A=const, equations (1) (2). Then the reflection coefficient can be expressed as
Therefore the equation (1) for the reflection coefficient R transforms into the following equation for the amplitude P:
where B1=BP0/A, C1=CP0/A, and rg cg0 /cg . This equation should be used in the AVOA technique instead of (1).
Note that cg can be expressed as cg =cm (θ,)/r, where r is a wave path from source to receiver, and cm depends on the direction of wave propagation (for isotropic media cm=const). In assuming a weak anisotropy, one may assume a weak dependence of geometrical spreading on incidence angle: cm≈ const for a given sourcetoreceiver line with azimuth f. Then
where z is the normalincidence ray path, and cg0 ,cm /z. It is the approximate formula for divergent correction.
For multilayered media, the expressions for divergence correction can be found in Newman (1973). A practical methodology for the Pwave geometricalspreading correction in layered azimuthally anisotropic media can be found in Xu & Tsvankin (2004).
Appendix B Formula for calculating symmetryaxis direction
There is a distribution of n Bj values ( B1, B2 ,., Bn ) given for n azimuth angles j (1 ,2, ..., n). The problem is to fit the function
to these Bj values (j=1, 2, ..., n).
In this function, there are three unknown quantities: a, b, and 0 which may be determined by the leastsquares method, by introducing the functional f, equation (9). The goal of AVOA analysis is to determine the symmetryaxis direction, or to estimate the angle 0. The condition for minimum of the functional f yields the system of equations (10), from which the simple trigonometric equation for 0 is derived by excluding parameters a and b:
where KF0U0A0, L G0U0B0, MC0A02, N D0A0B0, and P E0B02, with U0, A0, B0, C0, D0, E0, F0, and G0 are determined as the arithmetic mean values (for example,
Bakulin, A., V. Grechka and I. Tsvankin, 2000. Estimation of fracture parameters from reflection seismic data Part I: HTI model due to a single fracture set. Geophysics, 65, 17881802. [ Links ]
Chichinina T., V. Sabinin and G. Ronquillo Jarrillo, 2003. Numerical modeling of Pwave AVOA in media containing vertical fractures / Expanded abstract. The Sixth International Conference on Mathematical and Numerical Aspects of Wave Propagation "Waves' 2003", Finland, Springer, 897902. [ Links ]
Hall, S. A. and JM. Kendall, 2000. Constraining the interpretation of AVOA for fracture characterization: in Anisotropy 2000: Fractures, converted waves, and case studies. Soc. Expl. Geophys., 107144. [ Links ]
Mallik, S., K. L. Craft, L. J. Meister and R. E. Chambers, 1998. Determination of the principal directions of azimuthal anisotropy from Pwave seismic data. Geophysics, 63, 692706. [ Links ]
Newman, P., 1973. Divergence effects in a layered earth. Geophysics, 38, 481488. [ Links ]
Obolentseva, I. R. and V. Yu. Grechka, 1989. Ray method for anisotropic media (algorithms and codes). Institute of Geology and Geophysics Press, Novosibirsk, 226 p. (in Russian). [ Links ]
Rüger, A., 1998. Variation of Pwave reflectivity with offset and azimuth in anisotropic media. Geophysics, 63, 935947. [ Links ]
Sheriff, R. E. and L. P. Geldart, 1983. Exploration Seismology, Volume 2, Dataprocessing and interpretation. Cambridge University Press, Cambridge, 400 p. [ Links ]
Tsvankin I., 1997, Reflection moveout and parameter estimation for horizontal transverse isotropy. Geophysics, 62, 614629. [ Links ]
Xu, X. and I. Tsvankin, 2004. Geometricalspreading correction for Pwaves in layered azimuthally anisotropic media. 74th Ann. Internat. Mtg. Soc. of Expl. Geophys., 111114. [ Links ]
Zheng, Y., D. TodorovicMarinic and G. Larson, 2004. Seismic fracture detection: ambiguity and practical solution. 74th Ann. Internat. Mtg. Soc. of Expl. Geophys., 15751578. [ Links ]