Seismic azimuthal gradient estimation

Abstract

Method of estimating azimuthal amplitude gradient is disclosed. This method uses a correlation of seismic attributes within a sliding volume of data to obtain azimuthal gradient.

Claims

1. A method for processing 3-D seismic data comprising: a) forming, using a computing processor, common image location gathers of seismic traces as a function of time or depth, offset, and azimuth; b) performing a linear regression on a gather at each time or depth to generate a model that includes an intercept term A, an average gradient term B.sub.avg, an anisotropic gradient term B.sub.cos, and an anisotropic gradient term B.sub.sin; c) calculating joint correlations of B.sub.cos and B.sub.sin seismic traces over a plurality of rectangular prisms in lateral position and time or depth, each prism centered about its own analysis point; d) forming a vector whose rectangular coordinates are (M, R) for each analysis point, wherein M represent half of difference between mean-squared amplitudes of B.sub.cos and B.sub.sin within the prism, and R represent the average product of B.sub.cos and B.sub.sin within the prism; e) calculating angles 2 .sub.up or 2 .sub.down the vector makes with the M axis; f) taking one-half the rotation angles .sub.up or .sub.down as symmetry azimuths of vertical reservoir fractures; and g) determining azimuthal amplitude gradient based on the determined symmetry azimuth, wherein said common image location gathers of seismic traces describe subsurface heterogeneity of a reservoir.

2. The method of claim 1, wherein the intercept term A is independent of incident angle or azimuth.

3. The method of claim 1, wherein the average gradient term B.sub.avg depends only on squared sine of incident angle.

4. The method of claim 1, wherein the anisotropic term B.sub.cos depends on product of cosine of twice source-receiver azimuth and squared sine of incident angle.

5. The method of claim 1, wherein the anisotropic term B.sub.sin depends on product of sine of twice source-receiver azimuth and square sine of incident angle.

6. The method of claim 1, wherein length of the vector is taken to be squared magnitude of azimuthal anisotropy.

7. The method of claim 1, wherein the B.sub.cos and B.sub.sin traces are made analytic by combining the traces with (1) times their respective Hilbert transforms.

8. The method of claim 1, wherein the joint correlations are formed by allowing Rc to represent average complex product of B.sub.cos times conjugate of B.sub.sin within each rectangular prism, M.sub.c to represent half of difference between mean-squared magnitudes of B.sub.cos and B.sub.sin within each prism; and P.sub.c to represent half of sum of mean-squared magnitudes of B.sub.cos and B.sub.sin within each prism.

9. The method of claim 1, wherein the joint correlation between the B.sub.cos and B.sub.sin is calculated within variable-height windows that conform in depth to geologic horizons.

10. A method for processing 3-D seismic data comprising: a) forming common image location gathers of seismic traces as a function of time or depth, offset, and azimuth; b) performing a linear regression, using a computing processor, on a gather at each time or depth to generate a model that includes an intercept term A, an average gradient term B.sub.avg, an anisotropic gradient term B.sub.cos and an anisotropic gradient term B.sub.sin; c) calculating joint correlations of B.sub.cos and B.sub.sin seismic traces over a plurality of rectangular prisms in lateral position and time or depth, each prism centered about its own analysis point; d) forming a vector whose rectangular coordinates are (M, R) for each analysis point, wherein M represent half of difference between mean-squared amplitudes of B.sub.cos and B.sub.sin within the prism, and R represent the average product of B.sub.cos and B.sub.sin within the prism; e) calculating angles 2 .sub.up or 2 .sub.down the vector makes with the M axis; f) determining one-half the rotation angles .sub.up or .sub.down as symmetry azimuths of vertical reservoir fractures; g) determining azimuthal amplitude gradient based on the determined symmetry azimuth; and h) detecting one or more subsurface heterogeneity utilizing the determined azimuthal amplitude gradient, wherein said common image location gathers of seismic traces describe subsurface heterogeneity of a reservoir.

11. The method of claim 10, wherein the intercept term A is independent of incident angle or azimuth.

12. The method of claim 10, wherein the average gradient term B.sub.avg depends only on squared sine of incident angle.

13. The method of claim 10, wherein the anisotropic term B.sub.cos depends on product of cosine of twice source-receiver azimuth and squared sine of incident angle.

14. The method of claim 10, wherein the anisotropic term B.sub.sin depends on product of sine of twice source-receiver azimuth and square sine of incident angle.

15. The method of claim 10, wherein length of the vector is taken to be squared magnitude of azimuthal anisotropy.

16. The method of claim 10, wherein the B.sub.cos and B.sub.sin traces are made analytic by combining the traces with /(1) times their respective Hilbert transforms.

17. The method of claim 10, wherein the joint correlations are formed by allowing Rc to represent average complex product of B.sub.cos times conjugate of B.sub.sin within each rectangular prism, M.sub.c to represent half of difference between mean-squared magnitudes of B.sub.cos and B.sub.sin within each prism; and P.sub.c to represent half of sum of mean-squared magnitudes of B.sub.cos and B.sub.sin within each prism.

18. The method of claim 10, wherein the joint correlation between the B.sub.cos and B.sub.sin is calculated within variable-height windows that conform in depth to geologic horizons.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) A more complete understanding of the present invention and benefits thereof may be acquired by referring to the follow description taken in conjunction with the accompanying drawings in which:

