Polynomial amplitude versus azimuth inversion in horizontally transverse isotropic media, as tested on fractured coal seams in the Surat Basin*

Keywords: anisotropy, AVAZ, AVO, cracks, fractures, seismic, WAZ.

Amplitude versus azimuth (AVAZ) inversion in horizontally transverse isotropic (HTI) media is an effective method of characterising anisotropic variation within individual reflectors, as well as characterising fractures. Rüger (1998) introduced the reflectivity equation in HTI media assuming weak anisotropy. Several AVAZ methods of characterising azimuthal variation in seismic exist, such as assuming ellipticity (Rüger 1998) and using azimuthal Fourier transformations (Downton and Roure 2015).

In this paper, the P-wave reflectivity equation in HTI media is reformulated to allow for the simple separation of the isotropic component from the anisotropic component of the reflectivity by recognising that the P-wave reflectivity varies parabolically with respect to the squared cosine of the azimuth, cos² φ. Using this reformulation, the anisotropic component is then separated into the substituent elliptic and anelliptic anisotropic components that are extended to their respective AVO gradient terms. The seismic response of wet and dry crack densities is then defined (Bakulin et al. 2000). A method of AVO inversion to calculate individual Thomsen parameters is then proposed and error relationships between different techniques of characterising azimuthal anisotropy are approximated.

The inversion procedure presented herein is applied to a seismic survey from the Surat Basin and the inverted anisotropy is shown to match the anisotropy that was calculated from wire line logs to within 5%.

The P-wave reflectivity in HTI media, R_P^HTI assuming weak anisotropy is

where Z is the impedance, G is the shear modulus, β is the S-wave velocity, α is the P-wave velocity, i is the incident angle, φ is the azimuth and γ, δ and ε are Thomsen (1986) parameters (Rüger 1998). Furthermore a line over a parameter indicates the average of a variable about the reflecting interface and the superscript (v) indicates that the parameter occurs in the verticle plane. Eqn 1 is in a form where the isotropic component, R_P^iso(i) has been separated from the anisotropic component, R_P^ani(i, φ) of the reflectivity such that

The anisotropy has elliptic and anelliptic components (Alkhalifah and Tsvankin 1995) that can now be decoupled from the anisotropic component R_P^ani(i, φ). For ease, the ellipticity, E(i), and anellipticity, F(i), parameters are introduced and are defined as

The P-wave reflectivity in HTI media in terms of cos² φ can now be characterised by the polynomial

which is parabolic about cos² φ. Eqn 5 is in a form where pre-stack azimuthal seismic data can be inverted to find the isotropic, R_P^iso(i), elliptic E(i) and anelliptic, F(i) components of the azimuthal reflectivity provided the azimuth φ is known. The azimuth φ is defined as φ = φ_az – φ_sym (Rüger 1998) where φ_az is an arbitrary azimuth and φ_sym is the symmetry azimuth of the anisotropy. This can be calculated using the elliptic approximation (Rüger 1998) or by using a Fourier transform (Downton and Roure 2015).

To extract the parameters E(i), F(i) and R_P^iso(i) from seismic reflection data, the effect of the wavelet, w(t), must be removed. For zero phase data, the seismic reflection amplitude, A(t_m), recorded at a time, t_m is

Hence, the reflectivity over a horizon at t_m for a zero-phase wavelet is

and for an arbitrary wavelet is

where is a Fourier transform operator about time, t. Equations 7 and 8 can be extended to calculate R_P^iso(i), E(i) and F(i).

To account for the total amount of azimuthal anisotropy, the anisotropic response over all azimuths can be calculated by integrating over the domain Φ where φ∈Φ = [0, π), such that

The ratio of the anisotropic to isotropic reflectivity is therefore

and will be referred to as the anisotropy ratio. The use of this term does not require the reflectivity to be decoupled from the wavelet and can therefore be used to tie inverted seismic to measurements taken from wire line logs.

Using the AVO approximation (Shuey 1985), three gradients can be defined from Eqn 5

where B^iso and B^ani are the isotropic and anisotropic AVO gradients in HTI media, as introduced by Rüger (1998), respectively. Thus, the azimuthal AVO gradient equation in HTI media is

Parameters from Thomsen (1986), namely Δδ^(V), Δε^(V) and Δγ can be calculated by linearising B₁^ani and B₂^ani about tan² i. B₁^ani linearised about tan² i is defined by the gradient term

and the intercept term

which allows for the calculation of Δδ^(V). Δγ can be calculated using

where and are obtained from well log data, or perhaps from an isotropic AVO inversion if the bed thickness and data quality were adequate. Δε^(V) can be calculated via

Bakulin et al. (2000) applied established fracture theory to B^ani assuming weak anisotropy. This theory is applied to the parameters B₁^ani and B₂^ani assuming HTI media due to small penny-shaped cracks (Hudson 1980) which are overlain by an isotropic layer. In the case of wet cracks, the AVO gradients are approximated to

