2.8 Transformadas espectrales

Introducción¶

En las geociencias, buena parte de los datos con los que tratamos — ondas sísmicas, mareas oceánicas, variaciones de presión atmosférica — llega en forma de series de tiempo. Estas señales suelen contener información en múltiples frecuencias, y un análisis simple en el dominio del tiempo puede ocultar rasgos importantes de los datos.

Los espectrogramas son herramientas útiles para analizar y visualizar el contenido de frecuencia de las series de tiempo. Un espectrograma representa la densidad espectral de una señal a medida que cambia en el tiempo. Es especialmente útil en aplicaciones geocientíficas donde las señales son no estacionarias, es decir, su contenido de frecuencia cambia con el tiempo. Algunos ejemplos:

• Sismología de terremotos: donde las señales de baja y de alta frecuencia son relevantes en distintas etapas de un evento.
• Ciencias de la atmósfera: donde patrones periódicos como las mareas y las ondas tienen componentes de frecuencia distintivas.
• Percepción remota: donde distintos procesos, como las fluctuaciones de humedad del suelo, pueden exhibir frecuencias características a lo largo del tiempo.

Por qué importan los espectrogramas

• Análisis tiempo-frecuencia: los espectrogramas muestran cómo evoluciona en el tiempo el contenido de frecuencia de una señal, lo que los hace adecuados para estudiar datos no estacionarios.
• Extracción de características: los espectrogramas resaltan los eventos transitorios y los patrones de larga duración, útiles para detectar y caracterizar fenómenos geofísicos como sismos, deslizamientos de tierra u ondas atmosféricas.
• Análisis multiescala: muchos procesos naturales operan en distintas escalas de tiempo. Los espectrogramas permiten visualizar y extraer características a través de esas escalas.
• Visualización: ofrecen una visualización compacta e interpretable de series de tiempo complejas.

En esta sección transformamos los datos proyectándolos sobre una base de funciones. Las dos transformadas más usadas son la de Fourier y la de ondículas (wavelets).

Advertencia. Filtrar cualquier dato debe hacerse con cuidado. Los artefactos de filtrado pueden llevar a una malinterpretación completa. Trampas comunes:

• No todos los filtros preservan la causalidad: algunas señales pueden aparecer antes de los eventos y malinterpretarse como precursores.
• Filtrar sobre huecos en los datos introduce artefactos de alta frecuencia.
• Los efectos de borde al filtrar series de tiempo son difíciles de mitigar.

La lección cubre varios niveles y métodos para transformar los datos.

1. Transformadas de Fourier: 1D [Nivel 1]
2. Transformadas de Fourier: 2D [Nivel 3]
3. Espectrogramas [Nivel 2]
4. Transformadas de ondículas [Nivel 3]

🖥️ Diapositivas — Sesión 07 (mié 14 oct)

Primero descargamos los datos: sismogramas registrados en el Puget Sound para el sismo M8.2 de Chignik, Alaska, del 29 de julio de 2021. Consultamos el centro de datos FDSN con el cliente IRIS. Note que los servicios de datos de IRIS ahora los opera EarthScope, pero el nombre de cliente IRIS sigue funcionando.

El código de canal HHZ denota un sismómetro de banda ancha de alta ganancia (H), muestreado a 100 muestras por segundo (H), componente vertical (Z). Es un canal de banda ancha: registra el movimiento del suelo en una banda de frecuencias amplia, no en una sola frecuencia.

También descargamos los metadatos de la estación (la respuesta instrumental) y removemos la respuesta con remove_response(output="VEL"). Esto convierte las cuentas crudas del digitalizador en velocidad del suelo en m/s, de modo que las amplitudes tienen unidades físicas.

1. Transformadas de Fourier [Nivel 1]¶

Usamos el módulo scipy.fft para transformar las dos series de tiempo (sismo y ruido). El módulo más antiguo scipy.fftpack es heredado y ya no debe usarse.