(2) FIG. 1 shows a collection 3D seismic data over a thin reservoir having vertical fractures oriented perpendicular to symmetry azimuth.

(3) FIGS. 2A-2B shows side by side comparison of azimuthal gradient polarity (FIG. 2A) and symmetry azimuth (FIG. 2B).

(4) FIG. 3A shows computation of azimuths using point-wise method. FIG. 3B shows rotation angles for window-based statistical method.

(5) FIG. 4A-4B shows conversion of B.sub.cos (FIG. 4A) and B.sub.sin (FIG. 4B) traces in .sub.max (FIG. 4C) and .sub.up (FIG. 4D) azimuth traces.

(6) FIG. 5A shows a synthetic reflectivity gather that was processed by standard methods. FIG. 5B shows corresponding anisotropy parameters obtained by point-wise methods.

(7) FIGS. 6A-6H show ideal crossplots of B.sub.sin versus B.sub.cos for various symmetry azimuths and spins. FIG. 6A shows for .sub.sym=22.5. FIG. 6B shows for .sub.sym=0. FIG. 6C shows for .sub.sym=22.5. FIG. 6D shows for .sub.sym=45. FIG. 6E shows for .sub.sym=67.5. FIG. 6F shows for .sub.sym=90. FIG. 6G shows for .sub.sym=67.5. FIG. 6H shows for .sub.sym=45.

(8) FIGS. 7A-7D show down and up interpretations of same joint distribution of =B.sub.cos and b=B.sub.sin for .sub.up=30. FIG. 7A shows plot of a vs b. FIG. 7B shows plot of the down interpretation. FIG. 7C shows plot of R.sub.r vs M cos 2. FIG. 7D shows plot of the up interpretation.

(9) FIG. 8 shows an embodiment as described in the detailed description.

(10) FIG. 9 illustrates that joint correlations can be calculated within variable-height windows that conform in depth to geologic horizons [1, 2] and interpolated horizons [3]. Each window have fixed lateral widths, fixed vertical height [6] at window's analysis point [4, 5] about which the window is centered.

(11) FIGS. 10A-10B show various synthetic results using windowed statistical method (a) before and (b) after regions analysis.

(12) FIGS. 11A-11C show divisions of (R, M) plane.

DETAILED DESCRIPTION

(13) Turning now to the detailed description of the preferred arrangement or arrangements of the present invention, it should be understood that the inventive features and concepts may be manifested in other arrangements and that the scope of the invention is not limited to the embodiments described or illustrated.

(14) The following examples of certain embodiments of the invention are given. Each example is provided by way of explanation of the invention, one of many embodiments of the invention, and the following examples should not be read to limit the scope of the invention.

(15) The present invention provides tools and methods for estimating azimuthal amplitude gradient using correlation of seismic attributes within a sliding volume of seismic data (windowed statistical method). As used herein, azimuthal amplitude gradient represents a change in seismic amplitude with respect to source-receiver azimuth or orientation of source-receiver line. In one embodiment, a new method of performing AVAZ analysis is provided. One of the goals of this method is to obtain estimates of symmetry azimuth .sub.sym that is more stable in the presence of noise and small misalignments of data from one azimuth to the next and to resolve the ambiguity between gradient polarity and symmetry azimuth. These estimations can be used to detect subsurface heterogeneities such as, but not limited to, faults, fractures, stratigraphic discontinuities, and the like. Other advantages will be apparent from the disclosure herein.

(16) In combination with well information or prior geologic constraints, the estimation azimuthal amplitude gradient can infer the orientations and magnitudes of fractures away from the well. Such information is valuable for developing oil and gas resources, particularly where artificial fracturing of reservoir rock is required. Moreover, the present invention can utilize well control and/or prior geologic information to resolve fundamental ambiguity problem generally encountered in AVAZ analysis.

(17) Ambiguity of AVAZ Analysis

(18) FIG. 1 illustrates one typical geometry associated with seismic data collection. As shown in FIG. 1, 3D seismic data can be collected over a thin reservoir having fractures oriented perpendicular to the symmetry azimuth. In such a setup, 3D seismic data can be recorded over a variety of source-receiver azimuths. Ray paths in this figure all reflect off the top of a vertically fractured reservoir interval, which is modeled as being anisotropic, but having transverse isotropy with a horizontal axis of symmetry (HTI medium). To a first approximation, amplitude of a plane P-wave (R.sub.p) reflected from such an interface varies both with source-receiver azimuth and incidence angle with respect to a ray normal to the interface:
R.sub.p(,)=A+B()sin.sup.2(1)

(19) Constant term A is normal-incidence reflection coefficient, which equals half of the fractional change of acoustic impedance across the interface. The term A depends neither on the azimuth nor the incidence angle. Coefficient B() is the amplitude gradient, and it specifies how rapidly the reflected amplitude changes with incidence angle at any azimuth . The amplitude gradient has the following two components:
B()=B.sub.iso+B.sub.ani cos.sup.2(.sub.sym)(2)
In equation (2), .sub.sym represents symmetry azimuth of the fracture sets, as shown in FIG. 1, and B.sub.ani is the anisotropic gradient. The symmetry azimuth in this model is horizontal and perpendicular to the fractures. The B.sub.iso term is the amplitude gradient in the isotropy plane parallel to the fractures, when =.sub.sym90.

