Métodos potenciales de prospección, FCAG, 2024.
Leemos el dato sintético, el cual consiste de $200\times 200$ coordenadas de observación $x_i,y_i$ (m) y los valores de la anomalía magnética total (TFA) $\Delta T_i$ (nT). El dato tiene un nivel de ruido de unos 20 nT.
x, y, TFA = np.loadtxt("mag-sintetico.dat").T # x [m], y [m], T [nT]
nx, ny = 200, 200
x = np.reshape(x, (ny,nx))
y = np.reshape(y, (ny,nx))
TFA = np.reshape(TFA, (ny,nx))
...
plt.contourf(y, x, TFA, 100, cmap="inferno")
Graficamos el mapa de anomalía total TFA en unidades de nT. La anomalía de intensidad total sintética simula un borde continental pasivo. Se aprecia una intrusión ígnea sobre la corteza oceánica (región $y>0$), la delineación del límite de la corteza oceánica y continental ($y=0$) y un dique atravezando la parte superficial de la corteza continental ($y<0$).
Calculamos los números de onda angular adimensional $\kappa_x,\kappa_y$ y los números de onda angular dimensional $k_x,k_y$ que utilizaremos a lo largo de esta clase.
kx = np.fft.fftfreq(x.shape[0])
ky = np.fft.fftfreq(y.shape[0])
kx = (2*np.pi) * kx / dx # número de onda angular dimensional, ya que divido por dx.
ky = (2*np.pi) * ky / dy # número de onda angular dimensional, ya que divido por dy.
KX,KY = np.meshgrid(kx,ky)
Kr = np.sqrt(KX*KX+KY*KY)
Trabajamos también en el dominio de las frecuencias para el cálculo de las derivadas primeras espaciales. Para ello hacemos uso de las identidades:
\begin{align} \mathcal{F}\{\frac{{\partial}\phi}{{\partial}x}\} &= ik_x\mathcal{F}\{\phi\},\\ \mathcal{F}\{\frac{{\partial}\phi}{{\partial}y}\} &= ik_y\mathcal{F}\{\phi\}. \end{align}
¿Cómo podríamos reducir el efecto observado en los bordes del resultado?
Otra forma de proceder podría ser utilizando operadores de derivada hacia adelante. Para ello, hacemos el producto en el dominio de Fourier entre la TDF 2D de la cuadrícula y la respuesta en frecuencia de los operadores de derivada espacial hacia adelante en la dirección $x$ y en la dirección $y$:
Comparamos el RMS (root mean square) de los residuos de entre las derivadas obtenidas (por definición y con diferencias hacia adelante).
RMS derivadas en x: 0.0106 RMS derivadas en y: 0.0129
Comparamos los resultados sobre una linea que corta por el medio de la cuadrícula:
Ahora utilizamos operadorres de diferencias centrales en el dominio de Fourier.
Podemos verificar los resultados de las derivadas respecto a $x$ e $y$ con el uso apropiado de la subrutina np.gradient(). Esta subrutina aplica por defecto diferencias centrales en los puntos interiores de la cuadrícula y diferencias hacia adelante o atrás en los bordes.
grad_y, grad_x = np.gradient(TFA, edge_order = 1)
grad_y = grad_y / dy
grad_x = grad_x / dx
El RMS (root mean square) de los residuos de cada derivada (por definición vs. diferencias finitas) es:
RMS derivadas en x: 0.0094 RMS derivadas en y: 0.012
El RMS entre los resultados obtenidos del procesamiento en el dominio de Fourier con el operador de diferencias centrales y los dados por la subrutina de NumPy (en puntos interiores de la cuadrícula) es:
RMS derivadas en x: 1.5e-17 RMS derivadas en y: 2.1e-17
Comparamos los resultados sobre una linea que corta por el medio de la cuadrícula:
Aquí calculamos $\nabla \cdot$TFA. Nuevamente, utilizamos el dominio de Fourier (con la respuesta en frecuencia de operadores de diferencias centrales). Comparamos con los resultados que pueden obtenerse por medio de la subrutina numpy.gradient. El cálculo de una divergencia es de particular interés también para otras ramas de la geofísica. Por ejemplo, en el procesamiento de aplanamiento de cubos sísmicos poststack.
RMS residuo (interior de la cuadrícula): 2.5e-17
¿Cómo podríamos reducir el efecto observado en los bordes del resultado obtenido por medio del dominio de Fourier? Podemos extender límites del dato original rellenando con unas pocas muestras dadas por la reflexión de las datos en los bordes de la cuadrícula. Obtenemos:
RMS residuo: 0.000989687612511392
También podemos efectuar reflexión especular (mirrorring). Obtenemos:
RMS residuo: 0.000494843806255697
Calculamos el laplaciano de la cuadrícula con la respuesta en frecuencia del filtro $3\times 3$ de diferencias finitas. En este ejemplo empleamos la cuadrícula extendida por reflexión especular para procesar en el dominio de Fourier. Comparamos con el resultado de emplear la convolución lineal 2D entre el filtro laplaciano y la cuadrícula:
from scipy.signal import convolve2d
T = ... # TFA [nT]
kernel = np.array([0,-1,0,-1,4,-1,0,-1,0])
kernel = np.reshape(kernel,(3,3));
T_lap = convolve2d(T, kernel, mode="same");
El cálculo del laplaciano juega un rol fundamental en otras aplicaciones de la geofísica. Por ejemplo, en el procesamiento de aplanamiento de cubos sísmicos poststack; full-waveform inversion (FWI); phase unwrapping, y donde sea que encontremos una ecuación de Poisson.
[[ 0 -1 0] [-1 4 -1] [ 0 -1 0]]
RMS residuo (puntos interiores de la cuadrícula): 6.0767e-14
Como observamos los resultados lucen "ruidosos". El dato original tiene ruido aleatorio de 20 nT que contamina todo el espectro del número de onda. El ruido de alto número de onda ha sido exacerbado por el operador laplaciano. El laplaciano involucra filtros de derivada segunda, que son pasa altos en el dominio transformado.
Comparamos los resultados sobre una linea que corta por el medio de la cuadrícula:
Eso es todo por hoy.