La transformada de Fourier es una descomposición de la serie de tiempo sobre una base ortonormal de funciones coseno y seno. La transformada de Fourier de una serie de tiempo $f(t)$ (igual si la variable es el espacio $x$) es:

$\hat{F}(f) = \int_{-\infty}^\infty f(t) \exp^{-i2\pi ft} dt$

$\hat{F}(f)$ es el valor complejo de Fourier en la frecuencia $f$. La transformada de Fourier determina qué frecuencias dominan la serie de tiempo.

1.1 Nyquist¶

La transformada de Fourier que usamos en esta clase toma una serie de tiempo discreta de números reales. Si la serie abarca $T$ segundos con $N$ muestras regularmente espaciadas, el intervalo de muestreo es $dt = T/N$. La frecuencia más alta que puede resolverse en una serie de tiempo discreta, llamada frecuencia de Nyquist, está limitada por $dt$:

$F_{Nyq} = \frac{1}{2 dt}$

En la práctica, no se pueden restringir señales que varían más rápido que dos muestras de tiempo. Aquí $dt = 0.01$ s, así que $F_{Nyq} = 50$ Hz.

1.2 Incertidumbres¶

• La transformada de Fourier discreta produce una aproximación de la TF. Cuanto más corta la serie de tiempo, menos exacta la TF. Esto significa que la TF sobre ventanas de tiempo cortas es menos exacta.
• La TF supone (y requiere) la periodicidad de la serie, es decir, que la serie de tiempo finita/recortada se repetiría en el tiempo. Para forzarlo, aplicamos un taper (afinamiento gradual) a la serie de tiempo, de modo que el primer y el último punto sean iguales (a cero).

Consulte la documentación de ObsPy para conocer la función taper. Grafique los espectros de amplitud y de fase.

Notará arriba que los valores de fase están distribuidos aleatoriamente entre $-\pi$ y π. Podemos comprobarlo mostrando la distribución de los espectros de fase y de amplitud.

También podemos analizar las características espectrales de la serie de tiempo de ruido. Abajo:

1. calcule la transformada de Fourier
2. grafique los espectros de fase y de amplitud
3. grafique la distribución de los valores de fase y de amplitud

Superponga sus PDF (funciones de densidad de probabilidad).

Nota que sus diferencias estadísticas están en las colas de las distribuciones. Por lo tanto, métricas estadísticas como la mean (media) o la variance (varianza) pueden no ser discriminantes, pero la kurtosis (curtosis) podría serlo.

2. Transformadas de Fourier 2D [Nivel 3]¶

La transformada de Fourier 2D se aplica a una matriz 2D. Primero aplica una transformada de Fourier 1D a cada fila de la matriz, y luego aplica una transformada de Fourier 1D a cada columna de la matriz intermedia.

Las transformadas de Fourier 2D dan los coeficientes de Fourier que dominan una imagen. Esto puede usarse para filtrar los datos. Otra aplicación es comprimir los datos conservando unos pocos coeficientes en lugar de almacenar la imagen completa.

Practicaremos sobre un campo sintético con aspecto de topografía. La topografía real tiene un espectro de número de onda «rojo»: la mayor parte de la potencia se concentra en las longitudes de onda largas (las cordilleras — piense en la Sierra Madre Occidental o en los Andes), con cada vez menos potencia en las longitudes de onda cortas (la rugosidad de pequeña escala). La amplitud espectral decae aproximadamente como una ley de potencias del número de onda $k$, $|\hat{F}(k)| \propto k^{-\beta}$. Podemos construir un terreno fractal así directamente en el dominio de Fourier: extraer fases aleatorias, escalar las amplitudes por $k^{-\beta}$ y aplicar la transformada inversa.

Considere elevation como un conjunto de datos 2D. Podemos realizar una transformada 2D, que da un espectro en las dimensiones espaciales. Los ejes de la imagen transformada son números de onda (ciclos por km). Usamos fftshift para colocar el número de onda cero en el centro de la imagen.

