compute AVO attributes from seismic (v2.1)¶

Originally thought as a simple add-on to the AVO projection notebook to show how to compute AVO attributes from seismic data. It now forms the basis of a more complete document showing how to do AVO analysis on seismic data using in Python.

Revision history:

• 2018-04: changed seg-y library from obspy to segyio
• 2018-02: complete revision of the notebook. Added 2 more functions to calculate AVO Intercept and Gradient; analysis on actual 2D and 3D seismic data and not just maps, compare the results of different methods of calculation, AVO crossplots and rotations.
• 2016-03: updated to Python 3.x
• 2015-12: first version

import libraries¶

The main libraries to import are the usual numpy and matplotlib, plus bruges, segyio to load SEG-Ys and xarray to store and manipulate 3D seismic cubes. On top of that, a few scipy methods to grid and interpolated data.

I use here the rainbow colormap from matplotlib even if it's one of the crappy ones just to make the final results look similar to the defaults of the so-called "professional" software.

load data¶

Also see the seismic data and amplitude extraction notebooks for a somewhat more comprehensive overview of using seismic data in Python.

calculate AVO attributes¶

I (and everybody else) use Shuey's two-term AVO equation which basically computes a linear regression of near and far amplitudes against $\sin^2(\theta)$, where $\theta$ is the average angle of incidence related to each angle stack; I don't actually have this piece of information for the QSI angle stacks but I can guesstimate to be around 10 and 30 degrees.

This is Shuey's equation:

$$R(\theta) = I + G \, \sin(\theta)^2$$

It's like the simple $y = mx + c$ equation, right? So $m$ is the slope of the line, i.e. our AVO gradient ($G$) and $c$ is the value of the line back at $y=0$ or the AVO intercept ($I$).

The space where I work has values of $\sin(\theta)^2$ along the horizontal axis and actual seismic amplitudes on the vertical axis. So if I have two angle stacks I can calculate AVO intercept and gradient for each sample of every trace in the volumes because for each sample I have:

• the amplitude of Near ($y_1$)
• the amplitude of Far ($y_2$)
• the value of $\sin(\theta)^2$ for Near ($x_1$)
• the value of $\sin(\theta)^2$ for Far ($x_2$)

And the last two coordinates will be fixed for all the samples, because $\theta$ is fixed for Near and Far (in my example, 10 and 30 degrees respectively).

So once I have the coordinates for two points, I can solve the above equation like this:

$$y_1 = m \cdot x_1 + c$$ $$y_2 = m \cdot x_2 + c$$ $$y_2 - y_1 = m \cdot (x_2 -x_1)$$ $$m = \frac{y_2 - y_1}{x_2 -x_1}$$

And that's the slope of the line, while the intercept is:

$$c = y_1 - m \cdot x_1$$

So if we talk amplitudes and angles these simple equations can be written as:

$$G = \frac{FAR - NEAR}{\sin(30)^2 - \sin(10)^2}$$ $$I = NEAR - G \, \sin(10)^2$$

This is a function that can be coded like this in Python:

This function only works with two inputs. If I have multiple inputs, i.e. if I want to calculate AVO attributes by using for example Near, Mid and Far angle stacks I need to either use scipy.stats.linregress or use np.linalg.lstsq to compute the least-squares solution to a linear matrix equation (the following function is modified from work by Wayne Mogg, author of a series of brilliant OpenDTect plugins):

** WORK IN PROGRESS // THE FOLLOWING CELLS ARE NOT YET UPDATED **

I will extract amplitudes from these two volumes along a certain horizon (it will be fairly easy to use the method also on timeslices, 2D sections and 3D volumes too).

The original dataset also has AVO Intercept and Gradient at Top Heimdal so I can compare Intercept and Gradient calculated in Python to see if they make sense.

The following lines are to extract amplitudes at Top Heimdal from the two cubes (see also this other notebook for techniques on amplitude extraction):

Then I grid Near (An) and Far (Af) amplitude maps:

To calculate intercept and gradient I will use simple linear regression on near and far amplitudes (i.e., Shuey's equation approach). I have tried numpy.polyfit and scipy.stats.linregress which have to be run inside a for loop, then a function which I have discovered on stackexchange which relies on arrays computation and makes it so fast it's almost unbelievable.

I need to reshape near and far amplitudes as well as create angle of incidence arrays (this actually took me far more time to figure it out then all the rest):

Holy cow: on my work PC numpy.polyfit does the job in 26 seconds, scipy.stats.linregress in 60 and multiple_linregress in 0.01 seconds!!! (On my personal computer, a Macbook Pro 13" with Intel Core i5 2.6Ghz, the times are 13, 14 and 0.01 seconds; the third solution is so fast that it's tempting to try and calculate Intercept and Gradient on the entire volume).

Let's see if there are any differences between np.polyfit and multiple_linregress. A quick look at the distribution of intercept and gradient calculated with the two methods shows similar shapes and ranges (in the plot below, blue is np.polyfit and red is multiple_linregress).

And if we look at the resulting maps we also see that they're the same (I also show the differences between the intercept and gradient maps calculated with the two methods):

display¶

Displaying the original AVO maps with the new ones I have calculated I see similar ranges of values but different pictures; the main problem is that I don't really know anything about the provided intercept and gradient maps, so I will need to do compare the results obtained in Python with the output of a commercial application (e.g., Petrel, Hampson-Russell, etc).

---

Thanks to Evan Bianco for comments and suggestions.