(20) FIGS. 2A-2B illustrate the fundamental ambiguity problem with conventional AVAZ analysis. FIGS. 2A-2B plot reflection amplitude at a fixed incidence angle as a function of source-receiver azimuth. If B.sub.ani is positive, FIG. 2A applies. In this case, B.sub.iso must be positive and the symmetry azimuth .sub.sym is as is shown. If B.sub.ani is negative, FIG. 2B applies. Now B.sub.iso must be negative and the symmetry azimuth is separated by 90 from what it was before. However, the amplitudes as functions of azimuth are identical in the two cases. Because of this ambiguity, it is impossible to determine, from the amplitudes alone, which scenario applies.

(21) Fracture Detection by AVAZ Analysis Using the Standard Point-Wise Method

(22) In order to use a linear regression to extract the desired parameters from equation (2), namely B.sub.iso, B.sub.ani and .sub.sym, the terms can be linearized as a sum coefficients to be determined, each multiplied by a fixed basis function that is independent of the data:
B()=B.sub.avg+B.sub.cos cos 2+B.sub.sin sin 2(5)
where B.sub.avg is an average over all azimuths and where
B.sub.cos=B.sub.ani cos 2.sub.sym; B.sub.sin=B.sub.ani sin 2.sub.sym(6)

(23) Equation (5) is convenient in that three basis functions (1, cos , sin ) are mutually orthogonal over full range of azimuths. A vector of parameters
{right arrow over (X)}=[A,B.sub.avg,B.sub.cos,B.sub.sin].sup.T(7)
is estimated from linear regression using following equation:
{right arrow over (X)}(.sub.i{right arrow over (g)}.sub.i{right arrow over (g)}.sub.i.sup.T+l).sup.1(.sub.id.sub.i{right arrow over (g)}.sub.i)(8)
where I is the identity matrix, {right arrow over (g)}.sub.i is the vector of basis functions of the i-th seismic sample d.sub.i evaluated at its imaged point,
{right arrow over (g)}.sub.i=[1,sin.sup.2 .sub.i,cos 2.sub.i sin.sup.2 .sub.i,sin 2.sub.i sin.sup.2 .sub.i].sup.T(9)
and is a stability factory sometimes needed to constrain L-2 norm of {right arrow over (X)}.

(24) It can be seen from equation (5) that the azimuthal gradient is a vector quantity that takes two numbers to describe it at every reflection time and position. It can be expressed in Cartesian coordinates aligned with the seismic survey coordinates as it {right arrow over (u)}=(a, b), where a and b to be the B.sub.cos and B.sub.sin elements of the regression vector {right arrow over (X)}. It can also be expressed in polar coordinates as pointwise magnitude B.sub.mag=2{square root over (a.sup.2+b.sup.2)} and azimuth of the maximum signed amplitude gradient, .sub.max. This azimuth can then be found as half the counter-clockwise angle the {right arrow over (u)} vector makes with the -axis, as shown in FIG. 3A. FIG. 3B illustrates rotation angles for window-based statistical method. The azimuth of the minimum gradient, .sub.min, is half the counter-clockwise angle the vector {right arrow over (u)} makes with the -axis. In both cases, a clockwise angle represents a negative angle, so for in the example shown, .sub.min<0. If B.sub.ani>0 then .sub.sym=a.sub.max. If B.sub.ani<0 then .sub.sym=.sub.min. Since we do not know from B() alone whether B.sub.ani is positive or negative, an ambiguity arises in determining .sub.sym. However, in either case, the magnitude of the anisotropic gradient B.sub.mag=|B.sub.ani|, equals twice the length of the {right arrow over (u)} vector, and the amplitude gradient in the isotropy plane is given from
B.sub.iso=B.sub.avgB.sub.ani(10)

(25) Seismic data will generally have no zero-frequency component. Whatever reflections present in the earth are convolved by a zero-mean wavelet, which has equal positive and negative areas. For example, an interface having a positive B.sub.ani coefficient with a zero-phase wavelet would yield value correct of .sub.sym in central wavelet lobe, while an azimuth that is 90 opposite .sub.sym would appear in the negative side lobes of the wavelet. Since this literal interpretation of azimuth is unrealistic from a geologic perspective, it becomes critical to ensure that the polarity of B.sub.cos remains fixed through horizon extractions of .sub.max or else striping can occur.

(26) FIGS. 4A-4B illustrate this situation with perfect synthetic B.sub.cos (FIG. 4A) and B.sub.sin (FIG. 4B) traces. As shown, .sub.max (FIG. 4C) switches between two perpendicular directions as seismic polarity flips. When B.sub.sin is near zero and B.sub.cos<0, .sub.max becomes unstable but this is acceptable when understood that an azimuth of +90 is equivalent to 90 due to reciprocity. Thus, it is possible to obtain from the pointwise method a measure of azimuth that does not oscillate between two perpendicular directions. These polarity-independent azimuths can come in two types: spin-up and spin-down. Spin-up azimuths (.sub.up) are constrained to lie within 45, while spin-down (.sub.down) azimuths are constrained to lie outside this range, but inside the range 90. Spin-up azimuths can be derived from .sub.max azimuths from the following guidelines:
a. .sub.up=.sub.max if |.sub.max|<45
b. .sub.up=.sub.max90 if .sub.max45
c. .sub.up=.sub.max+90 if .sub.max45
A corresponding spin-up gradient (B.sub.up) can be defined as |{right arrow over (u)}| when |.sub.max|<45 and |{right arrow over (u)}| otherwise.