La energía se concentra cerca del centro del gráfico (números de onda bajos, longitudes de onda largas), tal como se diseñó. Ahora comprimiremos la imagen conservando solo los coeficientes de Fourier más grandes.

Ahora comparamos el conjunto de datos 2D original con los datos comprimidos por Fourier. Conservar el 1 % de los coeficientes es una razón de compresión de 100:1, y la reconstrucción aún captura la estructura de gran escala del terreno. El precio es la pérdida de la rugosidad de pequeña escala, que vive en los coeficientes de número de onda alto descartados.

3. Espectrogramas [Nivel 2]¶

En problemas dependientes del tiempo y multiescala, puede ser interesante extraer características de los datos a partir de la transformada de Fourier de tiempo corto (STFT, short time Fourier transform).

La STFT es una transformada de Fourier aplicada a ventanas cortas (traslapadas) para resolver las frecuencias en distintos momentos de la serie.

El espectrograma Zxx es una transformada de los datos originales. Es común usar espectrogramas como entrada de redes neuronales, en forma de arreglos 2D.

4. Transformada continua de ondículas [Nivel 3]¶

Introducción¶

La transformada continua de ondículas (CWT, Continuous Wavelet Transform) es una herramienta importante en el análisis de datos geocientíficos, en particular para el análisis tiempo-frecuencia de señales no estacionarias. Como el espectrograma, la CWT permite ver cómo varía en el tiempo el contenido de frecuencia de una señal. Sin embargo, la CWT ofrece mejor resolución en las distintas frecuencias, lo que la hace más adecuada para analizar señales con cambios de frecuencia transitorios o localizados.

La transformada de ondículas descompone una señal en versiones escaladas y desplazadas de una función oscilante pequeña conocida como la ondícula (wavelet). Esto hace que la CWT sea idónea para las aplicaciones geocientíficas, donde muchos fenómenos, como los sismos, las erupciones volcánicas y los patrones meteorológicos, pueden manifestarse en distintas escalas y frecuencias.

Figura: funciones de base de Fourier y de ondículas. Imagen de este artículo.

Existen muchas familias canónicas de ondículas. La diferencia entre familias es típicamente su forma, su compacidad y su suavidad. Por lo general, se elige una familia para la serie de tiempo específica. Las ondículas tienen energía finita y media cero.

Figura: familias de funciones de base de ondículas. Imagen de este artículo

La transformada de ondículas es:

$\hat{F}(a,b) = \frac{1}{\sqrt{a}} \int_{-\infty}^\infty f(t) \bar{\Psi} (\frac{t-b}{a}) dt$

donde $\bar{\Psi}$ es la ondícula madre escalada por un factor $a$ y trasladada/desplazada por $b$. En la transformada continua, $a$ y $b$ toman valores continuos. La transformada discreta de ondículas es la transformada de ondículas realizada sobre un número finito de escalas y desplazamientos.

La representación tiempo-escala de una serie de tiempo es un escalograma. Las escalas pueden convertirse en pseudofrecuencias: si $f_c$ es la frecuencia central de la ondícula, la escala es $a$ y el intervalo de muestreo es $dt$, entonces la pseudofrecuencia es $f_a = f_c/(a \, dt)$.

Por qué importan las transformadas continuas de ondículas¶

• Análisis multirresolución: la CWT captura tanto las tendencias de baja frecuencia y larga duración como los rasgos de alta frecuencia y corta duración. Esto es particularmente valioso en las geociencias, donde los procesos ocurren en distintas escalas temporales y espaciales.
• Datos no estacionarios: muchas señales geocientíficas son no estacionarias, es decir, sus propiedades estadísticas cambian con el tiempo. La CWT revela esas componentes de frecuencia variables en el tiempo.
• Detección de rasgos locales: la CWT es idónea para identificar eventos localizados, como ondas sísmicas, deslizamientos de tierra o perturbaciones atmosféricas, al analizar cómo cambia la señal tanto en tiempo como en frecuencia.
• Mejor resolución tiempo-frecuencia: la transformada de ondículas ofrece una localización en tiempo y frecuencia más precisa que métodos como la transformada de Fourier o el espectrograma, sobre todo para eventos de corta vida.

