Métodos potenciales de prospección, FCAG, 2024.
Muchos procesamientos de los métodos potenciales requieren conocer el campo de interés sobre una cuadrícula y a una misma altitud. El método de la capa equvialente fue propuesto por Dampney (1969) para interpolar y continuar datos gravimétricos. Esta técnica permite generar una cuadrícula de los datos magnéticos o gravitatiorios que han sido observados de forma irregular y a distintas altitudes, e.g., datos medidos en avión. $\renewcommand{\B}[1]{\mathbf{#1}}$
Como describe Blakely (1996), si puede encontrarse una distribución de fuentes que genera el campo observado, esa distribución de fuentes puede ser empleada para calcular el campo en cualquier punto por encima de las observaciones iniciales. El método de la capa equivalente es un proceso en dos pasos:
La distribución de las fuentes debe producir un campo que es armónico en el área de interés, que se desvanece hacia infinito y que reproduce el campo observado (Blakely, 1996).
Para implementar el método, la capa equivalente a la profundidad $h$ puede ser aproximada por un conjunto discreto de fuentes equivalentes, $m_j$, que en el caso magnético, consisten en dipolos magnéticos de volumen unitario.
Si ${\Delta T}_i$ es la anomalía magnética total observada en $(x_i,y_i,z_i)$, tenemos
\begin{equation} {\Delta T}_i =\sum_{j=1}^{N} m_j a_{ij}, \quad i=1,\ldots,N. \end{equation}Observamos que la anomalía modelada depende linealmente de las fuentes equivalentes. Podemos interpretar también que $m_j$ es la magnetización de la celda $j$ y $a_{ij}$ es el campo en el punto $i$ producido por la celda $j$ con magnetización unitaria (Blakely, 1996)
La estructura del 'kernel' $a_{ij}$ o función de Green en el caso considerado adopta la siguiente expresión:
\begin{equation} a_{ij} = a(x_i,y_i,z_i,x_j,y_j,z_j=h) = \frac{h-z_i}{\left[(x_i-x_j)^2+(y_i-y_j)^2+(z_i-h)^2\right]^{3/2}}. \end{equation}En la práctica, las coordenadas horizontales de la fuente $m_j$, $(x_j,y_j)$, se asumen iguales a las coordenadas de los puntos de observación $(x_i,y_i)$.
Si consideramos un sistema de coordenadas $x-y-z$ con $z>0$ hacia abajo (profundidad), $x>0$ hacia el Este e $y>0$ hacia el Norte, entonces a la hora de implementar este método, el signo de $h$ (profundidad de la capa equivalente) debe ser diferente al de las alturas de observación $z_i$.
Una vez obtenidas las fuentes equivalentes, es posible interpolar el campo en cualquier punto $(x,y,z)$ ya que:
\begin{equation} {\Delta T}(x,y,z) = \sum_{j=1}^{N} m_j\,a(x,y,z,x_j,y_j,h), \end{equation}lo que permite además efectuar continuación ascendente si $z>z_i \forall i$.
Matricialmente, la última ecuación puede escribirse para todos los puntos de interés como el sistema,
$$\mathbf{\Delta T} = \mathbf{A}\mathbf{m},$$siendo $\mathbf{A}$ la matriz de diseño, $\mathbf{m}$ el vector columna con las fuentes y $\mathbf{\Delta T}$ los valores de la anomalía modelada ${\Delta T}_i$ agrupados en un vector columna.
$\newcommand{\B}[1]{\mathbf{#1}}$ Las fuentes puntuales pueden obtenerse como la solución $\mathbf{m}^*$ dada por el método de mínimos cuadrados con regularización (a.k.a, Thikonov de orden cero):
\begin{equation} \mathbf{m}^*= \left( \mathbf{A}^{\mathrm{T}}\mathbf{A}+\lambda\mathbf{I} \right)^{-1} \mathbf{A}^{\mathrm{T}}\mathbf{\Delta T}^{\mathrm o}, \end{equation}donde $\mathbf{\Delta T}^{\mathrm o}$ son los datos observados, $\mathbf{I}$ la matriz identidad y $\lambda>0$ el parámetro de regularización.
$$\begin{pmatrix} {\B{A}} \\ \sqrt{\lambda}\B{I} \end{pmatrix} \B{m} =\begin{pmatrix} {\B{\Delta T}^{\mathrm o}} \\ \B{0} \end{pmatrix}.$$El problema anterior es equivalente a resolver por mínimos cuadrados al sistema aumentado
Muchas veces, para resolver el problema de cuadrados mínimos con regularización, empleamos programas que esperan entradas de la forma m = solve(A,DT). Podemos utilizar esta identidad para cargar la matriz y el vector aumentados correspondientes.
Resolver el sistema lineal de ecuaciones por medio de mínimos cuadrados con regularización involucra, por lo general, un costo computacional muy alto. Exiten numerosas publicaciones donde se detallan estrategias o algoritmos de inversión apropiados para manejar una gran cantidad de observaciones. A modo de ejemplo, Siqueira et al. (2017) proponen un método iterativo para el caso gravimétrico, que no requiere invertir el sistema de ecuaciones lineales para obtener $\mathbf{m}^{*}$.
Leemos el dato del archivo datos_eqlayer_magnetic.txt y graficamos adecuadamente.
data = np.loadtxt("datos_eqlayer_magnetic.txt")
x,y,z,tf = data # Norte [m], Este [m], altura de observación [m] , TFA [nT]
Número de observaciones 7054 Altura media de las observaciones (m): 541.8
Para ilustrar la metodología, seleccionamos una menor cantidad de puntos sobre toda el área prospectada.
# Selecciono "nselect" puntos del total, sobre toda el área prospectada:
nselect = 3000
tmp = np.random.choice(len(tf), nselect, replace=False)
x,y,z,tf = x[tmp],y[tmp],z[tmp],tf[tmp]
Número de observaciones 3000 Altura media de las observaciones (m): 542.4
Seleccionamos un valor para la profundidad de la capa equivalente $h$ y la altura $z^c$ a donde será interpolado el mapa luego de realizar la inversión.
h = 1000. # m
znew = 1500 # Altura de continuación (m)
Profundidad del plano equivalente = 1000.0 m Altura del mapa interpolado = 1500 m
def green_function(xp,yp,zp,xm,ym,zm):
distance = np.sqrt((xp-xm)**2 + (yp-ym)**2 + (zp-zm)**2);
return (zm-zp)/ distance**3
Cada elemento de $A$ se obtiene de:
A[i,j] = green_function(x[i],y[i],-z[i],x[j],y[j],+h); # Como h>0 (profundidad), entonces z<0 (alturas)
Cada fila de $A$ puede obtenerse directamente con
A[i,:] = green_function(x[i],y[i],-z[i],x[:],y[:],+h); # Como h>0 (profundidad), entonces z<0 (alturas)
lo que hace la construcción de $A$ mucho más rápida.
Matriz de diseño... Hecho.
Realizamos un preprocesamiento del dato previo a la inversión:
tf_mean = tf.mean()
tf_std = tf.std()
tf = (tf-tf_mean) / tf_std # normalización
Acondicionamos la matriz de diseño:
scale = A.std(axis=0)
P = np.diag(1/scale)
Obtenemos $\mathbf{m^*}$ (con preacondicionamiento):
lamb = 1.0
G = np.vstack([P @ A,lamb*I]) # matriz aumentada con A normalizada por izquierda
b = np.hstack([P @ DT,zeros] # vector aumentado con DT multiplicado por P por izquierda
m = np.linalg.lstsq(G,b,rcond=None)[0]
Si no realizamos el acondicionamiento del sistema multiplicando por izquierda con la matriz $\mathbf{P}$: $\mathbf{P}\mathbf{A}\mathbf{m}=\mathbf{P}\mathbf{\Delta T}$ , el resultado $\mathbf{m}$ no es aceptable para los datos considerados.
Parámetro de regularización: 1.0 Inversión ... Hecho.
Una forma sencilla de verificar si la profundidad de la capa equivalente $h$ propuesta es adecuada (y el parámetro de regularización seleccionado), consiste en verificar la bondad del ajuste. Para ello observamos la distribución del vector residuo $\B{r}=\B{\Delta T}^{\text{o}}-\B{\Delta T}$, entre los datos observados, $\B{\Delta T}^{\text{o}}$, y los calculados, $\B{\Delta T}=\B{A}\B{m^*}$.
tf_mod = A @ m; # dato modelado
res = tf - tf_mod; # residuos
res = res * tf_std + tf_mean # residuo en el rango de los datos observados
En esta práctica también visualizamos los datos simulados y los datos observados.
Valor medio de los residuos (nT): 26.5 Desviación estándar de los residuos (nT): 6.6 Valor medio de la TFA observada: 26.5 Desviación estándard de la TFA observada: 190.8
Generamos una cuadrícula donde interpolar la TFA:
npoints = 50
xnew = np.linspace(x.min(),x.max(),npoints)
ynew = np.linspace(y.min(),y.max(),npoints)
Interpolamos en la cuadrícula (xnew,ynew,znew) por medio de
...
for k in np.arange(len(m)):
tf_int[i,j] += m[k] * green_function(xnew[i],ynew[j],-znew,x[k],y[k],+h); # observar signos de "z" y "h"
O más rápido (utilizando el producto escalar):
...
tf_int[i,j] = m[:] @ green_function(xnew[i],ynew[j],-znew,x[:],y[:],+h); # observar signos de "z" y "h"
Número de puntos a interpolar 50 x 50 dx interpolado (m) = 961 dy interpolado (m) = 1591
Interpolando... Hecho.
Finalmente, volvemos al rango de dato original:
tf_int = tf_int * tf_std + tf_mean # regresando al rango original de las observaciones.
Finalmente, visualizamos el campo interpolado y continuado a la altura $z^c$ seleccionada.
Eso es todo por hoy.
Revisando el código que se encuentra en
https://github.com/fatiando/verde/blob/master/verde/base/least_squares.py
observamos un técnica para estabilizar el sistema de ecuaciones normales. En este código se comenta que las columnas de la matriz de diseño (matriz jacobiana) se escalan de forma que tengan varianza unitaria. Esta regularización se dice que permite utilizar un rango más amplio para ajustar el parámetro de regularización. La normalización es luego removida para el problema directo.
scale = A.std(axis=0) # std de cada columna
for j in np.arange(len(scale)):
A[:,j] /= scale[j]
...
m = m / scale
Este tema también es tratado en el libro de Max Meju (1995):
"It may also be desirable to scale each column of the matrix A by the root mean of squares value of the coefficients as it can speed up the convergence of the iterative process (Marquard, 1963). However, this should be done only when the elements of one row are markedly different from those in another row because of the round-off errors incurred during scaling".
Blakely, R.J., 1996. Potential theory in gravity and magnetic applications. Cambridge university press.
Cooper, G.R.J., 2000. Gridding gravity data using an equivalent layer. Computers & Geosciences, 26(2), pp.227-233.
Dampney, C.N.G., 1969. The equivalent source technique. Geophysics, 34(1), pp.39-53.
Meju, Max A., 1995. Geophysical Data Analysis: Understanding Inverse Problem Theory and Practice. Society of Exploration Geophysicists.
Siqueira, F.C., Oliveira Jr, V.C. and Barbosa, V.C., 2017. Fast iterative equivalent-layer technique for gravity data processing: A method grounded on excess mass constraint. Geophysics, 82(4), pp.G57-G69.
Esta notebook está basada en el siguiente ejemplo. Para lo cual se utiliza el paquete open-source harmonica, un proyecto del software de modelado e inversión en geofísica Fatiando a Terra:
Uieda, L., V. C. Oliveira Jr, and V. C. F. Barbosa (2013), Modeling the Earth with Fatiando a Terra, Proceedings of the 12th Python in Science Conference, pp. 91 - 98.