(27) Referring to FIG. 4D, .sub.up does not switch directions as the seismic polarity reverses. However, it is constrained to lie within 45. If the true symmetry azimuth is, for example, 60, .sub.up will come out to be 30. It is also evident that .sub.up is unstable when B.sub.cos is near zero and B.sub.sin is negative. This is acceptable since spin-up azimuths of 45 and 45 are equivalent. .sub.up is also unstable around zero-crossing of .sub.max. The spin-down azimuths can be obtained from the spin-up azimuths from the following:
a. .sub.down=.sub.up90 if .sub.up>0
b. .sub.down=.sub.up+90 if .sub.up<0
Synthetic Data Example

(28) To obtain a more realistic example, synthetic data for a 1-D model was generated. 141 synthetic traces with random source-receiver azimuths and offsets ranging from 0 to 3.1 km were selected. These data underwent azimuthal velocity analysis using trim statistics to align events to be as flat as possible at each time sample. The data were top muted to an angle of about 45 and a 20 db/s gain was applied. Every other one of the resulting traces are shown in FIG. 5A. Incidence angles vwere computed and regression was carried out. Average gradient B.sub.avg was destretched. Components of the estimate vector {right arrow over (X)} are shown as first four traces in FIG. 5B. Spin-up gradient and azimuth traces were then computed using the pointwise method.

(29) Comparing these results with theoretical ones computed from equation 4 and listed in Table 1, there is a qualitative (but not quantitative) agreement of the average and spin-up gradient estimates. The lack of quantitative agreement may be due to computation of spherical divergence and incidence angle, which assumes isotropic media. Estimates of spin-up and maximum-amplitude azimuths were very unstable, taking on virtually every possible value, depending precisely on where they were sampled. This instability may be due to slight time misalignments of the seismic data at different source-receiver azimuths and/or slightly amount of noise that is present in the synthetic data.

(30) Windowed Statistical Method

(31) The windowed statistical method of performing AVAZ analysis can obtain estimates of symmetry azimuth .sub.sym that is more stable in the presence of noise and small misalignments of the data from one azimuth to the next as well as resolve the ambiguity problem. For a single interface with ideal data, a crossplot of these B.sub.sin versus B.sub.cos estimates would form a straight line, as the band-limited B.sub.ani varies over positive and negative values. The slope of this line is affected by the symmetry azimuth, as shown in FIGS. 6A-6H. FIG. 6A shows for .sub.sym=22.5. FIG. 6B shows for .sub.sym=0. FIG. 6C shows for .sub.sym=22.5. FIG. 6D shows for .sub.sym=45. FIG. 6E shows for .sub.sym=67.5. FIG. 6F shows for .sub.sym=90. FIG. 6G shows for .sub.sym=67.5. FIG. 6H shows for .sub.sym=45. For example, if .sub.sym=0, equation (6) shows that B.sub.sin=0, and the line would have zero slope. If .sub.sym=45, then B.sub.cos=0 and the line is vertical. If .sub.sym=22.5, the B.sub.cos and B.sub.sin attributes are equal, and will plot along a 45 line. However, these crossplot slopes are not unique. A symmetry azimuth of 67.5 would also have a 45 slope.

(32) More realistically, a crossplot of B.sub.sin versus B.sub.cos over an arbitrary window of data would produce an approximate elliptical distribution of points, similar to ellipse shown in FIG. 7A. FIG. 7B shows plot of the down interpretation. FIG. 7C shows plot of R.sub.r vs M cos 2. FIG. 7D shows plot of the up interpretation. This window could form a rectangular prism centered about an analysis point, about which joint statistical parameters of the data would be attributed. Alternatively, the window would follow interpreted geological markers, as shown in FIG. 8. A weighting function could be constructed within the window, which gradually reduces emphasis placed on data samples near edge of the window.

(33) An input vector can be formed from B.sub.cos and B.sub.sin input traces:

(34) $\begin{array}{cc}\stackrel{.fwdarw.}{s}\left(t\right)=\left[\begin{array}{c}a\left(t\right)\\ b\left(t\right)\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{co}s}\left(t\right)\\ B_{sin}\left(t\right)\end{array}\right]& \left(11\right)\end{array}$
The joint correlation matrix of {right arrow over (s)} is the expected value of {right arrow over (s)}{right arrow over (s)}.sup.T averaged over a sliding time-space window. This correlation is the 22 matrix .sub.s, which is computed as

(35) $\begin{array}{cc}{}_{\stackrel{.fwdarw.}{s}}=E\left\{\stackrel{.fwdarw.}{s}{\stackrel{.fwdarw.}{s}}^T\right\}=\left(\begin{array}{cc}{}_a^2& R\\ R& {}_b^2\end{array}\right)& \left(12\right)\end{array}$
whose elements are slowly-varying functions of time and space. They are computed by convolving the data in a super-gather by the weighting functions described.

