Métodos potenciales de prospección, FCAG, 2024.
En el dominio transformado de Fourier, la aplicación de un filtro Butterworth viene dada por
En la expresión anterior,
es la transformada 2D de Fourier del filtro de paso bajo Butterworth $b(x,y,z_0;{\color{blue}n},{\color{blue}{k_c}})$ de orden $\color{blue}n$ y número de onda radial de corte $\color{blue}{k_{c}}$, siendo $k_r=||\mathbf{k}||=\sqrt{k_x^2+k_y^2}$ es el número de onda radial.
El filtro de Butterworth es muy utilizado en los métodos potenciales de prospección. Aplicar este filtro consiste en definir un orden ${\color{blue}n}$ y número de onda radial de corte $\color{blue}{k_{c}}$ adecuados (muchas veces por prueba y error), calcular la transformada discreta 2D de Fourier de la cuadrícula de observaciones ($\mathcal{F}[U]$), calcular $\mathcal{F}[b]$, hacer el producto y antitransformar.
Graficamos el espectro de amplitud del operador de paso bajo Butterworth $\mathcal{F}[b]$. A los fines ilustrativos, consideramos $n=8$ y ubicamos el número de onda angular adimensional de corte en $\kappa_c = \pi/2$.
kx = ... # número de onda angular adimensional en x
ky = ... # número de onda angular adimensional en y
kr = ... # número de onda radial adimensional
kc = np.pi/2 # número de onda angular adimensional de corte
n = 8 # orden del filtro Butterworth
filtro = 1.0/np.sqrt( 1. + (kr/kc)**n ) # TDF 2D del filtro Butterworth
$\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 para el cálculo de $g_z$. Luego escribimos una pequeña subrutina para conocer los rangos de la anomalía generada, de forma de controlar resultados y comparar.
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.
Nota. Para hacer evidente el efecto del filtro, sumamos ruido aleatorio a la anomalía. El ruido lo modelamos como una variable aleatoria $r$ que sigue una distribución Gaussiana de media $\mu=0$ y desviación estándard $\sigma$: $r \sim \mathcal{N}(\mu,\sigma)$. Para este ejemplo, ajustamos $\sigma$ a cierto porcentaje de la amplitud máxima de la anomalía observada.
mu, sigma = 0.,(5./100)*np.max(np.abs(anomalia)) # la desviación se ajusta al 5% de la amplitud máxima de la anomalía.
ruido = np.random.normal(loc=0,scale=sigma,size=anomalia.shape)
Una alternativa para generar ruido más realista es aplicar al ruido un pasa bajas o pasa bandas para evitar la situación de ruido blanco dada por
np.random.normal(). Nos decidimos por esta opción aquí.
gz mínimo [mGal]: -0.53 gz máximo [mGal]: 0.54 gz promedio [mGal]: -0.0
Elegimos valores de $n$ y $k_{rc}$ y aplicamos el filtro. Utilizamos $n=8$ y $\kappa_c=\pi/4$ en este ejemplo en particular. Otros parámetros son posibles.
nx,ny = anomalia.shape
kx = 2*np.pi*np.fft.fftfreq(nx,1) # adimensional, dx == 1
ky = 2*np.pi*np.fft.fftfreq(ny,1) # adimensional, dy == 1
kkx, kky = np.meshgrid(kx, ky)
# módulo del vector numero de onda
kr = np.sqrt(kkx*kkx +kky*kky)
Filtro Butterworth de paso bajo: orden 8 y corte en k_c angular adimensional de 0.79
gz mínimo [mGal]: -0.53 gz máximo [mGal]: 0.53 gz promedio [mGal]: -0.0
Representamos gráficamente los resultados e interpretamos. La anomalía residual la calculamos como la diferencia entre la anomalía observada y la anomalía filtrada.
¿Se observan en el residuo amplitudes importantes vinculadas a la anomalía?
Eso es todo por hoy.