Abstract
Some new computational techniques are suggested for estimating symmetry axis azimuth of fractures in the viscoelastic anisotropic target layer in the framework of QVOA analysis (Quality factor Versus Offset and Azimuth). The different QVOA techniques are compared using synthetic viscoelastic surface reflected data with and without noise. I calculated errors for these techniques which depend on different sets of azimuths and intervals of offsets. Superiority of the high-order “enhanced general” and “cubic” techniques is shown. The high-quality QVOA techniques are compared with one of the high-quality AVOA techniques (Amplitude Versus Offset and Azimuth) in the synthetic data with noise and attenuation. Results are comparable.
1. Introduction
Explicit numerical techniques of AVOA analysis (Amplitude Versus Offset and Azimuth) are well known for determining the direction of fractures (see [1]). They are based on the equation for approximating the reflection coefficient at the boundary between elastic anisotropic layers [2–4]. Implicit numerical AVOA techniques leading to the solution of nonlinear systems of equations have been proposed in [5], and their effectiveness has been shown therein.
However, the AVOA methods do not take into consideration the presence of additional viscous attenuation of seismic waves in target oil layers. In this context, the interest is considering for oil layers the function of quality factor , which characterizes the magnitude of this viscous attenuation. The equation for dependence of factor on the direction of fractures was proposed by Chichinina et al. in works [6, 7] and serves as the basis for QVOA techniques (Quality factor Versus Offset and Azimuth).
Generally speaking, the approximate analytical QVOA equation (see (1) below) is similar to the analogous AVOA equation. However, numerical techniques for its solutions can have their own characteristics. Therefore, to make a conclusion about the adequacy of QVOA methodology to determine the direction of the fractures, one should consider different numerical techniques and test them in a strict numerical experiment.
A number of efficient numerical QVOA techniques were proposed in [8]. The present paper proposes several additional techniques. The qualitative and quantitative properties of the techniques are investigated using numerical experiments. Also, I compare the results with the best AVOA technique from [5] with the use of synthetic seismograms with attenuation and noise.
2. Models of 3D Seismic Data
In the QVOA techniques (as in AVOA), one uses 3D seismic data obtained from the receivers located on the ground surface as if at the points of a rectangular grid. The symmetry axis azimuth of the target layer fractures is evaluated for a rectangle surrounding the point of the grid (for a bin), taking into account only those traces for which the midpoint (MP) is inside the bin. If there are not enough points for the accuracy of the calculations, the adjacent bins are combined into a superbin and calculations are carried out for it. Therefore, the preliminary stage of evaluation is extracting from the 3D data all traces whose MPs lay inside the superbin. QVOA techniques are applied to seismic traces selected for a single superbin.
In the numerical experiment below, the location of the receivers on the surface is more ideal: at the points of the polar grid. I have made this to apply the 2D numerical model for generating 3D seismograms by rotating it around the source.
3. Overview of QVOA Techniques
If fractured rocks with a preferential orientation of fractures are saturated with a fluid, then the fluid flows may lead to azimuthally varying attenuation of seismic waves. Let us define the attenuation factor (where is a complex velocity [7]) that corresponds to the definition of the quality factor , where is the wave energy and is the energy lost per cycle due to attenuation [9]. If we set , it is possible to derive the approximate equation [7] for the attenuation factor in the following form:where is the azimuth and the incidence angle is the angle relative to the normal vector in the target layer.
This approximation was made under the assumption that the term with is negligible.
According to [7], for HTI media,where is the symmetry axis azimuth of the fracture-strike direction.
Five numerical techniques for estimating the symmetry axis azimuth using QVOA were proposed in [8]. Three of them, the general technique (G), the truncated technique (T), and the sectored technique (S), are based on (1)-(2).
In technique G, the equation that is obtained by substituting (2) into (1) is solved by the least-squares method:
In technique T, a truncated version of (3) without the last term is solved.
Technique S solves in turn (1) and the first equation of (2). For this technique, the seismic data of the superbin should be separated by sectors. Equation (1) is solved for each sector independently to obtain the function , and it is assumed that in each sector —the middle azimuth of the sector.
The least-squares method is used to obtain in these three techniques, but it gives nonlinear systems of equations that have no analytical solution. To deal with the analytical solutions, the first equation of (2) can be replaced by the following approximate equation in techniques T and S:
As a result of this replacement, two additional techniques arise, approximate truncated technique (AT) and approximate sectored technique (AS), which are similar to the linear and sector AVOA methods (see [5]). They have the analytical solutions.
4. Generalized and Enhanced QVOA Techniques
The practice of applying the approximate techniques (AT and AS) showed that they provided more preferred results than corresponding T and S techniques, based on (2), for certain seismic data. Perhaps, this is due to the fact that (4) has one parameter more () than (2) and therefore a more general form than (2).
Let us try to replace (2) by the following generalized equations:
After formal substitution of (5) into (1), generalized equation is obtained instead of (3):where and .
Equation (6) is the same as the equation for the general method of AVOA [5], except for the replacement of the amplitude function in it by attenuation factor . Let us call this new QVOA technique based on (6) the generalized technique (Gd). The solution algorithm for it is given in Appendix.
To account for nonlinearity in the sectored technique, one can base it on (1), function , and the second equation of (5). Here, the least-squares solution gives a system of four nonlinear equations, which is simpler than technique Gd. Let us call this new technique SC—sectored on C.
As shown in Figure 2, the dependence of on can differ from the quadratic dependence prescribed by (6). There are two ways to take this discrepancy into account. The first is to formally add a term in (6). This leads to solving a nonlinear system of 11 equations in the least-squares method. The solution of this system is not given here due to its bulkiness. It can be derived by analogy with the Appendix. This new technique will be called cubic (C).
The second way is to improve the quality of the approximate equation (1). It was obtained in [7] from a more complex equation of the following form:
Substitution of generalized equations of the form (5) into (7) gives
Application of the least-squares method to (8) gives a nonlinear system of 12 equations, which is solved similarly to Gd (see Appendix). The improved technique thus obtained (which is an enhanced version of the generalized technique) will be called EGd—enhanced Gd.
It is also worth considering a truncated version of (8) for simplicity:
This leads to solving a system of 9 nonlinear equations in the least-squares method. The technique, based on (9), will be called the truncated enhanced Gd technique (TEGd).
I do not cite here the systems of equations for EGd and TEGd techniques, because they can be derived by analogy with the Appendix.
Finally, to complete the view, it is necessary to consider the substitution of theoretical equations (2) into (7). This gives two more techniques that will be called by analogy: the enhanced general technique (EG) and the truncated enhanced general technique (TEG).
All 12 numerical QVOA techniques for evaluation of the symmetry axis azimuth are listed in Table 1 for reference.
There are other possible methods based on (1), but to clarify the question about the applicability of the QVOA methodology, the 12 already presented ones should be enough to us.
5. Calculation of Attenuation Factor
Correctly determining the attenuation factor (or quality factor ) based on seismograms is very important for the accuracy of QVOA techniques. At the time, this issue has been given a lot of attention, and many different ways were invented. I adhere to the following procedure.
The process of attenuation of the wave amplitudes caused by dispersion in the medium is subjected to the law (see, e.g., [10]):where is the distance traversed by the ray, is time, and is the absorption coefficient.
The attenuation factor is associated with as follows [11]:where , is frequency, and is travel time.
The ratio of spectral amplitudes of the wave impulses reflected from the bottom () and from the top () of target layer with attenuation (see Figure 1) can be expressed from (10) and (11) as follows (see, e.g., [12]):where is a coefficient that includes reflection coefficients and geometrical spreading and is travel time within the target layer.
Under the assumption of weak dependence of on the frequency , (12) can be rewritten in the following form:where .
In order to reliably calculate , (13) must be solved for a sufficiently representative interval of frequencies, where the magnitude of is approximately constant (see Figure 1). The solution of (13) can be obtained by the least-squares method or the method of frequency shift [13] or simply taking the numerical derivative from with respect to . For example, using a 9-point derivative gives the following formula:where is the step of discretization and and are the boundary indexes of the interval.
The value of is calculated from (11), where .
The quality of the spectral ratio significantly affects the computed value . Therefore, some of the usual conversions, which visually enhance the shape of the impulse but actually spectrally distort it, should be avoided. Examples of operations that should be avoided include multiplication of impulse by window functions to reduce the influence of the time window borders and smoothing impulse by filters which was offered for AVOA in [5].
The boundaries of the time window of the impulse should be taken as widely as possible, while still considering the noise magnitude. This is a very sensitive parameter. In the following computational experiment, I selected the boundaries at the 0.15 level of the maximum of the envelope (i.e., by the formula ) in the case of the seismograms with noise and by the formula in the case of seismograms without noise. For calculating these envelopes of noisy impulses, it is necessary to smooth the impulses with filters.
To calculate the spectral interval , I applied the adaptive algorithm that calculated a minimum of the average second derivative from with respect to (of the mean curvature). In the experiment, the algorithm varied the length of spectral interval from till of the peak frequency of spectrum of the impulse from the top (), and the position of the spectral interval between and of .
Important in the calculation of is the value of the travel time of the ray inside the target layer. I applied the algorithm for calculating by the ray method from [14], which also yielded the value of the incidence angle . It uses a correlation function between the envelopes of the smoothed impulses from the top and bottom of target layer on the same trace. Unfortunately, this algorithm may give erroneous results in the case when the impulse shape remains substantially deformed after the smoothing.
6. Numerical Experiment for Comparison of the QVOA Techniques
The 12 techniques from Table 1 were compared in terms of the ability to give the most accurate value of in an anisotropic HTI medium.
I generated surface synthetic seismograms of P-reflections in a three-layer medium with anisotropic viscoelastic middle layer, using the modified method of [14] for 2D viscoelastic wave modeling. I derived 2D models of the anisotropic layer for different values of by rotating around axis the two-dimensional stiffness tensor (see [7, 15]) for the plane in anisotropic HTI layer. I set necessary anisotropic parameters for defining the tensor so and (see [15] for definitions). I set the following relaxation times needed for the simulation of the viscoelastic wave propagation: and sec (see [14] for definitions). They gave the values of the factor in the range of 16–26. One can see the example of calculated seismogram in [8].
Such replacement of the full 3D picture of waves by a “fan” of 2D pictures of waves is undoubtedly an approximation and should worsen the results of evaluation of . I came up with a way of correcting these results. It is described in the next section.
Velocity in layers from top to bottom was 3200, 4000, and 4800 m/sec (another variant for the third layer, 3200), was half of , densities were constant, and capacity of the top two layers was 1600 and 400 m. The source of explosive type generated one Ricker impulse of frequency of 30 Hz. The receivers were spaced every 100 m, starting from the source, and measured the -component of particle velocity. In each seismogram, there were 50 traces corresponding to 50 offsets, including zero. However, since the traces with a normal incidence angle should not depend on the azimuth, they were excluded from consideration.
Just as in [5], I calculated for different intervals of offsets from the minimum offset uniformly to the maximum one. In one set of intervals, the minimum offset was fixed at 2 (first offset), and the number of maximum offsets was varied from number 50 to number 3. In the other set of intervals, the maximum offset was fixed at the number 50, and the minimum offset was varied from number 2 to number 48. Of course, the maximum angle of incidence , corresponding to the maximum offset, and the minimum angle of incidence , corresponding to the minimum offset, were also accordingly changed in these sets of offsets.
I got different sets of synthetic seismograms for different azimuths. It is obvious that the sets of azimuths, which have a symmetry with respect to azimuth , will give better results than the sets without such symmetry, due to the repetition of the data of seismograms. Therefore, I selected two sets of azimuths: asymmetrical set , computationally the most unfortunate, and symmetrical set .
To better approximate to real data, I added a random normal Gaussian noise to synthetic seismograms—different for different seismograms. The maximum amplitude of the noise was 10% of the maximum amplitude of the impulse reflected from the top of target layer at the first trace of seismogram.
As the results did not depend significantly on the preset values of , I accepted for simplicity .
7. Correction of Algorithms on Data Nature
Replacement of 3D model of an anisotropic viscoelastic medium by a simplified model of rotation of 2D medium, which was done to save computer resources, without a doubt, should lead to a distortion of the simulation results and can give different dependence on for the coefficients of (1) than the theoretical dependence (2). To understand and correct this effect, let us consider the results of calculation of attenuation factor .
Figures 2 and 3 show the dependence of the calculated values of on for different values of from the asymmetric set of azimuths. Figure 2 is for variant (i.e., ). Figure 3 is for variant (i.e., ).
It is seen that these curves are more complex than quadratic dependence (see (1)). This can cause large errors in rough approximation techniques, such as T, S, AT, and AS.
In Figure 3, the curves show lower smoothness than in Figure 2, especially at far offsets, . However, the dependence of the curves on the azimuth is apparently not affected significantly, as we will see later by comparing Figures 7–9 with Figures 4–6.
(a)
(b)
(a)
(b)
Let us look more closely at Figure 2. Curves are gathering to the value stronger than to zero, and the average curve 45° is everywhere shifted from the middle of the bulb of curves (where it must be on straight sections of curves according to (2)) approximately on . In Figure 3, the nature of the condensations and displacements of curves are similar, though less pronounced. Since this property of data is not consistent with (2), we must admit the need in correction of (2) with respect to our synthetic seismic data.
The simplest way is to introduce a correcting shift in formulae (2) or in (5) if it is applied:where is the correction angle, which can be different for different methods and for different ratios and .
The idea was that the curve of error of on calculated value of has a minimum at the full 3D data at and, to obtain this minimum in quazi-3D data, it is needed to shift this curve in axis . If one had 3D but not 2D seismograms, one would not need to apply the shift at all. One can evaluate the value of shift by a simple enumeration of values of with finding the minimum of value of average error on intervals of offsets. However, it is a noneconomical procedure.
Therefore, the shift values in the range were pretested. It was found that value was suitable for techniques Gd, SC, C, EGd, and TEGd, and value was suitable for technique G. Value gave good results for techniques EG and TEG in the case of seismograms without noise, but in the case of seismograms with noise, it was well for TEG and for EG.
For techniques S, AS, T, and AT, I used for different variants different values of —those obtained with the enumerating procedure mentioned above. In principle, this way of definition of can be applied to all techniques (with excellent results), but it is too uneconomical.
In implicit (nonlinear) techniques, the shift has led not only to lower error of solutions but also to a more correct operation of the algorithm for determining a unique solution from the set of possible ones (see Appendix).
8. Results for Seismograms without Noise
In Figures 4–9, results for sets of seismograms without added noise are presented.
In Figures 4–6, the error of calculated values of in degrees (difference with the correct value 0) for the asymmetric set of azimuths is shown and . In Figures 4 and 5, the techniques in which error is formally above 4° at a substantial interval of offsets are presented. They can be called techniques of low quality. Figure 4 is for a fixed (the minimum, second offset), and Figure 5 is for a fixed (the maximum, 50th offset). Figure 6 presents the techniques of higher quality: case (a) for the fixed minimum (second) offset and case (b) for the fixed maximum (50th) offset.
In Figures 7–9, similar results are shown for variant .
From Figures 4–9, one can conclude that the techniques S, AS, T, AT, G, SC, EGd, EG, and TEG exhibit, in varying degrees, the cases and intervals of unstable errors, of big errors, and of small errors. In particular, T, G, SC, EGd, EG, and TEG may be sometimes included in the category of higher-quality methods (Figures 6 and 9).
The only techniques in which the error is small and stable for all sufficiently large intervals of offsets are Gd, C, and TEGd.
All techniques are characterized by instability of the errors at far intervals of offsets. This possibly is due to numerical errors at the top of target layer which were introduced by the finite-difference approximation used for obtaining seismograms, because at this interface at far offsets the angle of refraction reaches and exceeds the critical value.
For the symmetric set of azimuths, the error is 0.0000 for all techniques as in the variant and as in the variant in the absence of noise. This can be explained by the repetition of the data on symmetric azimuths and by good quality of the synthetic seismograms.
9. Results for Seismograms with Noise
The noise significantly affects the accuracy of the techniques. On the seismograms, one can see that the 10% noise in the impulses of reflections from the top at the first trace (the absolute magnitude of which was set the same for noise at all traces in the seismogram) really becomes the 100% noise at far traces of the seismogram in impulses of reflections from the bottom (see [5] for illustration) because of decreasing amplitude of impulse. This huge 100% noise, of course, introduces large errors in the computed quantities, spectra, incidence angles, and travel times, and thus in the value of .
Unfortunately, smoothing impulse by filters to decrease the noise did not improve the calculation of factor , as mentioned above. Therefore, this was not applied here.
As the noise was set different on different seismograms, the effect of repetition of data in the symmetric set of azimuths was missing, which deprived it of the advantages over the asymmetrical set.
In Figures 10 and 11, the results for the best three techniques and the asymmetric set of azimuths are presented.
(a)
(b)
(a)
(b)
The cubic technique (C) shows, perhaps, the best behavior. It gives the error of less than 5° and is quite stable across all incidence angles, almost in all graphs. Techniques TEG and EG are successfully competing with it.
Techniques Gd and SC, not shown here, give a little more error but are less stable. Other methods are less successful.
10. Comparison with the AVOA Method
It is not obvious that the AVOA methods will give good results in predicting , if applied to the seismograms with viscous attenuation, since they are based on the equation for elastic media. The seismograms used in [5] were “elastic” (without viscous attenuation). In this section, I use the seismograms with noise, obtained for the case of viscoelastic target layer (with attenuation) to test the best AVOA method from [5] and compare it with the best QVOA techniques.
Graphs of the amplitude function (which is used in the AVOA method instead of ), similar to Figure 2, show much smaller dependence on noise than the graphs for the attenuation factor , and therefore one can expect less accuracy of QVOA techniques. However, it is not quite so.
In Figures 12 and 13, I present the error of calculated for the best QVOA techniques (C, EG, and TEG) and AVOA General method (see [5]) for the seismograms with noise and the asymmetric set of azimuths. General method AVOA was used as for the top of target layer (marked with ) and for the bottom one (). Figure 12 shows the variant , and Figure 13 shows the variant . In the AVOA method, the correction was also used with the shift , and in the variant , the smoothing of amplitudes was made by the method of [5].
(a)
(b)
(a)
(b)
It should be noted that the results of QVOA techniques lay mostly between the results of the AVOA method for top and bottom. The error of , in general, is more than the error of , except Figure 12(b). This can be explained by an additional change of amplitude due to viscosity inside the target layer and thus by increase of the AVOA error. The corresponding results of AVOA in the seismograms without attenuation show a much smaller error for the bottom.
In this part of experiment, the results of AVOA and QVOA look comparable, except some instability of the AVOA results at the bottom.
11. Discussion and Conclusions
The different behavior of curves in the case (Figure 2) and in the case (Figure 3) may be linked to the achievement of the critical angle of refraction. However, nevertheless, the evaluation of the symmetry axis azimuth with different techniques is qualitatively in agreement with these curves. The distortion of computed is typical at far incidence angles, which determines the area of instability of the error of most techniques at , or . It can be seen in Figures 4–9.
Significant nonlinearity of curves in (Figures 2 and 3) is the cause of large errors of truncated techniques (T and AT) and the line sector techniques (S and AS). Techniques SC and G have smaller errors, since they take into account the nonlinearity.
The techniques of higher quality, Gd, C, EG, TEG, EGd, and TEGd, show better results because they are based on nonlinear equations; the higher the nonlinearity of the base equation (the cubic technique C) is, the better the results are.
It is interesting to note that techniques C, EG, and TEG have the most stable error on the offsets (Figures 4–9).
The presence of noise requires smoothing seismograms. Unfortunately, smoothing by filters [5] suppresses not only the high frequencies of noise but also the low frequencies of impulse. The remaining distortion of the impulse shape is a source of errors when calculating the spectrum, spectral ratio, spectral interval, travel time, and, ultimately, the attenuation factor. Such a large number of errors introduced into the QVOA techniques can give unstable and incorrect calculation of the symmetry axis azimuth, if the amplitude of the noise is large. As such, I gave up smoothing the impulses in the calculation of the spectra.
In Figures 10(a) and 11(a), the area of instability is present when . Large errors at small can be explained from (3). In (3), depends on the product . Error in gives error in which is greater when is less.
For the seismograms with noise, the best results are given with techniques which gave the most stable error in the absence of noise: C, EG, and TEG.
A comparison of these techniques with the best AVOA method at viscous attenuation (Figures 12 and 13) shows that they are generally a little worse than , that is, the general method of AVOA applied to the top of target layer. But, in some cases, they are slightly better than . This means that the AVOA method is not the best for the viscoelastic target layer if it is applied to the bottom of the layer, where, despite large opportunities for errors, the best QVOA techniques are not worse than it. One can say that the amplitudes (or rather the reflection coefficients) in the presence of attenuation are not described quite well at the bottom of target layer by Rüger’s equation [3] which is the base for methods [5].
The obtained superiority of techniques EG and TEG based on the improved equation of Chichinina [7] over other techniques (except cubic) indicates that the theory of QVOA proposed in [7] is correct, since it has received here experimental confirmation.
As can be seen from Appendix, for nonlinear techniques, there is a problem of choosing one of several solutions of the nonlinear system of equations derived with the least-squares method. This problem is linked with the well-known 90° ambiguity of solution due to dependence (2) on . For seismograms without noise, this ambiguity can be resolved well by applying the empirical criterion described in Appendix. But, in the case of noise, the criterion does not give reliable choice in some cases, especially for areas of instability.
In such unreliable conditions of noisy seismograms, good strategy for evaluating the symmetry axis azimuth is the use of different techniques together: for example, C, Gd, EG, and TEG QVOA techniques together with G and L (see [5]) AVOA methods at the top of target layer.
Applied to real field data, suggested QVOA and AVOA techniques can provide less reliable results. Real data contain more of the wave interference than synthetic data do, and it is almost impossible to correctly separate one from the other superimposed waves. Such distortion in the low frequencies of the impulses can lead to unpredictable results.
Appendix
Details of the Generalized QVOA Technique (Gd)
The functional of the square error of equation (6) is ( is number of offsets in the set):
According to the least-squares method, finding the minimum of the functional, one can derive a linear system of equations for fixed from (A.1):where , , , , , , , , , , , , , , , , , , , and .
The condition of equality to zero of the derivative of the functional with respect to gives another equation of the following form:where , , , , , , , , , , , , , and .
The nonlinear equation (A.3) is solved by the method of bisecting interval within the region (−90°, 90°). Provided for each intermediate value , system (A.2) is recalculated.
As a rule, one should obtain more than one solution of equation (A.3). Selecting a single correct value of should correspond to the minimum of the minima of . However, this criterion is insufficient, because has about the same minimal value as at and as at . This 90° ambiguity of solution is a trigonometric property of dependence ((2) or (5)) on . To resolve the ambiguity, one should use additional observations and dependencies.
I have determined an additional criterion for separating the solutions from a consideration of Figure 2. When increasing , the curves of on are as if compressed toward the beginning. This means that if is the average value of over the offset interval, then has a maximum at and a minimum at . This criterion is suitable for all techniques.
For data with noise, this criterion behaves well in C, EG, and TEG techniques out of areas of instability.
Conflicts of Interest
The author declares that he has no conflicts of interest.