Métodos potenciales de prospección, FCAG, 2024.
Generamos el vector de posiciones de estaciones $\mathbf{x}$, alturas de estaciones $\mathbf{z}$, y de mediciones de anomalía de aire libre (AAL) observadas $\mathbf{dg}$.
# Constantes:
G = 6.674e-11 # unidades: m3 / ( kg s2 )
CB0 = 2. * np.pi * G
CB0 = CB0 * ... # unidades: [mGal] para z en [m] y densidad en [kg/m3].
# Datos:
x = np.array([0.0,106.6,180.0,260.0,593.3]) # posiciones de las estaciones en el perfil: [m]
z = np.array([421.0,431.0,434.5,426.5,409.0]) # alturas de las estaciones: [m]
dg = np.array([51.35,52.40,52.72,51.39,50.13]) # AAL: [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("$A_{AL}$")
plt.scatter(x,dg)
plt.xlabel("$x$ [m]")
plt.ylabel("$A_{AL}$ [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.
def correlation_coefficient(x,y):
r = ...
return r
den = [1000,...,3600] # kg/m3
r = np.zeros(len(den)) # Un coeficiente por cada valor de densidad, para graficar.
cb = np.zeros((len(den),len(dg))) # Se guardan las curvas para cada densidad para luego graficar.
for j in range(len(den)):
# Corrección de Bouguer para den[j]
cb[j,:] = dg - 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 = [dg - 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] = 1000.0 => r = 0.93 σ [kg/m3] = 1500.0 => r = 0.84 σ [kg/m3] = 2000.0 => r = 0.55 σ [kg/m3] = 2100.0 => r = 0.43 σ [kg/m3] = 2200.0 => r = 0.29 σ [kg/m3] = 2300.0 => r = 0.13 σ [kg/m3] = 2350.0 => r = 0.04 σ [kg/m3] = 2400.0 => r = -0.05 σ [kg/m3] = 2500.0 => r = -0.22 σ [kg/m3] = 2670.0 => r = -0.47 σ [kg/m3] = 2800.0 => r = -0.61 σ [kg/m3] = 3000.0 => r = -0.75 σ [kg/m3] = 3500.0 => r = -0.9
Mínimo de la correlación: 0.04 Densidad [km/m3]: 2350.0
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.0 => r* = 0.04 ------------------------------------------
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] = 2371.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] = 2370.88
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 = 0.947
Interpretamos este resultado así: suponer un modelo lineal representa una reducción del 95% 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 libería NumPy:
from numpy.linalg import lstsq
a, b = lstsq(...)[0]