(36) A unitary matrix, U.sup., can be found which transforms input signal vector {right arrow over (s)} into a new predicted anisotropy signal
{right arrow over (s)}=U.sup.{right arrow over (s)}=[ab].sup.T(13)
whose components a and b are uncorrelated with each other. The stronger prediction component, a, referred to as the background while the weaker component, b, is referred to as the anomaly. The correlation matrix .sub.{right arrow over (s)} of the predicted signal is given by

(37) $\begin{array}{cc}\begin{array}{c}{}_{\stackrel{.fwdarw.}{s^{}}}=E\left\{s^{}{s}^{}\right\}\\ =E\left\{U^{}s{s}^{}U\right\}\\ =U^{}{}_{\stackrel{.fwdarw.}{s}}U\\ =U^{-1}{}_{\stackrel{.fwdarw.}{s}}U\\ =\left(\begin{array}{cc}{}_a^2& R^{}\\ R^{}& {}_b^2\end{array}\right)\end{array}& \left(14\right)\end{array}$

(38) It is required that the unitary transform U.sup. make .sub.a.sub.b and R=0.

(39) Real Solutions

(40) If a(t) and b(t) are real seismic traces, with no quadrature (imaginary) component, joint statistics can be computed as follows:

(41) $\begin{array}{cc}{}_a^2=1/N\underset{j=1}{\overset{N}{.Math.}}{a}_j^2& \left(15\right)\\ {}_b^2=1/N\underset{j=1}{\overset{N}{.Math.}}{b}_j^2& \left(16\right)\\ R=1/N\underset{j=1}{\overset{N}{.Math.}}{a}_j{b}_j& \left(17\right)\end{array}$

(42) where the sums are over a sliding window. The statistics P is defined to be the average power in a and b, and M to be half the difference in powers between a and b:
P=(.sub.a.sup.2+.sub.b.sup.2)(18)
M=(.sub.a.sup.2.sub.b.sup.2)(19)

(43) A coordinate rotation by an angle of can be obtained by pre-multiplying the signal vector {right arrow over (s)} by the unitary matrix

(44) $\begin{array}{cc}{U}^{}=\left(\begin{array}{cc}\cos & -\sin \\ \sin & \cos \end{array}\right)& \left(20\right)\end{array}$
where 2 is the angle (M, R) vector makes with the Maxis. As shown in FIG. 8 and in FIG. 3b:
2=a tan 2[R,M]2n(21)
where a tan 2 [y, x] is the clockwise angle from the x-axis to the vector (y, x) and where n is any integer. If n is an even number, is said to be a spin-up solution. If it is odd, it is said to be a spin-down solution. Since all solutions having the same parity are equivalent, n=0 can be used for the spin-up solution and n=sign (.sub.up) for the spin-down solution. These choices ensure that both up and down rotation angles will always be between 180 as shown in FIGS. 7A-7D. More generally, it is always possible to choose a rotation angle in the range of a range of .sub.p90, where .sub.p is an arbitrary preferred azimuth opposite the branch line shown in FIGS. 7A-7D.
Complex Solutions

(45) The real solution can suffer from several limitations. Slight misalignments in the B.sub.cos and B.sub.sin traces can lead to large variations in the computed azimuths. Due to the narrow-band nature of seismic data, these time differences can manifest themselves as approximate phase difference between these traces. To overcome this difficulty, equation (11) can be generalized to:

(46) $\begin{array}{cc}\stackrel{.fwdarw.}{s}\left(t\right)=\left[\begin{array}{c}a\left(t\right)\\ b\left(t\right)\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{co}s}\left(t\right)\\ B_{sin}\left(t\right)\end{array}\right]+i\left[\begin{array}{c}\mathrm{Hi}\left\{B_{\mathrm{co}s}\left(t\right)\right\}\\ \mathrm{Hi}\left\{B_{\mathrm{si}n}\left(t\right)\right\}\end{array}\right]& \left(22\right)\end{array}$
Where Hi{ } denotes the 90 phase-rotated signal (also known as the Hilbert transform). The complex correlation matrix

(47) $\begin{array}{cc}{}_{\stackrel{.fwdarw.}{s}}=E\left\{\stackrel{.fwdarw.}{s}{\stackrel{.fwdarw.}{s}}^{}\right\}=\left(\begin{array}{cc}{}_a^2& R\\ R^{*}& {}_b^2\end{array}\right)& \left(23\right)\\ \mathrm{where}& \\ {}_a^2=1/N\underset{j=1}{\overset{N}{.Math.}}{.Math.a_j.Math.}^2;{}_b^2=1/N\underset{j=1}{\overset{N}{.Math.}}{.Math.b_j.Math.}^2;R=1/N\underset{j=1}{\overset{N}{.Math.}}{a}_j{b}_j^{*}=R_r+{\mathrm{iR}}_i& \left(24\right)\end{array}$
and where a.sub.j and b.sub.j are samples of the continuous functions a(t) and b(t). In the complex case, the matrix U.sup., which diagonalizes .sub.{right arrow over (S)} and makes R=0 in equation (14), is not unique. However, a preferred embodiment chooses a unitary matrix of the form

