Source code for bruges.reflection.reflection

"""
Various reflectivity algorithms.

:copyright: 2018 Agile Geoscience
:license: Apache 2.0
"""
from functools import wraps
from collections import namedtuple
import warnings

import numpy as np

from bruges.rockphysics import moduli
from bruges.rockphysics import anisotropy
from bruges.util import deprecated


[docs]def critical_angles(vp1, vp2, vs2=None): """Compute critical angle at an interface Given the upper (vp1) and lower (vp2) velocities at an interface. If you want the PS-wave critical angle as well, pass vs2 as well. Args: vp1 (ndarray): Upper layer P-wave velocity. vp2 (ndarray): Lower layer P-wave velocity. vs2 (ndarray): Lower layer S-wave velocity (optional). Returns: tuple: The first and second critical angles at the interface, in degrees. If there isn't a critical angle, it returns np.nan. """ ca1 = ca2 = np.nan if vp1 < vp2: ca1 = np.degrees(np.arcsin(vp1/vp2)) if (vs2 is not None) and (vp1 < vs2): ca2 = np.degrees(np.arcsin(vp1/vs2)) return ca1, ca2
[docs]def reflection_phase(reflectivity): """ Compute the phase of the reflectivity. Returns an array (or float) of the phase, in positive multiples of 180 deg or pi rad. So 1 is opposite phase. A value of 1.1 would be +1.1 \times \pi rad. Args: reflectivity (ndarray): The reflectivity, eg from `zoeppritz()`. Returns: ndarray: The phase, strictly positive """ ph = np.arctan2(np.imag(reflectivity), np.real(reflectivity)) / np.pi ph[ph == 1] = 0 ph[ph < 0] = 2 + ph[ph < 0] return ph
[docs]def acoustic_reflectivity(vp=None, rho=None, axis=0, mode='same'): """ The acoustic reflectivity, given Vp and RHOB logs. You can also pass either `vp` or `rho` and get velocity or density reflectivity; if you only have an impedance log or dataset, pass it as `vp`. For offset reflectivity, use `reflectivity()`. Args: vp (ndarray): The P-wave velocity. At least one of Vp and Rho is required. rho (ndarray): The bulk density. At least one of Vp and Rho is required. If you only have Impedance, pass it in as Vp and leave this as None. It will use a value of 1 to compute impedance and you'll get velocity reflectivity. axis (int): The dimension along which to compute the reflectivity. mode (str): 'same' means the output will be the same shape. 'valid' means that only strictly valid reflectivities are computed so the result will be one sample shorted in the `axis` dimension. Returns: ndarray: The reflectivity coefficient series. """ if (vp is None) and (rho is None): raise ValueError("You must pass either `vp` or `rho`.") elif vp is None: vp = np.ones_like(rho) elif rho is None: rho = np.ones_like(vp) if axis < 0: axis = vp.ndim + axis if axis > 0: imp = np.moveaxis(vp * rho, axis, 0) else: imp = vp * rho if mode == 'same': pad_width = [(0, 1)] + (imp.ndim - 1) * [(0, 0)] imp = np.pad(imp, pad_width=pad_width, mode='edge') rc = (imp[1:] - imp[:-1]) / (imp[1:] + imp[:-1]) if axis > 0: return np.moveaxis(rc, 0, axis) else: return rc
[docs]def reflectivity(vp, vs, rho, theta=0, method='zoeppritz_rpp', axis=0, mode='same'): """ Offset reflectivity, given Vp, Vs, rho, and offset. For acoustic reflectivity, use `acoustic_reflectivity()`. Computes 'upper' and 'lower' intervals from the three provided arrays, then passes the result to the specified method to compute reflection coefficients. For acoustic reflectivity, either use the `acoustic_reflectivity()` function, or call `reflectivity()` passing any log as Vs, e.g. just give the Vp log twice (it won't be used anyway): reflectivity(vp, vp, rho) For anisotropic reflectivity, use either `anisotropy.blangy()` or `anisotropy.ruger()`. Args: vp (ndarray): The P-wave velocity; float or 1D array length m. vs (ndarray): The S-wave velocity; float or 1D array length m. rho (ndarray): The density; float or 1D array length m. theta (ndarray): The incidence angle; float or 1D array length n. axis (int): The dimension along which to compute the reflectivity. mode (str): 'same' means the output will be the same shape. 'valid' means that only strictly valid reflectivities are computed so the result will be one sample shorted in the `axis` dimension. method (str): The reflectivity equation to use; one of: - 'scattering_matrix': scattering_matrix - 'zoeppritz_element': zoeppritz_element - 'zoeppritz': zoeppritz - 'zoeppritz_rpp': zoeppritz_rpp - 'akirichards': akirichards - 'akirichards_alt': akirichards_alt - 'fatti': fatti - 'shuey': shuey - 'bortfeld': bortfeld - 'hilterman': hilterman Notes: - scattering_matrix gives the full solution - zoeppritz_element gives a single element which you specify - zoeppritz returns RPP element only; use zoeppritz_rpp instead - zoeppritz_rpp is faster than zoeppritz or zoeppritz_element - You can also pass a function with the same API as these built-in functions. Returns: ndarray. The result of running the specified method on the inputs. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m-1 array (for float inputs and one angle), or an m-1 x n array (for array inputs and an array of angles). """ methods = { 'scattering_matrix': scattering_matrix, 'zoeppritz_element': zoeppritz_element, 'zoeppritz': zoeppritz, 'zoeppritz_rpp': zoeppritz_rpp, 'akirichards': akirichards, 'akirichards_alt': akirichards_alt, 'fatti': fatti, 'shuey': shuey, 'bortfeld': bortfeld, 'hilterman': hilterman, } func = methods.get(method.casefold(), method) if axis < 0: axis = vp.ndim + axis if axis > 0: vp_ = np.moveaxis(np.asanyarray(vp, dtype=float), axis, 0) vs_ = np.moveaxis(np.asanyarray(vs, dtype=float), axis, 0) rho_ = np.moveaxis(np.asanyarray(rho, dtype=float), axis, 0) else: vp_, vs_, rho_ = vp, vs, rho if mode == 'same': pad_width = [(0, 1)] + (vp.ndim - 1) * [(0, 0)] vp_ = np.pad(vp_, pad_width=pad_width, mode='edge') vs_ = np.pad(vs_, pad_width=pad_width, mode='edge') rho_ = np.pad(rho_, pad_width=pad_width, mode='edge') rc = func(vp_[:-1], vs_[:-1], rho_[:-1], vp_[1:], vs_[1:], rho_[1:], theta) if axis > 0: return np.moveaxis(rc, 0, axis) else: return np.squeeze(rc)
[docs]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): vp1 = np.asanyarray(vp1, dtype=float) vs1 = np.asanyarray(vs1, dtype=float) + 1e-12 # Prevent singular matrix. rho1 = np.asanyarray(rho1, dtype=float) vp2 = np.asanyarray(vp2, dtype=float) vs2 = np.asanyarray(vs2, dtype=float) + 1e-12 # Prevent singular matrix. rho2 = np.asanyarray(rho2, dtype=float) new_shape = [-1] + vp1.ndim * [1] theta1 = theta1.reshape(*new_shape) if (np.nan_to_num(theta1) > np.pi/2.).any(): raise ValueError("Incidence angle theta1 must be less than 90 deg.") return func(vp1, vs1, rho1, vp2, vs2, rho2, theta1, **kwargs) return wrapper
[docs]def preprocess(func): """ Decorator to preprocess arguments for the reflectivity equations. 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 some arguments transformed. """ @wraps(func) def wrapper(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, **kwargs): # Interpret tuple for theta1 as a linspace. if isinstance(theta1, tuple): if len(theta1) == 2: start, stop = theta1 theta1 = np.linspace(start, stop, num=stop+1) elif len(theta1) == 3: start, stop, step = theta1 steps = (stop / step) + 1 theta1 = np.linspace(start, stop, num=steps) else: raise TypeError("Expected 2 or 3 parameters for theta1 expressed as range.") # Convert theta1 to radians and complex numbers. theta1 = np.radians(theta1).astype(complex) return func(vp1, vs1, rho1, vp2, vs2, rho2, theta1, **kwargs) return wrapper
[docs]@preprocess @vectorize def scattering_matrix(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0): """ Full Zoeppritz solution, considered the definitive solution. Calculates the angle dependent p-wave reflectivity of an interface between two mediums. Originally written by: Wes Hamlyn, vectorized by Agile. Returns the complex reflectivity. Args: vp1 (float): The upper P-wave velocity. vs1 (float): The upper S-wave velocity. rho1 (float): The upper layer's density. vp2 (float): The lower P-wave velocity. vs2 (float): The lower S-wave velocity. rho2 (float): The lower layer's density. theta1 (ndarray): The incidence angle; float or 1D array length n. Returns: ndarray. The exact Zoeppritz solution for all modes at the interface. A 4x4 array representing the scattering matrix at the incident angle theta1. """ theta1 *= np.ones_like(vp1) p = np.sin(theta1) / vp1 # Ray parameter. theta2 = np.arcsin(p * vp2) # Trans. angle of P-wave. phi1 = np.arcsin(p * vs1) # Refl. angle of converted S-wave. phi2 = np.arcsin(p * vs2) # Trans. angle of converted S-wave. # Matrix form of Zoeppritz equations... M & N are matrices. M = np.array([[-np.sin(theta1), -np.cos(phi1), np.sin(theta2), np.cos(phi2)], [np.cos(theta1), -np.sin(phi1), np.cos(theta2), -np.sin(phi2)], [2 * rho1 * vs1 * np.sin(phi1) * np.cos(theta1), rho1 * vs1 * (1 - 2 * np.sin(phi1) ** 2), 2 * rho2 * vs2 * np.sin(phi2) * np.cos(theta2), rho2 * vs2 * (1 - 2 * np.sin(phi2) ** 2)], [-rho1 * vp1 * (1 - 2 * np.sin(phi1) ** 2), rho1 * vs1 * np.sin(2 * phi1), rho2 * vp2 * (1 - 2 * np.sin(phi2) ** 2), -rho2 * vs2 * np.sin(2 * phi2)]]) N = np.array([[np.sin(theta1), np.cos(phi1), -np.sin(theta2), -np.cos(phi2)], [np.cos(theta1), -np.sin(phi1), np.cos(theta2), -np.sin(phi2)], [2 * rho1 * vs1 * np.sin(phi1) * np.cos(theta1), rho1 * vs1 * (1 - 2 * np.sin(phi1) ** 2), 2 * rho2 * vs2 * np.sin(phi2) * np.cos(theta2), rho2 * vs2 * (1 - 2 * np.sin(phi2) ** 2)], [rho1 * vp1 * (1 - 2 * np.sin(phi1) ** 2), -rho1 * vs1 * np.sin(2 * phi1), - rho2 * vp2 * (1 - 2 * np.sin(phi2) ** 2), rho2 * vs2 * np.sin(2 * phi2)]]) M_ = np.moveaxis(np.squeeze(M), [0, 1], [-2, -1]) A = np.linalg.inv(M_) N_ = np.moveaxis(np.squeeze(N), [0, 1], [-2, -1]) Z_ = np.matmul(A, N_) return np.transpose(Z_, axes=list(range(Z_.ndim - 2)) + [-1, -2])
[docs]def zoeppritz_element(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, element='PdPu'): """ Returns any mode reflection coefficients from the Zoeppritz scattering matrix. Pass in the mode as element, e.g. 'PdSu' for PS. Wraps scattering_matrix(). Returns the complex reflectivity. Args: vp1 (float): The upper P-wave velocity. vs1 (float): The upper S-wave velocity. rho1 (float): The upper layer's density. vp2 (float): The lower P-wave velocity. vs2 (float): The lower S-wave velocity. rho2 (float): The lower layer's density. theta1 (ndarray): The incidence angle; float or 1D array length n. element (str): The name of the element to return, must be one of: 'PdPu', 'SdPu', 'PuPu', 'SuPu', 'PdSu', 'SdSu', 'PuSu', 'SuSu', 'PdPd', 'SdPd', 'PuPd', 'SuPd', 'PdSd', 'SdSd', 'PuSd', 'SuSd'. Returns: ndarray. Array length n of the exact Zoeppritz solution for the specified modes at the interface, at the incident angle theta1. """ elements = np.array([['PdPu', 'SdPu', 'PuPu', 'SuPu'], ['PdSu', 'SdSu', 'PuSu', 'SuSu'], ['PdPd', 'SdPd', 'PuPd', 'SuPd'], ['PdSd', 'SdSd', 'PuSd', 'SuSd']]) Z = scattering_matrix(vp1, vs1, rho1, vp2, vs2, rho2, theta1).T return np.squeeze(Z[np.where(elements == element)].T)
[docs]def zoeppritz(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0): """ Returns the PP reflection coefficients from the Zoeppritz scattering matrix. Wraps zoeppritz_element(). Returns the complex reflectivity. Args: vp1 (float): The upper P-wave velocity. vs1 (float): The upper S-wave velocity. rho1 (float): The upper layer's density. vp2 (float): The lower P-wave velocity. vs2 (float): The lower S-wave velocity. rho2 (float): The lower layer's density. theta1 (ndarray): The incidence angle; float or 1D array length n. Returns: ndarray. Array length n of the exact Zoeppritz solution for the specified modes at the interface, at the incident angle theta1. """ return zoeppritz_element(vp1, vs1, rho1, vp2, vs2, rho2, theta1, 'PdPu')
[docs]@preprocess @vectorize def zoeppritz_rpp(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0): """ Exact Zoeppritz from expression. This is useful because we can pass arrays to it, which we can't do to scattering_matrix(). Dvorkin et al. (2014). Seismic Reflections of Rock Properties. Cambridge. Returns the complex reflectivity. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. Returns: ndarray. The exact Zoeppritz solution for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ p = np.sin(theta1) / vp1 # Ray parameter theta2 = np.arcsin(p * vp2) phi1 = np.arcsin(p * vs1) # Reflected S phi2 = np.arcsin(p * vs2) # Transmitted S a = rho2 * (1 - 2 * np.sin(phi2)**2.) - rho1 * (1 - 2 * np.sin(phi1)**2.) b = rho2 * (1 - 2 * np.sin(phi2)**2.) + 2 * rho1 * np.sin(phi1)**2. c = rho1 * (1 - 2 * np.sin(phi1)**2.) + 2 * rho2 * np.sin(phi2)**2. d = 2 * (rho2 * vs2**2 - rho1 * vs1**2) E = (b * np.cos(theta1) / vp1) + (c * np.cos(theta2) / vp2) F = (b * np.cos(phi1) / vs1) + (c * np.cos(phi2) / vs2) G = a - d * np.cos(theta1)/vp1 * np.cos(phi2)/vs2 H = a - d * np.cos(theta2)/vp2 * np.cos(phi1)/vs1 D = E*F + G*H*p**2 rpp = (1/D) * (F*(b*(np.cos(theta1)/vp1) - c*(np.cos(theta2)/vp2)) \ - H*p**2 * (a + d*(np.cos(theta1)/vp1)*(np.cos(phi2)/vs2))) return np.squeeze(rpp)
[docs]@preprocess @vectorize def akirichards(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ The Aki-Richards approximation to the reflectivity. This is the formulation from Avseth et al., Quantitative seismic interpretation, Cambridge University Press, 2006. Adapted for a 4-term formula. See http://subsurfwiki.org/wiki/Aki-Richards_equation. Returns the complex reflectivity. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The Aki-Richards approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta2 = np.arcsin(vp2/vp1*np.sin(theta1)) drho = rho2-rho1 dvp = vp2-vp1 dvs = vs2-vs1 meantheta = (theta1+theta2) / 2.0 rho = (rho1+rho2) / 2.0 vp = (vp1+vp2) / 2.0 vs = (vs1+vs2) / 2.0 # Compute the coefficients w = 0.5 * drho/rho x = 2 * (vs/vp1)**2 * drho/rho y = 0.5 * (dvp/vp) z = 4 * (vs/vp1)**2 * (dvs/vs) # Compute the terms term1 = w term2 = -1 * x * np.sin(theta1)**2 term3 = y / np.cos(meantheta)**2 term4 = -1 * z * np.sin(theta1)**2 if terms: fields = ['term1', 'term2', 'term3', 'term4'] AkiRichards = namedtuple('AkiRichards', fields) return AkiRichards(np.squeeze([term1 for _ in theta1]), np.squeeze(term2), np.squeeze(term3), np.squeeze(term4) ) else: return np.squeeze(term1 + term2 + term3 + term4)
[docs]@preprocess @vectorize def akirichards_alt(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ This is another formulation of the Aki-Richards solution. See http://subsurfwiki.org/wiki/Aki-Richards_equation Returns the complex reflectivity. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The Aki-Richards approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta2 = np.arcsin(vp2/vp1*np.sin(theta1)) drho = rho2-rho1 dvp = vp2-vp1 dvs = vs2-vs1 theta = (theta1+theta2)/2.0 rho = (rho1+rho2)/2.0 vp = (vp1+vp2)/2.0 vs = (vs1+vs2)/2.0 # Compute the three terms term1 = 0.5 * (dvp/vp + drho/rho) term2 = (0.5*dvp/vp-2*(vs/vp)**2*(drho/rho+2*dvs/vs)) * np.sin(theta)**2 term3 = 0.5 * dvp/vp * (np.tan(theta)**2 - np.sin(theta)**2) if terms: fields = ['term1', 'term2', 'term3'] AkiRichards = namedtuple('AkiRichards', fields) return AkiRichards(np.squeeze([term1 for _ in theta1]), np.squeeze(term2), np.squeeze(term3) ) else: return np.squeeze(term1 + term2 + term3)
[docs]@preprocess @vectorize def fatti(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ Compute reflectivities with Fatti's formulation of the Aki-Richards equation, which does not account for the critical angle. See Fatti et al. (1994), Geophysics 59 (9). Real numbers only. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The Fatti approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta1 = np.real(theta1) drho = rho2-rho1 rho = (rho1+rho2) / 2.0 vp = (vp1+vp2) / 2.0 vs = (vs1+vs2) / 2.0 dip = (vp2*rho2 - vp1*rho1)/(vp2*rho2 + vp1*rho1) dis = (vs2*rho2 - vs1*rho1)/(vs2*rho2 + vs1*rho1) d = drho/rho # Compute the three terms term1 = (1 + np.tan(theta1)**2) * dip term2 = -8 * (vs/vp)**2 * dis * np.sin(theta1)**2 term3 = -1 * (0.5 * np.tan(theta1)**2 - 2 * (vs/vp)**2 * np.sin(theta1)**2) * d if terms: fields = ['term1', 'term2', 'term3'] Fatti = namedtuple('Fatti', fields) return Fatti(np.squeeze(term1), np.squeeze(term2), np.squeeze(term3) ) else: return np.squeeze(term1 + term2 + term3)
[docs]@preprocess @vectorize def shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False, return_gradient=False): """ Compute Shuey approximation with 3 terms. http://subsurfwiki.org/wiki/Shuey_equation Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. return_gradient (bool): Whether to return a tuple of the intercept and gradient (i.e. the second term divided by sin^2(theta)). Returns: ndarray. The Aki-Richards approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta1 = np.real(theta1) drho = rho2-rho1 dvp = vp2-vp1 dvs = vs2-vs1 rho = (rho1+rho2)/2.0 vp = (vp1+vp2)/2.0 vs = (vs1+vs2)/2.0 # Compute three-term reflectivity r0 = 0.5 * (dvp/vp + drho/rho) g = 0.5 * dvp/vp - 2 * (vs**2/vp**2) * (drho/rho + 2 * dvs/vs) f = 0.5 * dvp/vp term1 = r0 term2 = g * np.sin(theta1)**2 term3 = f * (np.tan(theta1)**2 - np.sin(theta1)**2) if return_gradient: fields = ['intercept', 'gradient'] Shuey = namedtuple('Shuey', fields) return Shuey(np.squeeze(r0), np.squeeze(g)) elif terms: fields = ['R0', 'Rg', 'Rf'] Shuey = namedtuple('Shuey', fields) return Shuey(np.squeeze([term1 for _ in theta1]), np.squeeze(term2), np.squeeze(term3) ) else: return np.squeeze(term1 + term2 + term3)
[docs]@deprecated('Please use shuey() instead.') def shuey2(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0): """ Compute Shuey approximation with 2 terms. Wraps `shuey()`. Deprecated, use `shuey()` instead. """ r, g, _ = shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1, terms=True) return r + g
[docs]@deprecated('Please use shuey() instead.') def shuey3(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ Compute Shuey approximation with 3 terms. Wraps `shuey()`. Deprecated, use `shuey()` instead. """ return shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1)
[docs]@preprocess @vectorize def bortfeld(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ Compute Bortfeld approximation with three terms. http://sepwww.stanford.edu/public/docs/sep111/marie2/paper_html/node2.html Real numbers only. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The 3-term Bortfeld approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta1 = np.real(theta1) drho = rho2-rho1 dvp = vp2-vp1 dvs = vs2-vs1 rho = (rho1+rho2)/2.0 vp = (vp1+vp2)/2.0 vs = (vs1+vs2)/2.0 k = (2 * vs/vp)**2 rsh = 0.5 * (dvp/vp - k*drho/rho - 2*k*dvs/vs) # Compute three-term reflectivity term1 = 0.5 * (dvp/vp + drho/rho) term2 = rsh * np.sin(theta1)**2 term3 = 0.5 * dvp/vp * np.tan(theta1)**2 * np.sin(theta1)**2 if terms: fields = ['term1', 'term2', 'term3'] Bortfeld = namedtuple('Bortfeld', fields) return Bortfeld(np.squeeze([term1 for _ in theta1]), np.squeeze(term2), np.squeeze(term3) ) else: return np.squeeze(term1 + term2 + term3)
[docs]@deprecated('Please use bortfeld() instead.') def bortfeld2(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ The 2-term Bortfeld approximation for ava analysis. Wraps `shuey()`. Deprecated, use `bortfeld()` instead. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The 2-term Bortfeld approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta1 = np.radians(theta1) theta2 = np.arcsin(vp2/vp1*np.sin(theta1)) term1 = 0.5 * np.log((vp2*rho2*np.cos(theta1)) / (vp1*rho1*np.cos(theta2))) svp2 = (np.sin(theta1)/vp1)**2 dvs2 = (vs1**2-vs2**2) term2 = svp2 * dvs2 * (2+np.log(rho2/rho1)/np.log(vs2/vs1)) if terms: return term1, term2 else: return (term1 + term2)
[docs]@deprecated('Please use bortfeld() instead.') def bortfeld3(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): return bortfeld(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1)
[docs]@preprocess @vectorize def hilterman(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False): """ Not recommended, only seems to match Zoeppritz to about 10 deg. Hilterman (1989) approximation from Mavko et al. Rock Physics Handbook. According to Dvorkin: "arguably the simplest and a very convenient [approximation]." At least for small angles and small contrasts. Real numbers only. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. terms (bool): Whether or not to return a tuple of the terms of the equation. The first term is the acoustic impedance. Returns: ndarray. The Hilterman approximation for P-P reflectivity at the interface. Will be a float (for float inputs and one angle), a 1 x n array (for float inputs and an array of angles), a 1 x m array (for float inputs and one angle), or an n x m array (for array inputs and an array of angles). """ theta1 = np.real(theta1) ip1 = vp1 * rho1 ip2 = vp2 * rho2 rp0 = (ip2 - ip1) / (ip2 + ip1) pr2, pr1 = moduli.pr(vp2, vs2), moduli.pr(vp1, vs1) pravg = (pr2 + pr1) / 2. pr = (pr2 - pr1) / (1 - pravg)**2. term1 = rp0 * np.cos(theta1)**2. term2 = pr * np.sin(theta1)**2. if terms: fields = ['term1', 'term2'] Hilterman = namedtuple('Hilterman', fields) return Hilterman(np.squeeze(term1), np.squeeze(term2)) else: return np.squeeze(term1 + term2)
[docs]def blangy(vp1, vs1, rho1, vp2, vs2, rho2, d1=0, e1=0, d2=0, e2=0, theta1=0): """Implements the Blangy equation with the same interface as the other reflectivity equations. Wraps bruges.anisotropy.blangy(), which you may prefer to use directly. Args: vp1 (ndarray): The upper P-wave velocity; float or 1D array length m. vs1 (ndarray): The upper S-wave velocity; float or 1D array length m. rho1 (ndarray): The upper layer's density; float or 1D array length m. vp2 (ndarray): The lower P-wave velocity; float or 1D array length m. vs2 (ndarray): The lower S-wave velocity; float or 1D array length m. rho2 (ndarray): The lower layer's density; float or 1D array length m. d1 (ndarray): The upper delta; float or 1D array length m. e1 (ndarray): The upper epsilon; float or 1D array length m. d2 (ndarray): The lower delta; float or 1D array length m. e2 (ndarray): The lower epsilon; float or 1D array length m. theta1 (ndarray): The incidence angle; float or 1D array length n. Returns: ndarray. The Blangy approximation for P-P reflectivity at the interface. Wraps `anisotropy.blangy()`. """ _, anisotropic = anisotropy.blangy(vp1, vs1, rho1, d1, e1, # UPPER vp2, vs2, rho2, d2, e2, # LOWER theta1) return anisotropic