where and e is the crack density, which is defined as the number of cracks, N, per unit volume, V, multiplied by the average radius, r, of the penny-shaped cracks () (Molotkov and Bakulin 1997; Tsvankin and Grechka 2011).

In the case of dry cracks, the AVO gradients are approximated to

and could be used to distinguish between gas and liquid-filled cracks (Bakulin et al. 2000).

Under the elliptic assumption, the P-wave reflectivity of HTI media is simplified to

The long and short axes values of an inverted ellipse are assumed to be approximately the same as R_P^HTI(i, 0) and R_P^ani(i, 0) respectively (Rüger 1998). Thus, the approximate error associated with the elliptic approach is a function of the anelliptic component of the reflectivity such that

The near offset AVO gradient equation in HTI media is

and assumes that the reflectivity of an HTI interface varies elliptically with azimuth, φ, where (Rüger 1998, 2001). Hence, the error associated with this approximation is

The Fourier method fully models P-wave reflection in HTI media as described by Rüger (1998). If the incident angle, i, is known, then there are four independent terms in Eqn 5. However, when Eqn 1 undergoes a Fourier transform, there are six independent terms (Pšenčík and Martins 2001; Downton and Roure 2015)

Because of this, a minimum of six azimuthal amplitude measurements are needed to define Eqn 23 and, this is therefore a more complex approach than the method presented herein.

A high-order Fourier coefficient term, r₄, can be calculated from Eqn 23 (Downton and Roure 2015) that corresponds to an anelliptic term, Δη^(V) /16 where (Alkhalifah and Tsvankin 1995; Downton and Roure 2015). However, if weak anisotropy is assumed, as is the case of the parametrisation of Eqn 23 (Pšenčík and Martins 2001), then the anelliptic term in weak anisotropy is (Bakulin et al. 2000). This is also evident as the anelliptic component of Eqn 1 is . Because ε^(V) does not occur in the elliptic or isotropic components of Rüger’s equation, the Fourier method does not easily separate the elliptic component of the anisotropy from the anelliptic component because it is a component of all Fourier coefficients the terms r₀, r₂ and r₄ which are introduced in (Downton and Roure 2015). In contrast, E and B₁^ani exactly represent elliptic components of the reflectivity and AVO gradient in HTI media respectively, making the methodology presented herein a robust technique of calculating the elliptic and anelliptic anisotropic components.

The r₂ coefficient (which is introduced by Downton and Roure (2015)) approximates the AVO gradient, B^ani, whereby, as reported by Downton and Roure (2015)

The error associated with this approximation with regard to the true elliptic component of the AVO gradient, B₁^ani, is of the form

This term should not be used to approximate the elliptic gradient B^ani because the associated error as can be reduced from Eqn 22 is

which is the magnitude of the anelliptic term and therefore renders Eqn 24 ineffective.

Data from an open file 3D seismic survey in the Surat Basin were inverted to demonstrate the effectiveness of these techniques. The equations were applied by extracting seismic amplitude measurements from six azimuth volumes along two horizons in the survey, centred around a well in which density and full-wave sonic logs were acquired. The azimuth volumes were stacked using a restricted reflection angle of incidence of 37.5°.

To compare and assess the effectiveness of the equations, Thomsen’s parameters in HTI media were also calculated from the well logs using the method described by Kremor and Amrouch (2017). A comparison of the results for the well log measurements and the equations in this paper is given in Table 1. For both horizons, the error between the two anisotropy estimates is less than 5%.

Table 1.  Anisotropy ratio calculation
The parameters obtained from wire line logs are shown below along with the anisotropy ratio that was calculated using these parameters. The anisotropy ratio, calculated from a seismic survey in the Surat Basin, is shown, as well as the error between this and the result obtained from the wire line logs

The authors recognise that the quality of the results are highly dependent upon the preservation of signal amplitudes of the seismic data during processing. Although these appear to have been adequate for these data at this well location, it is likely that the wider application of this technique will require considerable cooperation with processing houses and appropriate selection of seismic horizons that are free of multiple energy and other noise.

A method of categorising seismic anisotropy caused by HTI reflectors that can easily be separated into its isotropic, elliptic and anelliptic components is presented. This was then extended to AVO and the calculation of Thomsen (1986) parameters about an interface that were then related to penny-shaped cracks. The ratio between the anisotropic and isotropic components of the reflectivity calculated using this method has been tied to a well to within 5% error. This method affords the geoscientist more interpretable attributes than what currently available AVAZ methods for HTI media can produce and will therefore be useful in fractured and unconventional reservoirs. Furthermore, the application of this method over the elliptic method will result in a reduced amount of error incorporated into the inverted parameters and is a simpler approach than the Fourier method.