(48) $\begin{array}{cc}{U}^{}=\left(\begin{array}{cc}{e}^{-i{}_{\mathrm{co}s}}& -e^{-i{}_{\mathrm{si}n}}\\ e^{-i{}_{\mathrm{si}n}}& e^{-i{}_{\mathrm{co}s}}\end{array}\right)& \left(25\right)\end{array}$
Note that it is a function of two angles, and . The latter represents half the phase difference between B.sub.cos and B.sub.sin. To diagonalize .sub.{right arrow over (S)}
2=a tan 2[R.sub.isign(R.sub.r),|R.sub.r|]2m(26)
2=a tan 2[|R|sign(R.sub.r),M]2n(27)
The result is an exact complex solution. If the exact anisotropic power .sub.e is denoted to be the length of the (M, R) vector, then eigenvalues {.sub.a, .sub.b} of the predicted correlation matrix .sub.{right arrow over (S)} (which will match those of the original correlation matrix .sub.{right arrow over (S)}) are P+ and P, respectively. The Cauchy-Schwarz inequality establishes that |R|.sup.2.sup.2.sub.a.sup.2.sub.b, so
.sub.e={square root over (|R|.sup.2+M.sup.2)}{square root over (.sub.a.sup.2.sub.b.sup.2+M.sup.2)}=P(28)
If the coefficient of anisotropy C.sub.a is denoted to be .sub.e/P, then value of C.sub.a lies between 0 and 1. The magnitude of anisotropy B.sub.mag is {square root over (2.sub.e)}.

(49) For m=0, value is restricted to lie within range 45. This means that actual phase difference that is tolerated is 90. If actual phase difference gradually exceeds this, it will result in a sudden discontinuity in the estimate of . This is a branch line in the direction. Since such discontinuities are generally undesirable, an alternative solution can be expressed as
2=sin.sup.1(R.sub.i/P)(29)
2=a tan 2[R.sub.r,M cos 2]2n(30)

(50) These solutions yield identical results as the exact solution given noise-free data, but gradually drive the estimate of the phase difference to zero as the noise level increases. Phase differences between the a and b traces approaching 45 masquerade as isotropic zones. Without determining spin boundaries in the phase direction, it is difficult or impossible to eliminate discontinuities in the a and b gradients. For this reason, a modified definition of anisotropic power may be preferable:
.sub.r={square root over (R.sub.r.sup.2+M.sup.2 cos.sup.2 2)}={square root over (R.sub.r.sup.2+M.sub.r.sup.2)}(31)
where M.sub.r=M cos 2 is the robust power difference in a and b, taking into account possible phase differences between them. It becomes clear that .sub.r.sub.eP. The spin-up azimuth .sub.up is found to be .sub.up and the spin-down azimuth .sub.down is .sub.down=.sub.upsign(.sub.up).

(51) FIG. 10A shows the results of this analysis on the synthetic gather of FIG. 5. The first two traces, B.sub.cos and B.sub.sin, represent azimuthal gradient in Cartesian coordinates aligned along the seismic survey axes. They come directly from the least-squares analysis. Next two traces, B.sub.up and B.sub.anom, represent the azimuthal gradient in Cartesian coordinates along and perpendicular to the spin-up symmetry axis (.sub.up). The B.sub.anom has a smaller rms magnitude than B.sub.up. The next trace (.sub.r) is the robust magnitude of anisotropy. The trace shows the estimated phase differences between B.sub.cos and B.sub.sin. This can range in values from 32 to 23.

(52) Referring to FIG. 10A, the last trace shows the estimated spin-up azimuth (.sub.up) obtained from .sub.up estimated with a 30 ms time window. It is constrained to lie within 45 with an azimuth of 45 being equivalent to an angle of +45. Events 3 and 4 are supposed to have a symmetry azimuth of 45, which is the branch angle that transitions between spin-up and spin-down. These spin transitions are visible as sharp 90 azimuth jumps at 544 and 594 ms. B.sub.up and B.sub.anom would reveal abrupt polarity flips at these times although their amplitudes at these times are low. These sudden flips in azimuth and polarity are symptomatic of inherent ambiguity of AVAZ analysis, which are further alleviated by regions analysis, described next. FIG. 10B shows the synthetic results using windowed statistical method after regions analysis.

(53) Regions Analysis

(54) The rotation angle and the symmetry azimuth a are obtained from the windowed statistic method by computing the angle of the {right arrow over (S)}=(M.sub.r, R.sub.r) statistics vector shown in the upper right portion FIG. 7A. If length .sub.r of this vector is sufficiently large, then an incremental change in either R.sub.r or M.sub.r or both will lead to an incremental change in the rotation angle, except at the branch line (shown in thick bold). Since the symmetry azimuth is half the rotation angle, the angle 2 of the statistic vector should be determined to within 360 by tracking how many times the vector rotates completely around the circle without plunging into the incoherent zone near the origin. This can be done by subdividing the (M.sub.r, R.sub.r) plane into a circular incoherent zone (Zone 0) centered around the origin and 2N.sub.z coherent zones outside this circle as shown in FIG. 11A. Radius of this circle can be chosen to be kP.sub.s where k is a positive number less than one (minimum coefficient of anisotropy) and P.sub.s is a stabilized average power in the B.sub.cos and B.sub.sin traces averaged over the correlation window:
P.sub.S=(.sub.a.sup.2,.sub.b.sup.2)+[other term](32)
The 3-D seismic volume can then be subdivided into a set of contiguous interconnected regions, each of which resides within a single zone in the (M.sub.r, R.sub.r) plane (FIGS. 11B and 11C). Other terms may be included in equation (32) to further subdivide coherent regions into smaller sub-regions to make them less interconnected with one another. This will afford greater flexibility in assigning spins.

