Métodos potenciales de prospección, FCAG, 2024.
Copiamos los vectores con las posiciones de estaciones $\mathbf{x}$, las alturas de estaciones $\mathbf{z}$, y las mediciones de la gravedad $\mathbf{g}$.
# Datos:
x = np.array([...]) # posiciones de las estaciones en el perfil: [m]
z = np.array([...]) # alturas de las estaciones: [m]
g = np.array([...]) # mediciones de gravedad: [mGal]
Graficamos,
plt.figure()
plt.subplot(1,2,1)
plt.title("Estaciones")
plt.plot(x,z)
plt.scatter(x,z)
plt.xlabel("$x$ [m]")
plt.ylabel("$z$ [m]")
plt.subplot(1,2,2)
plt.title("$g_{z}$")
plt.scatter(x,g)
plt.xlabel("$x$ [m]")
plt.ylabel("$g_{z}$ [mGal]")
plt.show()
Dada la anomalía de aire libre ($A_{AL}) $ en una estación $j$ a una altura $z_j$, se calcula la anomalía de Bouguer ($A_{B}$) con distintas densidades $\color{blue} \sigma$ en el término de la corrección de Bouguer,
donde $G=6.674\times 10^{-11}$ m$^3$/ kg s$^2$ es la constante de gravitación universal de Newton.
El método de Nettleton considera que la densidad que visualmente genera el perfil menos correlacionado con la topografía es una buena aproximación a la densidad del rasgo topográfico.
En este ejemplo calculamos, para cada $\sigma$, el coeficiente de correlación entre la topografía $\mathbf{z}$ y la anomalía $A_{B}$. De esta forma podemos guiarnos al momento de elegir la densidad que representa el rasgo topográfico de interés. El valor de densidad elegido debería generar un correlación próxima a cero.
Recordamos que el coeficiente de correlación entre dos series $\mathbf{x},\mathbf{y}$ está dado por
$$r = \frac{\sum_k \left(x_k-\overline{x}\right)\left(y_k-\overline{y}\right)}{\sqrt{\sum_k \left(x_k-\overline{x}\right)^2\sum_k\left(y_k-\overline{y}\right)^2 }} .$$
Podemos verificar el valor del coeficiente de correlación usando apropiadamente la subrutina numpy.corrcoefdel paquete NumPy, o también correlation de la librería scipy.spatial.correlation.
# Constantes:
G = 6.674e-11 # [m^3 / ( kg s^2 )]
CB0 = (2. * np.pi * G)* 1e5 # [mGal] para z en [m] y densidad en [kg/m3].
CAL = 0.3086 # [mGal / m]
def correlation_coefficient(x,y):
# subrutina para el cálculo del coeficiente de correlación
r = ...
return r
gamma = 981000. # [mGal]
aal = ... - gamma # Anomalía de aire libre para el método de Nettleton
# Método de Nettleton
den = [2000,...,3600] # valores de densidades plausibles del rasgo topográfico [kg/m3]
r = np.zeros(len(den)) # Un coeficiente por cada valor de densidad, para graficar.
cb = np.zeros((len(den),len(aal))) # Se guardan las curvas para cada densidad para luego graficar.
# Corrección de Bouguer para den[j]:
for j in range(len(den)):
cb[j,:] = aal - CB0 * den[j] * z # [mGal] con CB0 = 2*np.pi*G
# Coeficiente de correlación para den[j]:
r[j] = correlation_coefficient(z,cb[j,:])
Otra forma de realizar lo mismo (por compresión):
cb = [aal - CB0 * den_k * z for den_k in den]
r = [correlation_coefficient( z , cb_k ) for cb_k in cb]
El método visual se reduce entonces a:
# Método analítico de la correlación
index_min = np.argmin(np.abs(r)) # índice donde la correlación es mínima
print("Mínimo de la correlación: ", r[index_min])
print("Densidad [km/m3]: ", den[index_min])
=========================================== Método de Nettleton =========================================== σ [kg/m3] = 2300 => r = 0.99 σ [kg/m3] = 2305 => r = 0.99 σ [kg/m3] = 2310 => r = 0.99 σ [kg/m3] = 2315 => r = 0.99 σ [kg/m3] = 2320 => r = 0.98 σ [kg/m3] = 2325 => r = 0.97 σ [kg/m3] = 2330 => r = 0.96 σ [kg/m3] = 2335 => r = 0.92 σ [kg/m3] = 2340 => r = 0.82 σ [kg/m3] = 2345 => r = 0.44 σ [kg/m3] = 2350 => r = -0.39 σ [kg/m3] = 2355 => r = -0.8 σ [kg/m3] = 2360 => r = -0.92 σ [kg/m3] = 2365 => r = -0.95 σ [kg/m3] = 2370 => r = -0.97 σ [kg/m3] = 2375 => r = -0.98 σ [kg/m3] = 2380 => r = -0.99 σ [kg/m3] = 2385 => r = -0.99 σ [kg/m3] = 2390 => r = -0.99 σ [kg/m3] = 2395 => r = -0.99
Mínimo de la correlación: -0.39 Densidad [km/m3]: 2350
Nota. Para simplificar el análisis de los gráficos, tanto la topografía como la anomalía han sido escalados. En este ejemplo se tomó la siguiente normalización: $ \hat{x} = \left(x - \overline{x}\right)/\: \text{std}(x)$.
Del análisis visual y los resultados numéricos es posible resolver la densidad buscada:
index_min = np.argmin(np.abs(r)) # índice donde la correlación es mínima
print("Mínimo de la correlación: ", r[index_min])
print("Densidad [km/m3]: ", den[index_min])
------------------------------------------ Densidad seleccionada: ------------------------------------------ σ* [kg/m3] = 2350 => r* = -0.39 ------------------------------------------
Notemos que el valor de la densidad seleccionada pertenece al juego de valores que utilizamos para evaluar la corrección de Bouguer. Es decir, la respuesta depende de los valores de densidad que hemos considerado para la búsqueda exhaustiva del mínimo.
Dada $\sigma_0$, resolvemos analíticamente la corrección ${\color{blue} \Delta \color{blue} \sigma}$ que cancela al coeficiente de correlación entre la anomalía de Bouguer y la topografía $\mathbf{z}$:
$${\color{blue} \Delta \color{blue} \sigma} \quad\text{tal que}\quad r\left(\:\mathbf{z},A_B\left(\mathbf{z},\sigma_0+{\color{blue} \Delta \color{blue} \sigma}\right)\:\right) = 0. $$Como se detalla en la teoria de la cátedra, este problema tiene una solución analítica. De la solución analítica obtenemos la corrección para determinar la densidad apropiada del rasgo topográfico:
Notemos que ${\color{blue} \Delta \color{blue} \sigma}$ puede ser tanto negativa (corrección por exceso), como positiva (corrección por defecto). La solución viene luego dada por
$$\sigma^* = \sigma_0+{\color{blue} \Delta \color{blue} \sigma},$$con $\sigma_0$ un valor inicial de densidad dado por el usuario.
El método se programa como subrutina y luego obtenemos la correción drho para un valor rho_0 inicial.
def minima_corr(AB,z):
...
return drho
rho_0 = ... # valor inicial de la densidad [km/m3]
cb = dg - CB0 * rho_0 * z # anomalía de Bouguer [mGal]
drho = minima_corr(cb,z) # corrección de la densidad [km/m3]
rho = rho_0 + drho # densidad estimada [km/m3]
print(rho)
=========================================== Método de mínima correlación =========================================== σ*[kg/m3] = 2348.0 => r* = -0.0
El resultado analítico es coherente con el obtenido por el método de Nettleton. En este caso llegamos a la solución en un sólo paso y el valor obtenido es independiente del valor inicial de la densidad $\sigma_0$.
Este método involucra una regresión lineal entre los valores de la anomalía de aire libre y la altura topográfica: $\mathbf{dg} = a\mathbf{z}+\mathbf{b}$. De la regresión podemos obtener la pendiente del ajuste, $a^*$, y con ella determinar el valor buscado de la densidad $\color{blue}\sigma$. En la teoría de la Cátedra, se describe como arribar a la relación empleada para resolver el problema:
=========================================== Método de la anomalía de aire libre =========================================== Densidad estimada [kg/m3] = 2347.67
Por último, aplicamos el método de Siegert, en su versión modificada analizada en la Teoría de la cátedra. Dados los valores de gravedad $\mathbf{g}$ y las alturas de las estaciones $\mathbf{h}$, se obtienen las diferencias $\Delta \mathbf{g}, \Delta \mathbf{h}$. Con ellas se calcula la expresión de la constante de la corrección combinada de aire libre y Bouguer, $\mathcal{K}$. Esto resulta en la expresión analítica:
Esta expresión para $\mathcal{K}$ resulta de obtener la pendiente de la recta de ajuste del gráfico $\Delta \mathbf{g}$ vs. $\Delta \mathbf{h}$ por mínimos cuadrados.
Calculada la constante de la corrección combinada $\mathcal{K}$, se obtiene ${\color{blue}{\sigma_b}}$ de la relación:
# Método de Siegert
dg = ... # diferencias dg = [ g[1] - g[0], ... ]
dh = ... # diferencias dh = [ h[1] - h[0], ... ]
K = ... # constante de la corrección combinada
sigma_b = ... # valor estimado del rasgo topográfico
=========================================== Método de Siegert =========================================== Densidad estimada [kg/m3] = 2346.13
Eso es todo por hoy.
$\newcommand{\mean}[1]{<\!{#1}\!>}$$\newcommand{\B}[1]{\mathbf{#1}}$
Sólo por curiosidad podemos calcular el estadístico $R^2$ entre el modelo lineal obtenido $a^{*}\mathbf{z} + \mathbf{b}^{*}$, $\mathbf{b}^{*}_j = b^{*}$, y el modelo que consiste en considerar el valor medio de las mediciones $\mean{\B{g}}$:
$$R^2 = \dfrac{\textrm{SSR}(\mathbf{g},\mean{\B{g}})-\textrm{SSR}(\B{g},a^*\B{z}+\B{b}^*)}{\textrm{SSR}(\B{g},\mean{\B{g}})},$$donde $\textrm{SSR}(\B{x},\B{y})$ es la suma de los cuadrados de la diferencias entre los elementos de $\B{x}$ e $\B{y}$.
def SSR(x,y):
return np.sum(np.square(x-y))
El ajuste resulta en un R2 = 1.0
Interpretamos este resultado así: suponer un modelo lineal representa una reducción del 100% del error de ajuste de suponer el valor medio como modelo de las observaciones.
En esta notebook empleamos el resultado analítico para obtener la pendiente del ajuste lineal. Puede verificarse el resultado mediante la subrutina scipy.stats.linregress() de la libreria SciPy:
from scipy.stats import linregress
a, b, _, _, _ = linregress(...)
También es posible utilizar numpy.linalg.lstsq() de la librería NumPy:
from numpy.linalg import lstsq
a, b = lstsq(...)[0]