Métodos potenciales de prospección, FCAG, 2024.
Cargamos 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 posee además ruido aditivo de amplitud 20 nT.
x, y, T = 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))
T = np.reshape(T, (ny,nx))
...
dx = (x.max()-x.min())/(nx-1) # intervalo de muestreo en X.
dy = (y.max()-y.min())/(ny-1) # intervalo de muestreo en Y.
...
print(f"dx ={dx} [m], dy = {dy} [m]") # imprimo por pantalla parámetros de interés.
...
plt.contourf(y, x, T, 100, cmap="inferno") # visualización.
dx = 502.51 [m], dy = 502.51 [m]
Graficamos el mapa de anomalía total $\Delta T$.
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 atravesando la parte superficial de la corteza continental ($y<0$).
Procesamos la cuadrícula de anomalías en el dominio de Fourier para calcular el mapa de derivada primera vertical $\frac{{\partial}\Delta T}{{\partial}z}$ y derivada segunda vertical $\frac{{\partial}^2\Delta T}{{\partial}z^2}$. Para ello, utilizamos las identidades $$\mathcal{F}\{\dfrac{{\partial}^{n}\phi}{{\partial}z^n}\}=|k|^{n}\mathcal{F}\{\phi\},$$ para $n=1,2$.
Estas propiedades se cumplen para un campo potencial $\phi$. La derivada primera vertical puede obtenerse de considerarse la continuación ascendente. La derivada segunda vertical puede deducirse como consecuencia de ser $\phi$ un potencial armónico: $\Delta \phi = 0$.
Para realizar el procesamiento, pueden realizarse los siguientes pasos:
kx = ... # número de onda en X dimensional [1/m]
ky = ... # número de onda en Y dimensional [1/m]
kr = ... # número radial de onda dimensional [1/m]
order = 1
T_F = np.fft.fft2(T)
dTdz_F = T_F * kr ** order
dTdz = np.real(np.fft.ifft2(dTdz_F)) # derivada primera vertical [nT/m]
...
dTdz_2 = ... # derivada segunda vertical [nT/m2]
Visualizamos los resultados.
El mapa ${\partial_z}{\Delta}T$ presenta una mejor delineación de los cuerpos. También observamos que la delineación de ${\partial^2_z}{\Delta}T$ es más localizada, pero el mapa es mucho más sensible al ruido presente en $\Delta T$. Por ejemplo, el dique que atraviesa la corteza continental queda enmascarado por el ruido en el mapa ${\partial^2_z}{\Delta}T$.
Calculamos el atributo de gradiente total, magnitud del gradiente o magnitud,
$$\text{gt}(x,y)=\sqrt{{\partial_x}{\Delta}T^2+{\partial_y}{\Delta}T^2+{\partial_z}{\Delta}T^2},$$operando en el dominio de Fourier. Para ello aplicamos las relaciones
\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\}; \\ \mathcal{F}\{\frac{{\partial}\phi}{{\partial}z}\} &= |k|\mathcal{F}\{\phi\}. \end{align}
¿Cómo podríamos reducir el efecto observado en los bordes del resultado?
Utilizamos padding con la técnica de inverse reflect o derivate mirror (Smith et al., 2022) para obtener continuidad en las derivadas primeras espaciales sobre los bordes de la cuadrícula original y reducir el rizado observado en las aplicaciones anteriores.
Nota. Si bien aquí solo aplicamos el padding, es recomendable también aplicar un taper al dato luego del padding.
Observamos el padding reflect inverse en una señal 1D y en el dato original.
True
Implementamos un simple taper a partir de una ventana triangular.
Eso es todo por hoy.