Ejemplo sencillo del filtro Butterworth¶

Métodos potenciales de prospección, FCAG, 2024.

En el dominio transformado de Fourier, la aplicación de un filtro Butterworth viene dada por

\begin{equation} \mathcal{F}[U^f(x,y,z_0)]=\mathcal{F}[b]\mathcal{F}[U(x,y,z_0)], \end{equation}

siendo

\begin{equation} \mathcal{F}[b]= \dfrac{1}{\sqrt{(1+ (k_{r}/{\color{blue}{k_c}})^{\color{blue} n}}}, \end{equation}

la transformada 2D de Fourier del filtro de paso bajo Butterworth $b(x,y,z_0;{\color{blue}n},{\color{blue}{k_c}})$ con 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}$, la norma del vector número de onda.

El filtro de Butterworth es muy utilizado en los métodos potenciales de prospección. Aplicar este filtro consiste en definir un orden $n$ y número de onda radial de corte $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.

El operador de paso bajo Butterworth¶

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       = ...
ky       = ...
kr       = ...
krc      = np.pi/2 # número de onda angular adimensional de corte
n        = 8       # orden del filtro
filtro   = 1.0/np.sqrt( 1. + (kr/krc)**n )

Aplicación: dato sintético¶

$\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 aquí al 5% de la amplitud máxima de la anomalía.
ruido    = np.random.normal(loc=0,scale=sigma,size=anomalia.shape)
gz mínimo   [mGal]:  -0.58
gz máximo   [mGal]:  0.57
gz promedio [mGal]:  -0.0

Procesamiento¶

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.

  • Si utilizamos número de onda adimensionales, calculamos el número de onda radial adimensional $k^a_r$ y probamos el filtro para números de onda adimensional de corte $\kappa_c$ entre 0 y $\pi$. Por lo general, podemos comenzar por $\kappa_c=\pi/2$ y reducimos el valor según el análisis del resultado obtenido y del residuo.
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)
  • De utilizar números de onda (con unidades), entonces calculamos el número de onda radial $k_r$ y probamos el filtro para números de onda de corte $k_c$ entre 0 y $\text{min}\{k^\text{Nyquist}_x,k^\text{Nyquist}_y\}$. En este caso, podriamos comenzar por $k_c=\text{min}\{k_x^\text{Nyquist},k_y^\text{Nyquist}\}/2$ y reducir el valor analizando el resultado obtenido y su residuo.
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.52
gz promedio [mGal]:  -0.0

Resultados¶

Representamos gráficamente los resultados e interpretamos.

¿Se observan en el residuo amplitudes importantes vinculadas a la anomalía?

Eso es todo por hoy.

Referencias¶

  • Advances in Gravity and Magnetic Processing and Interpretation, J. Derek Fairhead, EAGE 2015
  • Manual de aplicaciones de O@sis Mont@j