Métodos potenciales de prospección, FCAG, 2023.
En primer lugar, vamos a generar una anomalía gravimétrica sintética que utilizaremos para simular mediciones realizadas en el campo. Para este ejemplo, el modelo directo es la anomalía gravimétrica vertical producida en superficie por un cilindro 2D horizontal y homogéneo, de radio $R$, contraste de densidad $\Delta\sigma$ y con eje centrado en $(x_0,z_0)$:
La expresión es lineal en $\Delta\sigma$ y no lineal en $R$, $x_0$ y $z_0$. La coordenada de cada estación de observación viene dada por $x$.
def model(x,m):
...
density_constrast, radius, x0, z0 = m
...
return gz # [mGal]
El vector columna ($M\times 1$) con $M=4$, $\mathbf{m}=(m_1,m_2,m_3,m_4)^T = (\Delta\sigma,R,x_0,z_0)^T$, contiene los parámetros físicos y geométricos del modelo y lo llamaremos vector de parámetros del modelo.
Los parámetros que definen el dato sintético los agrupamos como elementos del vector $\mathbf{m}^*=(\Delta\sigma^*,R^*,x_0^*,z_0^*)^T$:
m_true = ...
Ahora tomamos una serie de $N$ estaciones, cada una ubicada en la posición $x_i$, con $i=0,\cdots,N-1$. Los valores de la anomalía en estas $N$ estaciones las agrupamos en un vector columna $Ṇ\times 1$, $\mathbf{d}^\text{o}$, que llamaremos vector de datos observados.
x = ...
d0 = model(x, m = m_true)
A continuación graficamos para visualizar el dato observado.
plt.plot(x, d0,...)
Una vez obtenidas las observaciones, buscaremos estimar los parámetros del modelo propuesto por medio de la minimización de una función objetivo. La función objetivo es una medida de la norma euclídea del vector de las diferencias entre el vector de los datos observados $\mathbf{d}^\text{o}$ y el vector con los datos calculados $\mathbf{d}^\text{c}$ por medio de un modelo $\mathbf{d}^c=\mathbf{f}(\mathbf{m})$.
es la función del problema directo, en este caso el modelo de la anomalía gravimetrica dada por un cilindro 2D.
En este ejemplo sencillo, utilizaremos como modelo a la misma expresión que utilizamos para generar el dato observado sintético.
Construimos la función objetivo como
\begin{equation} q(\Delta\sigma,R,x_0,z_0) = q(\mathbf{m}) := \frac{1}{2}||\,\mathbf{d}^\text{o}-\mathbf{d}^\text{c}\,||^2= \frac{1}{2}||\,\mathbf{d}^\text{o}-\mathbf{f}(\mathbf{m})\,||^2= \frac{1}{2}||\,\mathbf{e}\,||^2. \end{equation}El vector $\mathbf{e}(\mathbf{m})=\mathbf{d}^\text{o} - \mathbf{d}^\text{c}$ se conoce como vector de error o desajuste.
Calculamos la función objetivo y graficamos para distintos valores de $\Delta\sigma$, dejando el resto de los parámetros fijos (en los valores dados en $\mathbf{m}^*$).
def funcion_de_costo(error)
q = error.T @ error
return q / 2
Procedemos a buscar los parámetros que minimizan la función objetivo por medio de un algoritmo de descenso de gradiente.
Nos decidimos a probar el siguiente método:
El gradiente de la función de costo respecto de los parámetros del modelo viene dado por
La expresión $A=\frac{\partial\mathbf{f}(\mathbf{m})}{\partial\mathbf{m}}$ es la matriz jacobiana ($N\times M$) del modelo respecto de los parámetros evaluada en $\mathbf{m}$. Elemento por elemento es $A_{ij}=\frac{\partial}{\partial m_j}f(\mathbf{m})_i$.
A[i,:] está asociada a la estación $x_i$,A[:,j] contiene la derivada de la función del problema directo respecto del parámetro $m_j$.En el problema considerado la matriz jacobiana tiene dimensiones $N{\times}4$, siendo $N$ el número de estaciones de medición y $M=4$ el número de parámetros a invertir:
Dado que \begin{equation*} f(x;\mathbf{m} )\equiv g_z(x;\Delta\sigma,R,x_0,z_0) = 2\pi{G}\Delta\sigma\,R^2 \frac{z_0}{(x-x_0)^2+z_0^2}, \end{equation*}
las derivadas parciales del modelo respecto de los parámetros son:
\begin{align}{} \frac{\partial}{\partial \Delta\sigma}\,f(\mathbf{m})_{\color{blue}i} & = 2\pi{G}\,R^2 \frac{z_0}{({\color{blue}{x_i}}-x_0)^2+z_0^2}; \\ \frac{\partial}{\partial R}\,f(\mathbf{m})_{\color{blue}i} & = 4\pi{G}\Delta\sigma\,R\frac{z_0}{({\color{blue}{x_i}}-x_0)^2+z_0^2}; \\ \frac{\partial}{\partial x_0}\,f(\mathbf{m})_{\color{blue}i} & = 4\pi{G}\Delta\sigma\,R^2 \frac{({\color{blue}{x_i}}-x_0)\,z_0}{\left(({\color{blue}{x_i}}-x_0)^2+z_0^2\right)^2};\\ \frac{\partial}{\partial z_0}\,f(\mathbf{m})_{\color{blue}i} & = 2\pi{G}\Delta\sigma\,R^2 \frac{({\color{blue}{x_i}}-x_0)^2-z_0^2}{\left(({\color{blue}{x_i}}-x_0)^2+z_0^2\right)^2}. \end{align}Evaluando las derivadas parciales en $\mathbf{m}$ y para todos los $x=x_i$ construimos la matriz jacobiana. La matriz jacobiana nos permite luego evaluar el gradiente de la función de costo $q$.
def jacobiano(x,m):
...
dfdm1 = ...
dfdm2 = ...
dfdm3 = ...
dfdm4 = ...
...
A[:,0], A[:,1], A[:,2], A[:,3] = dfdm1,dfdm2,dfdm3,dfdm4
return A
def gradiente_de_la_funcion_de_costo(A, error):
dq = - A.T @ error
return dq
La expresión $\mathbf{e}\,\frac{\partial\mathbf{f}({\mathbf{m}})}{\partial\mathbf{m}}=A^T\mathbf{e}$, se conoce como producto vector por jacobiano (vector-jacobian product). Existen formas de calcular este producto sin tener que programar la matriz jacobiana explícitamente (es decir, sin programar las derivadas parciales del modelo). Esto puede realizarse, por ejemplo, con una aproximación por diferencias finitas, o por medio de la diferenciación automática.
Programamos a continuación el método de gradient descent.
Escribimos el método como una subrutina:
def GradientDescent(m0, epsilon, maxiter, ...):
...
m = m - epsilon * gradiente_de_la_funcion_de_costo(m,...)
return m
Elegimos parámetros iniciales del modelo directo e hiperparámetros ($\color{green}\epsilon$, maxiter) del algoritmo de inversión,
m0 = ...
epsilon = ...
maxiter = ...
y procedemos a realizar la inversión:
m = GradientDescent(m0, epsilon, maxiter,...)
vector inicial m0 : [ 100. 500. 28000. 500.] vector obtenido m : [ 614. 989. 30000. 1500.] vector de parámetros reales m*: [ 600. 1000. 30000. 1500.] Valor de q(m0) : 311.99666566304074 Valor de q(m*) : 1.1417897117250778e-08 Número de iteraciones : 1000
Analizamos los resultados por medio de una serie de gráficos.
Notemos que interpretando el dato observado, podríamos iniciar la inversión para $\Delta\sigma>0$ y con $x_0 \approx 30$ km. Un análisis de los cambios de curvatura del dato observado permitiría estimar valores iniciales para $R$ y $z_0$.
Por último, el algoritmo fallará en hallar el mínimo si, por ejemplo, tomamos $x_0=0$. El análisis de la topografía de la función de costo nos puede indicar la razón de ello.
Optamos por emplear un algoritmo de minimización de gradientes conjugados no lineal del paquete SciPy scipy.optimize.fmin_cg para resolver el problema planteado. Notamos que partimos de otra solución inicial $\mathbf{m}_0$ para obtener un resultado similar al método anterior desarrollado en esta práctica. De otra forma, el método arriba a otra solución admisible debido a la ambiguedad en invertir el producto $\Delta\sigma R^2$.
Podemos utilizar el algoritmo aproximando el gradiente de la función de costo por diferencias finitas mediante la opción fprime=None:
from scipy import optimize
def cost(m):
...
return q
m0 = ...
m = optimize.fmin_cg(f = cost, x0 = m0, fprime = None) # gradiente de la función de costo por diferencias finitas
...
Si proveemos la subrutina para el cálculo del gradiente de la función de costo, entonces hacemos lo siguiente:
def gradcost(m):
...
return dq
m = optimize.fmin_cg(f = cost, x0 = m0, fprime = gradcost) # gradiente dado por subrutina
...
Warning: Desired error not necessarily achieved due to precision loss.
Current function value: 0.000000
Iterations: 25
Function evaluations: 1006
Gradient evaluations: 199
SciPy
vector inicial m0 : [ 470. 500. 28000. 500.]
vector obtenido m : [ 668. 948. 30000. 1500.]
vector de parámetros reales m*: [ 600. 1000. 30000. 1500.]
Valor de q(m0) : 279.6118688934704
Valor de q(m*) : 1.781775439413038e-07
Número de iteraciones : 25
Realizamos un análisis similar utilizando el software de aprendizaje automático Tensorflow. En breve: definimos el modelo y minimizamos una función de costo por medio de un algoritmo de optimización basado en el descenso del gradiente. Para calcular el gradiente de la función de costo, Tensorflow utiliza diferenciación automática.
import tensorflow as tf
def modelGZ(tf.keras.Model):
def __init__(self, **kwargs):
...
def call(self, x):
...
return gz
loss = tf.losses.mean_squared_error # función de costo
optimizer = tf.keras.optimizers.SGD(learning_rate=1000.) # algoritmo de optimización
model = modelGZ() # generamos una instancia del modelo
model.compile(optimizer = optimizer, loss = loss)
model.fit(x,gz, epochs=150, batch_size=len(x),verbose=0) # inversión
...
Costo actual: 0.000000
Parámetros estimados: rho [kg/m3] = 831.3364 R [m] = 849.54663 x0 [m] = 30000.0 z0 [m] = 1499.9998
Como se describe en la teoria, sabemos que la inversión en los métodos potenciales es ambigua, en el sentido que una multitud de modelos pueden explicar el dato observado. Es más, una estimación de la ambiguedad presente en el resultado de la inversión es quizás tan importante como el modelo obtenido.
Como ejemplo de la ambiguedad en la inversión en los métodos potenciales, partimos de la misma solución inicial que la utilizada en el método de gradient descent. Una vez obtenida la solución $\mathbf{m}$, calculamos el producto $\Delta\sigma R^2$ para $\mathbf{m}$ y para los parámetros reales del modelo $\mathbf{m}^*$. Claramente, la solución sobreesestima el valor del constraste de densidad, y compensa desestimando el valor del radio del cilindro.
Warning: Desired error not necessarily achieved due to precision loss.
Current function value: 0.000000
Iterations: 24
Function evaluations: 782
Gradient evaluations: 154
SciPy
vector inicial m0 : [ 100. 500. 28000. 500.]
vector obtenido m : [ 2316. 509. 30000. 1500.]
vector de parámetros reales m*: [ 600. 1000. 30000. 1500.]
ΔσR^2 (parámetros reales en m*) : 600000000.0 [kg/m]
ΔσR^2 (parámetros obtenidos en m): 599990127.7215952 [kg/m]
Eso es todo por hoy.
Meju, M. A., Geophysical data analysis, 1994
Goodfellow et. al, Deep learning Book, 2016