(55) Valid boundaries between regions are those that do not skip across Zone 0. These would include, for example, a boundary between Zones 1 and 1 but not between Zones 1 and 2. If too many boundaries are not valid, the (M.sub.r, R.sub.r) plane may have to be subdivided again using a larger radius for the incoherent zone. Chains of coherent regions that are connected through valid boundaries are then assembled and recursively inspected for parity violations. A parity violation occurs if any circular pathway through the 3-D volume passes through an odd number of valid coherent boundaries. For each parity violation, one of the boundaries can be designated as weakest. The weakest boundary could be one having smallest surface area or smallest average coherence. Conflicts may be resolved or avoided by marking the smaller of the two regions across the weakest boundary as being incoherent. If this eliminates too much of the 3-D volume, an alternative is to eliminate the weakest boundary although this can cause that boundary to retain its abrupt change in polarity of B.sub.back and B.sub.anon and cause the estimation of symmetry azimuth to suddenly shift by 90.

(56) Starting from an arbitrary region in the chain, initial spins are consistently assigned to the coherent regions contained in the chain so that a reversal in spin occurs only when a valid spin-flip boundary is crossed from either direction. Spin-flip boundaries are shown as dashed lines in FIGS. 11A-11C. They represent a transition between coherent zones N.sub.z and N.sub.z. Initial spins can then be converted into final spins by reversing all the spins in each coherent chain according to the following hierarchy: 1. If one or more pseudo-wells (actual or hypothetical) with known spins penetrate the regions of coherent chain and if a majority of these spins contradict the initial spins, then all of the spins of the chain should be reversed. If the majority of the well spins agree with the initial spins, then the final spins are set to match the initial spins. 2. If the final spins of a chain were not determined by pseudo-wells but one or more seismic horizons with known spins or B.sub.back polarities are known or conjectured, then if a majority of these spins or polarities conflict with those that were initially assigned then all of the spins of the chain should be reversed. On the other hand, if a majority of these spins or polarities match those that were initially assigned then the final spins should be set equal to the initial spins. 3. If the final spins of a chain were determined by neither pseudo-wells nor horizons then the final spins of the chain are either set to match or oppose the initial spins so as to make the majority of spins (based on volume) to match a predetermined regional spin.

(57) FIG. 9 illustrates that joint correlations can be calculated within variable-height windows that conform in depth to geologic horizons [1, 2] and interpolated horizons [3]. Each window have fixed lateral widths, fixed vertical height [6] at window's analysis point [4, 5] about which the window is centered. FIG. 10B shows the results of this analysis on the synthetic gather of FIG. 5A. Only two coherent zones were used for this example, corresponding to FIG. 10A (N.sub.z==1). The traces labeled BPI is the Branch Proximity Indicator given by

(58) $\begin{array}{cc}\mathrm{BPI}=\frac{\sin 2\left({}_{\mathrm{up}}-{}_p\right)}{{\sin }^{2}X+{\cos }^2\left({}_{\mathrm{up}}-{}_p\right)}& \left(33\right)\end{array}$
In this case, the preferred azimuth .sub.p=0 and the width of the indicator was X=1. This indicator is used to determine when branch line opposite the preferred azimuth (180 in this case) is crossed. It will suddenly switch polarity when this happens (i.e., at spin-reversal points 544 and 594 ms). The trace labeled CA represents the Coefficient of Anisotropy. This is always a positive number between 0 and 1, which denotes the normalized anisotropy magnitude, CA=.sub.r/P.sub.s. Additional term included in equation (32) was 5% of the average AVO gradient B.sub.avg smoothed over a 1000 ms window. This was done so as to cause the CA to fall to zero in the no-data zones between events. Minimum threshold of anisotropy (k) was set to 0.1 which caused each reflector to be placed into a separate coherent region. The trace labeled .sub.sym is the estimate of the symmetry azimuth. Once the azimuth is estimated, it can be used as a component of the azimuthal gradient in its direction (B.sub.ani) and in its perpendicular direction (B.sub.anom).

(59) In comparing the right-most traces of FIGS. 10A and 10B, it should be apparent that the windowed statistical method followed by regions analysis produced a more stable estimate of fracture symmetry azimuth than the convention point-by-point method. Furthermore, the results of the former method were accurate to within 2% of their correct values. Zones which cannot be computed due to lack of data are clearly marked in black.

(60) In closing, it should be noted that the discussion of any reference is not an admission that it is prior art to the present invention, especially any reference that may have a publication date after the priority date of this application. At the same time, each and every claim below is hereby incorporated into this detailed description or specification as additional embodiments of the present invention.

