# Source code for bruges.rockphysics.anisotropy

```"""
Anisotropy effects.

Backus anisotropy is from thin layers.

Hudson anisotropy is from crack defects.

"""
from collections import namedtuple
#from functools import wraps

import numpy as np

from bruges.rockphysics import moduli
from bruges.util import moving_average

[docs]def backus_parameters(vp, vs, rho, lb, dz):
"""
Intermediate parameters for Backus averaging. This is expected to be a
private function. You probably want backus() and not this.

Args:
vp (ndarray): P-wave interval velocity.
vs (ndarray): S-wave interval velocity.
rho (ndarray): Bulk density.
lb (float): The Backus averaging length in m.
dz (float): The depth sample interval in m.

Returns:
tuple: Liner's 5 intermediate parameters: A, C, F, L and M.

Notes:
Liner, C (2014), Long-wave elastic attenuation produced by horizontal
layering. The Leading Edge, June 2014, p 634-638.

"""
lam = moduli.lam(vp, vs, rho)
mu = moduli.mu(vp, vs, rho)

# Compute the layer parameters from Liner (2014) equation 2:
a = rho * np.power(vp, 2.0)  # Acoustic impedance

# Compute the Backus parameters from Liner (2014) equation 4:
A1 = 4 * moving_average(mu*(lam+mu)/a, lb/dz, mode='same')
A = A1 + np.power(moving_average(lam/a, lb/dz, mode='same'), 2.0)\
/ moving_average(1.0/a, lb/dz, mode='same')
C = 1.0 / moving_average(1.0/a, lb/dz, mode='same')
F = moving_average(lam/a, lb/dz, mode='same')\
/ moving_average(1.0/a, lb/dz, mode='same')
L = 1.0 / moving_average(1.0/mu, lb/dz, mode='same')
M = moving_average(mu, lb/dz, mode='same')

BackusResult = namedtuple('BackusResult', ['A', 'C', 'F', 'L', 'M'])
return BackusResult(A, C, F, L, M)

#
# def vectorize(func):
#     """
#     Decorator to make sure the inputs are arrays. We also add a dimension
#     to theta to make the functions work in an 'outer product' way.
#
#     Takes a reflectivity function requiring Vp, Vs, and RHOB for 2 rocks
#     (upper and lower), plus incidence angle theta, plus kwargs. Returns
#     that function with the arguments transformed to ndarrays.
#     """
#     @wraps(func)
#     def wrapper(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, **kwargs):
#         vp = np.asanyarray(vp, dtype=float)
#         vs = np.asanyarray(vs, dtype=float) + 1e-12  # Prevent singular matrix.
#         rho = np.asanyarray(rho, dtype=float)
#         lb = np.asanyarray(lb, dtype=float).reshape((-1, 1))
#         dz = np.asanyarray(dz, dtype=float)
#         return func(vp, vs, rho, lb, dz)
#     return wrapper

[docs]def backus(vp, vs, rho, lb, dz):
"""
Backus averaging. Using Liner's algorithm (2014; see Notes).

Args:
vp (ndarray): P-wave interval velocity.
vs (ndarray): S-wave interval velocity.
rho (ndarray): Bulk density.
lb (float): The Backus averaging length in m.
dz (float): The depth sample interval in m.

Returns:
namedtuple: the smoothed logs: vp, vs, plus rho. Useful for computing
other elastic parameters at a seismic scale.

Notes:
Liner, C (2014), Long-wave elastic attenuation produced by horizontal
layering. The Leading Edge, June 2014, p 634-638.

"""
# Compute the Backus parameters:
A, C, F, L, M = backus_parameters(vp, vs, rho, lb, dz)

# Compute the vertical velocities from Liner (2014) equation 5:
R = moving_average(rho, lb/dz, mode='same')
vp0 = np.sqrt(C / R)
vs0 = np.sqrt(L / R)

BackusResult = namedtuple('BackusResult', ['Vp', 'Vs', 'rho'])
return BackusResult(Vp=vp0, Vs=vs0, rho=R)

[docs]def backus_quality_factor(vp, vs, rho, lb, dz):
"""
Compute Qp and Qs from Liner (2014) equation 10.

Args:
vp (ndarray): P-wave interval velocity.
vs (ndarray): S-wave interval velocity.
rho (ndarray): Bulk density.
lb (float): The Backus averaging length in m.
dz (float): The depth sample interval in m.

Returns:
namedtuple: Qp and Qs.
"""
vp0, vs0, _ = backus(vp, vs, rho, lb, dz)

ptemp = np.pi * np.log(vp0 / vp) / (np.log(vp0 / vp) + np.log(lb/dz))
Qp = 1.0 / np.tan(ptemp)

stemp = np.pi * np.log(vs0 / vs) / (np.log(vs0 / vs) + np.log(lb/dz))
Qs = 1.0 / np.tan(stemp)

BackusResult = namedtuple('BackusResult', ['Qp', 'Qs'])
return BackusResult(Qp=Qp, Qs=Qs)

[docs]def thomsen_parameters(vp, vs, rho, lb, dz):
"""
Liner, C, and T Fei (2006). Layer-induced seismic anisotropy from
full-wave sonic logs: Theory, application, and validation.
Geophysics 71 (6), p D183–D190. DOI:10.1190/1.2356997

Args:
vp (ndarray): P-wave interval velocity.
vs (ndarray): S-wave interval velocity.
rho (ndarray): Bulk density.
lb (float): The Backus averaging length in m.
dz (float): The depth sample interval in m.

Returns:
namedtuple: delta, epsilon and gamma.

"""
A, C, F, L, M = backus_parameters(vp, vs, rho, lb, dz)

delta = ((F + L)**2.0 - (C - L)**2.0) / (2.0 * C * (C - L))
epsilon = (A - C) / (2.0 * C)
gamma = (M - L) / (2.0 * L)

ThomsenParameters = namedtuple('ThomsenParameters', ['δ', 'ε', 'γ'])
return ThomsenParameters(delta, epsilon, gamma)

[docs]def dispersion_parameter(qp):
"""
Kjartansson (1979). Journal of Geophysical Research, 84 (B9),
4737-4748. DOI: 10.1029/JB084iB09p04737.
"""
return np.arctan(1/qp) / np.pi

[docs]def blangy(vp1, vs1, rho1, d1, e1, vp0, vs0, rho0, d0, e0, theta):
"""
Blangy, JP, 1994, AVO in transversely isotropic media-An overview.
Geophysics 59 (5), 775-781. DOI: 10.1190/1.1443635

Provide Vp, Vs, rho, delta, epsilon for the upper and lower intervals,
and theta, the incidence angle.

:param vp1: The p-wave velocity of the upper medium.
:param vs1: The s-wave velocity of the upper medium.
:param rho1: The density of the upper medium.
:param d1: Thomsen's delta of the upper medium.
:param e1: Thomsen's epsilon of the upper medium.

:param vp0: The p-wave velocity of the lower medium.
:param vs0: The s-wave velocity of the lower medium.
:param rho0: The density of the lower medium.
:param d0: Thomsen's delta of the lower medium.
:param e0: Thomsen's epsilon of the lower medium.

:param theta: A scalar [degrees].

:returns: the isotropic and anisotropic reflectivities in a tuple. The
isotropic result is equivalent to Aki-Richards.

TODO
Use rocks.
"""
lower = {'vp': vp0,
'vs': vs0,
'rho': rho0,
'd': d0,
'e': e0,
}

upper = {'vp': vp1,
'vs': vs1,
'rho': rho1,
'd': d1,
'e': e1,
}

# Redefine theta
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 = isotropic + D - E

BlangyResult = namedtuple('BlangyResult', ['isotropic', 'anisotropic'])
return BlangyResult(isotropic, anisotropic)

[docs]def ruger(vp1, vs1, rho1, d1, e1, vp2, vs2, rho2, d2, e2, theta):
"""
Coded by Alessandro Amato del Monte and (c) 2016 by him

Rüger, A., 1997, P -wave reflection coefficients for transversely
isotropic models with vertical and horizontal axis of symmetry:
Geophysics, v. 62, no. 3, p. 713–722.

Provide Vp, Vs, rho, delta, epsilon for the upper and lower intervals,
and theta, the incidence angle.

:param vp1: The p-wave velocity of the upper medium.
:param vs1: The s-wave velocity of the upper medium.
:param rho1: The density of the upper medium.
:param d1: Thomsen's delta of the upper medium.
:param e1: Thomsen's epsilon of the upper medium.

:param vp0: The p-wave velocity of the lower medium.
:param vs0: The s-wave velocity of the lower medium.
:param rho0: The density of the lower medium.
:param d0: Thomsen's delta of the lower medium.
:param e0: Thomsen's epsilon of the lower medium.

:param theta: A scalar [degrees].

:returns: anisotropic reflectivity.

"""
vp = np.mean([vp1, vp2])
vs = np.mean([vs1, vs2])
z = np.mean([vp1*rho1, vp2*rho2])
g = np.mean([rho1*vs1**2, rho2*vs2**2])
dvp = vp2-vp1
z2, z1 = vp2*rho2, vp1*rho1
dz = z2-z1
dg = rho2*vs2**2 - rho1*vs1**2
dd = d2-d1
de = e2-e1
A = 0.5*(dz/z)
B = 0.5*(dvp/vp - (2*vs/vp)**2 * (dg/g) + dd) * np.sin(a)**2
C = 0.5*(dvp/vp + de) * np.sin(a)**2 * np.tan(a)**2
R = A+B+C

return R

[docs]def crack_density(porosity, aspect):
"""
Returns crack density from porosity and aspect ratio, phi and alpha
respectively in the unnumbered equation between 15.40 and 15.41 in
Dvorkin et al. 2014.

Args:
porosity (float): Fractional porosity.
aspect (float): Aspect ratio.

Returns:
float: Crack density.
"""
if porosity >= 1:
porosity /= 100.

return 3 * porosity / (4 * np.pi * aspect)

[docs]def hudson_delta_M(porosity, aspect, mu, lam=None, pmod=None):
"""
The approximate reduction in compressional modulus M in the direction
normal to a set of aligned cracks. Eqn 15.40 in Dvorkin et al (2014).

Args:
porosity (float): Fractional porosity, phi.
aspect (float): Aspect ratio, alpha.
mu (float): Shear modulus, sometimes called G.
lam (float): Lame's first parameter.
pmod (float): Compressional modulus, M.

Returns:
float: M_inf - M_0 = \Delta c_11.
"""
epsilon = crack_density(porosity, aspect)
if lam:
return epsilon * (lam**2 / mu) * 4*(lam + 2*mu)/(3*lam + 3*mu)
else:
return (4*epsilon/3) * ((pmod - 2*mu)**2 / mu) * (pmod/(pmod-mu))

[docs]def hudson_delta_G(porosity, aspect, mu, lam=None, pmod=None):
"""
The approximate reduction in shear modulus G (or mu) in the direction
normal to a set of aligned cracks. Eqn 15.42 in Dvorkin et al (2014).

Args:
porosity (float): Fractional porosity, phi.
aspect (float): Aspect ratio, alpha.
mu (float): Shear modulus, sometimes called G.
lam (float): Lame's first parameter, lambda.
pmod (float): Compressional modulus, M.

Returns:
float: M_inf - M_0 = \Delta c_11.
"""
epsilon = crack_density(porosity, aspect)
if lam:
return epsilon * mu * 16*(lam + 2*mu)/(9*lam + 12*mu)
else:
return (16*mu*epsilon/3) * pmod / (3*pmod - 2*mu)

[docs]def hudson_quality_factor(porosity, aspect, mu, lam=None, pmod=None):
"""
Returns Q_p and Q_s for cracked media. Equations 15.41 and 15.43 in
Dvorkin et al. (2014).

Args:
porosity (float): Fractional porosity, phi.
aspect (float): Aspect ratio, alpha.
mu (float): Shear modulus, sometimes called G.
lam (float): Lame's first parameter, lambda.
pmod (float): Compressional modulus, M.

Returns:
float: Q_p
float: Q_s
"""
Qp = 2*mu / hudson_delta_M(porosity, aspect, mu, lam, pmod)
Qs = 2*mu / hudson_delta_G(porosity, aspect, mu, lam, pmod)
return Qp, Qs

[docs]def hudson_inverse_Q_ratio(mu=None,
pmod=None,
pr=None,
vp=None,
vs=None,
aligned=True):
"""
Dvorkin et al. (2014), Eq 15.44 (aligned) and 15.48 (not aligned). You must
provide one of the following: `pr`, or `vp` and `vs`, or `mu` and `pmod`.

Args:
mu (float): Shear modulus, sometimes called G.
pmod (float): Compressional modulus, M.
pr (float): Poisson's ratio, somtimes called v.
vp (ndarray): P-wave interval velocity.
vs (ndarray): S-wave interval velocity.
aligned (bool): Either treats cracks as alligned (Default, True)
or assumes defects are randomly oriented (False)

Returns:
float: 2Q_s^-1
"""
if pr is not None:
x = (2 - 2*pr) / (1 - 2*pr)
elif (vp is not None) and (vs is not None):
x = vp**2 / vs**2
elif (mu is not None) and (pmod is not None):
x = pmod / mu
else:
raise TypeError("You must provide pr or (vp and vs) or (mu and pmod)")

if aligned:
return 0.25 * (x - 2)**2 * (3*x - 2) / (x**2 - x)
else:
a = 2*x / (3*x - 2)
b = x / 3*(x - 1)
return 1.25 * ((x - 2)**2 / (x - 1)) / (a + b)
```