Uso en sismología¶

En los datos sísmicos, los distintos tipos de ondas sísmicas (ondas P, ondas S y ondas superficiales) ocurren en frecuencias y duraciones distintas. Al aplicar la CWT, los sismólogos pueden detectar esas distintas fases de onda y sus tiempos de llegada precisos, que importan para la caracterización de los sismos.

Ejemplo en Python: transformada continua de ondículas de datos sísmicos¶

Usamos el paquete PyWavelets (pywt). Las funciones más antiguas scipy.signal.cwt y scipy.signal.morlet2 fueron removidas de SciPy, así que pywt es ahora la herramienta estándar. Elegimos la ondícula de Morlet compleja 'cmor1.5-1.0' (ancho de banda 1.5, frecuencia central 1.0), una elección común para datos sísmicos porque su forma compleja da tanto la amplitud como la fase.

pywt.cwt trabaja en escalas. Elegimos las frecuencias que queremos resolver (0.1 a 40 Hz, espaciadas logarítmicamente) y luego las convertimos en escalas con la relación $a = f_c / (f \, dt)$, que es lo que calcula pywt.frequency2scale. Para mantener liviano el cálculo, aplicamos la CWT a los primeros 20 minutos del registro del sismo, que contienen las ondas P, S y superficiales.

El escalograma muestra la onda P llegando primero con energía de alta frecuencia, seguida por la onda S y las ondas superficiales de periodo largo por debajo de 0.1-0.5 Hz. Compare esto con el espectrograma STFT de arriba: la transformada de ondículas resuelve mejor las ondas superficiales de baja frecuencia porque su ventana de tiempo se adapta a cada frecuencia.

Ventajas de la CWT en las geociencias:¶

• Detección de eventos localizados en el tiempo: la CWT identifica eventos geofísicos de corta vida, como la llegada de las ondas sísmicas durante un sismo.
• Fenómenos multiescala: procesos naturales como la actividad tectónica y las mareas oceánicas operan en un rango amplio de escalas temporales y espaciales. La CWT puede analizar señales de esos procesos usando un rango amplio de escalas.
• Detección de bordes: la transformada de ondículas es particularmente buena para detectar cambios o discontinuidades en las señales, como el inicio abrupto de las ondas sísmicas o las fronteras entre distintas capas geofísicas.

Casos de uso en las geociencias:¶

• Alerta sísmica temprana: al aplicar la CWT, los geocientíficos pueden detectar con más exactitud la primera llegada de las ondas sísmicas, que importa para los sistemas de alerta temprana como el SASMEX mexicano.
• Detección de deslizamientos de tierra: la CWT ayuda a identificar las señales de alta frecuencia indicativas de deslizamientos, ya que estos eventos suelen caracterizarse por ráfagas cortas de energía.
• Tremor volcánico: la CWT puede usarse para analizar las señales de tremor volcánico — como las que registran las redes de monitoreo del Popocatépetl —, que suelen ser complejas y exhibir rasgos tanto de corta como de larga duración.

Conclusión¶

La transformada continua de ondículas (CWT) es una herramienta útil para el análisis tiempo-frecuencia en las geociencias. Su capacidad de resolver señales a través de distintas escalas y frecuencias la hace idónea para estudiar procesos no estacionarios. Bibliotecas de Python como pywt, obspy y matplotlib permiten a los geocientíficos aplicar la CWT a sus datos y extraer conocimiento de fenómenos geofísicos complejos.

Las transformadas tiempo-frecuencia consumen tiempo de cómputo en un flujo de trabajo. Comparemos el costo de la CWT y de la STFT sobre el mismo segmento de 20 minutos.

La STFT es mucho más barata que la CWT. De estas transformadas podemos extraer características estadísticas similares.