(61) Although the systems and processes described herein have been described in detail, it should be understood that various changes, substitutions, and alterations can be made without departing from the spirit and scope of the invention as defined by the following claims. Those skilled in the art may be able to study the preferred embodiments and identify other ways to practice the invention that are not exactly as described herein. It is the intent of the inventors that variations and equivalents of the invention are within the scope of the claims, while the description, abstract and drawings are not to be used to limit the scope of the invention. The invention is specifically intended to be as broad as the claims below and their equivalents.

REFERENCES

(62) All of the references cited herein are expressly incorporated by reference. The discussion of any reference is not an admission that it is prior art to the present invention, especially any reference that may have a publication data after the priority date of this application. Incorporated references are listed again here for convenience:

(63) 1. Al-Shuhail, A. A., 2007, Fracture-porosity inversion from P-wave AVOA data along 2D seismic lines: An example from the Austin chalk of southeast Texas: Geophysics, 72, no. 1, B1-B7.

(64) 2. Davidson, M., Swan, H., Sil, S., Howell, J., Olson, R., and Zhou, C., 2011, A robust workflow for detecting azimuthal anisotropy: SEG, Expanded Abstracts, 30, no. 1, 259-263.

(65) 3. Davison, C., Ratcliffe, A., Grion, S., Johnston, R., Duque, C., and Maharramov, M., 2011, Azimuthal AVO analysis: Stabilizing themodel parameters: SEG, Expanded Abstracts, 30, no. 1, 330-334.

(66) 4. Downton, J., and Gray, D., 2006, AVAZ parameter uncertainty estimation: SEG, Expanded Abstracts, 25, no. 1, 234-238.

(67) 5. Downton, J., Roure, B., and Hunt, L., 2011, Azimuthal Fourier coefficients: CSEG Recorder, 36, no. 10, 22-36.

(68) 6. Gray, D., and Head, K., 2000a, Fracture detection in Manderson field: A 3-D AVAZ case history: The Leading Edge, 19, no. 11, 1214-1221.

(69) 7. Gray, D., Roberts, G., and Head, K., 2002, Recent advances in determination of fracture strike and crack density from P-wave seismic data: The Leading Edge, 21, no. 3, 280-285.

(70) 8. Gray, F. D., Zellou, A., Todorovic-Marinic, D., and Boerner, S., 2005, Reservoir fracture characterization: U.S. Pat. No. 6,928,367.

(71) 9. Jenner, E., and Williams, M. C., 2004, System for estimating azimuthal variations in seismic data: U.S. Pat. No. 6,681,184.

(72) 10. Jenner, E., 2002, Azimuthal AVO: Methodology and data examples: The Leading Edge, 21, no. 8, 782-786.

(73) 11. Mallick, S., and Frazer, L. N., 1987, Practical aspects of reflectivity modeling: Geophysics, 52, no. 10, 1355-1364.

(74) 12. Mallick, S., Chambers, R. E., and Gonzalez, A., 1996, Method for determining the principal axes of azimuthal anisotropy from P-wave data: U.S. Pat. No. 5,508,973.

(75) 13. Mallick, S., Craft, K. L., Meister, L. J., and Chambers, R. E., 1998, Determination of the principal directions of azimuthal anisotropy from P-wave seismic data: Geophysics, 63, no. 2, 692-706.

(76) 14. Rger, A., and Tsvankin, I., 1997, Using AVO for fracture detection: Analytic basis and practical solutions: The Leading Edge, 16, no. 10, 1429-1434.

(77) 15. Rger, A., 1995, P-wave reflection coefficients for transversely isotropic media with vertical and horizontal axis of symmetry: SEG, Expanded Abstracts, 14, 278-281.

(78) 16. Rger, A., 1997, Plane wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry: Geophysics, 62, no. 3, 713-722.

(79) 17. Rger, A., 2000, Variation of P-wave reflectivity with offset and azimuth in anisotropic media: SEG Books, 20, 277-289.

(80) 18. Swan, H. W., 1990, Amplitude versus offset measurement errors in a finely layered medium: Geophysics, 56, no. 1, 41-49.

(81) 19. Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, no. 10, 1954-1966.

(82) 20. Tsvankin, I., 1996, Pwave signatures and notation for transversely isotropic media: An overview: Geophysics, 61, no. 2, 467-483.

(83) 21. Wang, J., Zheng, Y., and Perz, M., 2007, VVAZ vs. AVAZ: Practical implementation and comparison of two fracture-detection methods: SEG, Expanded Abstracts, 26, no. 1, 189-193.

(84) 22. Whitcombe, D. N., Dyce, M., McKenzie, C. J. S., and Hoeber, H., 2004, Stabilizing the AVO gradient: SEG, Expanded Abstracts, 23, no. 1, 232-235.

(85) 23. Xu, X., and Tsvankin, I., 2006, Azimuthal AVO analysis with anisotropic spreading correction: A synthetic study: The Leading Edge, 25, no. 11, 1336-1342.

(86) 24. Zheng, Y., and Larson, G., 2004, Seismic fracture detection: Ambiguity and practical solution: SEG, Expanded Abstracts, 23, no. 1, 1575-1578.