## Anisotropic AVO¶

Trying to reproduce some figures in Blangy, JP, 1994, AVO in tranversely isotropic media—An overview. Geophysics 59 (5), 775–781. This notebook goes with a blog post on Agile Geoscience.

The key equation is:

$$A = \frac12 \left[ \frac{\Delta\rho}{\rho} + \frac{\Delta V_\textrm{P}}{V_\textrm{P}} \right]$$$$B = 2 \frac{V_\textrm{S}^2}{V_\textrm{P}^2} \times \left[ \frac{\Delta\rho}{\rho} + 2 \frac{\Delta V_\textrm{S}}{V_\textrm{S}} \right] \sin^2\! \theta$$$$C = \frac12 \left( \frac{\Delta V_\textrm{P}}{V_\textrm{P}} \right) \tan^2\!\theta$$$$D = \frac12 \Delta\delta \sin^2\!\theta$$$$E = \frac12 (\Delta \delta - \Delta \epsilon ) \sin^2\!\theta \, \tan^2\!\theta$$$$A_\textrm{RP} = A - B + C + D - E$$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
# Test data from Table 1
type1 = {'shale':      {'vp':3300., 'vs':1700., 'rho':2350., 'd':0.15, 'e':0.30},
'sand_gas':   {'vp':4200., 'vs':2700., 'rho':2350., 'd':0.00, 'e':0.00},
'sand_water': {'vp':4200., 'vs':2100., 'rho':2450., 'd':0.00, 'e':0.00},
}
type2 = {'shale':      {'vp':2896., 'vs':1402., 'rho':2250., 'd':0.15, 'e':0.30},
'sand_gas':   {'vp':3322., 'vs':2215., 'rho':2000., 'd':0.00, 'e':0.00},
'sand_water': {'vp':3322., 'vs':1402., 'rho':2250., 'd':0.00, 'e':0.00},
}
type3 = {'shale':      {'vp':2307., 'vs':1108., 'rho':2150., 'd':0.15, 'e':0.30},
'sand_gas':   {'vp':1951., 'vs':1301., 'rho':1950., 'd':0.00, 'e':0.00},
'sand_water': {'vp':1951., 'vs': 930., 'rho':2200., 'd':0.00, 'e':0.00},
}

In [3]:
print type2['sand_gas']['d'] - type2['shale']['d'] # delta-delta
print type2['sand_gas']['e'] - type2['shale']['e'] # delta-epsilon

-0.15
-0.3

In [4]:
def blangy(interface, case='gas', i=0):

upper = interface['shale']
lower = interface['sand_'+case]

trans_angle = np.arcsin(np.sin(inc_angle) * lower['vp']/upper['vp'])
theta = 0.5 * (inc_angle + trans_angle)

vp  = (upper['vp'] + lower['vp'])/2.0
vs  = (upper['vs'] + lower['vs'])/2.0
rho = (upper['rho'] + lower['rho'])/2.0

dvp = lower['vp'] - upper['vp']
dvs = lower['vs'] - upper['vs']
drho = lower['rho'] - upper['rho']
dd = lower['d'] - upper['d']
de = lower['e'] - upper['e']

A = 0.5 * (drho/rho + dvp/vp)
B = 2.0 * (vs**2 / vp**2) * ((drho/rho + 2 * dvs/vs)) * np.sin(theta)**2
C = 0.5 * (dvp/vp) * np.tan(theta)**2
D = 0.5 * dd * np.sin(theta)**2
E = 0.5 * (dd - de) * np.sin(theta)**2 * np.tan(theta)**2

isotropic = A - B + C
anisotropic = A - B + C + D - E

return isotropic, anisotropic, (A, -B, C, D, -E)

In [5]:
r0 = blangy(type3, case='gas', i=0)
r0[1]

Out[5]:
-0.1323878151886263
In [6]:
r10 = blangy(type3, i=10)
r10[1]

Out[6]:
-0.1402341169325664

## Reproduce Figure 2¶

We will try to reproduce Figure 2 from the Blangy paper.

In [7]:
from IPython.display import Image
Image(filename='blangy_fig2_rearrange.png', width=800)

Out[7]: