Métodos potenciales de prospección, FCAG, 2024.
En el dominio transformado de Fourier, la continuación analítica de una función armónica $U(x,y,z)$ observada sobre una superficie horizontal $z=z_0$ a $U(x,y,z_0+\color{blue}{\Delta z})$ viene dada por
\begin{equation} \mathcal{F}[U(x,y,z_0+\color{blue}{\Delta z})]=\mathcal{F}[\psi_c]\;\mathcal{F}[U(x,y,z_0)], \end{equation}
con
\begin{equation} \mathcal{F}[\psi_c]=e^{-\color{blue}{\Delta z}||\mathbf{k}||}, \end{equation}
la transformada bidimensional de Fourier del operador de continuación analítica y $||\mathbf{k}||=\sqrt{k_x^2+k_y^2}$, la norma del vector número de onda.
Realizar la continuación analítica consiste principalmente en calcular la TDF 2D del dato observado ($\mathcal{F}[U]$) y la expresión del operador de continuación en el dominio del número de onda ($\mathcal{F}[\psi_c]$), hacer el producto y antitransformar.
En los métodos potenciales la función armónica puede ser la componente vertical de la atracción gravitatoria $g_z(x,y,z_0)$ o la anomalía magnética escalar de intensidad total (TFA) $\Delta T(x,y,z_0)$.
Graficamos el espectro de amplitud del operador de continuación ascendente $\mathcal{F}[\psi_r]=e^{-\Delta h ||\mathbf{k}||}$, $\Delta h>0$. Observamos que el filtro toma valores entre $(0,1]$ debido a que estamos tratando con un exponencial decreciente. El valor máximo se alcanza para $||\mathbf{k}||=0.$
El espectro de amplitud del filtro podemos obtenerlo por medio del siguiente código,
# Números de onda
kx = 2*np.pi*np.fft.fftfreq(nx,dx) # número de onda en x [1/m]
ky = 2*np.pi*np.fft.fftfreq(ny,dy) # número de onda en y [1/m]
kkx, kky = np.meshgrid(kx, ky) # cuadrícula
kr = np.sqrt(kkx*kkx +kky*kky) # número de onda radial [1/m] (cuadrícula)
# Filtro
h = ... # altura de continuación [m].
filtro = np.exp(-h*kr) # TDF 2D del filtro de continuación ascendente (si h>0).
Graficamos el espectro de amplitud del operador de continuación descendente $\mathcal{F}[\psi_r]=e^{-\Delta h ||\mathbf{k}||}$ con , $\Delta h<0$. Observamos ahora que el filtro toma valores en $[1,+\infty)$ debido a que estamos tratando con un exponencial de argumento positivo. Es decir, toma valores grandes según $||\mathbf{k}||$ aumenta. El valor mínimo se alcanza para $||\mathbf{k}||=0.$
$\newcommand{\vec}[1]{\mathbf{#1}}$ $\newcommand{\versor}[1]{\hat{\mathbf{#1}}}$ $\newcommand{den}{\Delta\sigma}$ $\newcommand{G}{\mathcal{G}}$ $\newcommand{M}{\dfrac{4}{3}\pi\den R^3}$ $\newcommand{cte}{\dfrac{4}{3}\pi\G\den R^3}$
La componente vertical de la atracción gravitatoria debida a una esfera con contraste de densidad $\Delta\sigma$ viene dada por:
\begin{equation} g_z(\vec{r},\vec{r}_0) = \vec{g}\cdot\versor{k} = \cte\frac{z-z_0}{\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}^{\phantom{.}3}}. \end{equation}donde $\mathbf{r}=(x,y,z)$ son las coordenadas del punto de observación, $\mathbf{r}_0=(x_0,y_0,z_0)$ las coordenadas del centro de la esfera, $R$ el radio de la esfera y $r^2=||\vec{r}-\vec{r}_0||_2^2=(x-x_0)^2+(y-y_0)^2+(z-z_0)^2$.
Generamos un dato sintético de dos cuerpos esféricos con cierto contraste de densidad. Para ello, escribimos una subrutina gz_esfera() para el cálculo de $g_z$.
def gz_esfera(P,rho,R,P0):
G = 6.67e-11 # [m3 / kg s2]
cte = (4./3.)*np.pi*G
convert = 1e5 # m/s2 to mGal
x,y,z = P
x0,y0,z0 = P0
d2 = (x-x0)**2+(y-y0)**2+(z-z0)**2
g = cte * rho * (R**3) * (z-z0)/(np.sqrt(d2)**3)
return g * convert
Tomamos parámetros para generar un modelo de dos esferas y elegimos parámetros de campaña adecuados.
constraste [kg/m3] : 500.0 radio [m] : 200.0 profundidad [m] : -400.0 posición de los centros [m]: 240.0 -240.0
Calculamos la anomalía y luego visualizamos,
g1 = gz_esfera(...) # anomalía del cuerpo 1 [mGal]
g2 = gz_esfera(...) # anomalía del cuerpo 2 [mGal]
anomalia = g1 + g2 # mGal
gz mínimo [mGal]: 0.02 gz máximo [mGal]: 0.91 gz promedio [mGal]: 0.19
Aplicamos la continuación analítica al dato sintético recién calculado. Podemos tomar $\Delta z>0$ para hacer la continuación ascendente. También podemos considerar $\Delta z<0$ para realizar la continuación descendente, con recaudos.
kr = ...
h = ...
A = np.fft.fft2(anomalia) # TDF 2D del dato
filtro = np.exp(-h*kr) # TDF 2D del filtro
cont = np.real( np.fft.ifft2(A*filtro) ) # continuación
continuación para dz [m] : 100.0
gz mínimo [mGal]: 0.06 gz máximo [mGal]: 0.67 gz promedio [mGal]: 0.19
A fines de comparar el resultado de aplicar el filtro, calculamos la respuesta analítica para $z=z_0+\Delta z$:
gz mínimo [mGal]: 0.03 gz máximo [mGal]: 0.66 gz promedio [mGal]: 0.17
Representamos gráficamente los resultados e interpretamos.
Por último, analizamos el error entre la continuación analítica $\mathbf{g^c_z}$ para distintos $z$ y el resultado analítico dado por la expresión $\mathbf{g_z}$ de las esferas.
Observamos que la continuación descendente es inestable, mientras que la continuación ascendente es estable y arroja resultados con errores aceptables.
Eso es todo por hoy.