Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Generación de ruido sintético y datos sintéticos para las geociencias

1. Introducción

El ruido en los datos geocientíficos proviene del entorno, del instrumento y de la actividad humana. Si podemos simular el ruido, podemos construir conjuntos de datos controlados con una verdad de referencia conocida. Eso permite probar filtros, evaluar detectores contra puntos de comparación y aumentar los conjuntos de entrenamiento para el aprendizaje automático.

Esta lección cubre el ruido sintético en dos contextos:

  1. Series de tiempo 1D: registros sísmicos, presión atmosférica, temperatura superficial del mar.
  2. Campos geoespaciales 2D: topografía, humedad del suelo, mapas de temperatura.

Terminamos con un ejemplo sismológico realista y una discusión sobre cuándo los datos sintéticos son admisibles en la investigación.

Tipos de ruido

  • Ruido blanco: densidad espectral de potencia constante en todas las frecuencias. Típico del ruido instrumental.
  • Ruido rosa: densidad espectral de potencia proporcional a 1/f1/f. Común en sistemas naturales como los registros climáticos.
  • Ruido gaussiano: ruido cuya distribución de amplitudes es normal. El ruido blanco suele extraerse de una gaussiana.
  • Ruido espacialmente correlacionado: ruido que varía suavemente en el espacio, como los artefactos atmosféricos o de sensor en los mapas.

A lo largo de este cuaderno usamos el generador aleatorio moderno de numpy, con semilla fija por reproducibilidad. La semilla importa: un conjunto de datos sintético que no puede regenerarse no es reproducible.

🖥️ Diapositivas — Sesión 09 (lun 19 oct)

2. Modelos de ruido en 1D y 2D

2.1 Ruido blanco y ruido rosa (1D)

El ruido blanco es una secuencia de extracciones gaussianas independientes. El ruido rosa se construye en el dominio de la frecuencia: imponemos un espectro de amplitud proporcional a 1/f1/\sqrt{f} (de modo que la potencia vaya como 1/f1/f), sorteamos fases aleatorias y transformamos de regreso al dominio del tiempo. Usamos np.fft.irfft, que toma el espectro solo en las frecuencias no negativas y devuelve una señal real por construcción.

<Figure size 1000x600 with 2 Axes>

El ruido blanco fluctúa rápidamente, sin memoria de una muestra a la siguiente. El ruido rosa deambula: las frecuencias bajas concentran la mayor parte de la potencia, así que la serie deriva en escalas de tiempo largas.

Aplicaciones. El ruido blanco simula el autorruido instrumental en los registros sísmicos o atmosféricos. El ruido rosa imita la variabilidad natural, dominada por las frecuencias bajas en muchos procesos geofísicos.

2.2 Ruido espacialmente correlacionado (2D)

El ruido geoespacial rara vez es independiente de un píxel al siguiente. La interferencia atmosférica y la deriva de los sensores varían suavemente en el espacio. Una manera simple de simularlo: generar un campo 2D de ruido blanco y luego suavizarlo con un filtro gaussiano. El ancho del filtro fija la longitud de correlación.

<Figure size 1200x600 with 4 Axes>

Aplicaciones. El ruido espacialmente correlacionado simula la distorsión atmosférica en las imágenes satelitales, los artefactos de interpolación en los productos en malla y los errores suaves de sensor en los mapas geofísicos.

3. Eventos transitorios sintéticos

Ahora construimos una serie de tiempo que contiene un evento transitorio más ruido. Este es el ingrediente básico de un problema de clasificación binaria: ¿contiene una ventana señal o solo ruido?

3.1 La señal del evento: la ondícula de Ricker

La ondícula (wavelet) de Ricker, también llamada ondícula de sombrero mexicano, es un modelo estándar de una fuente sísmica impulsiva. Es la segunda derivada negativa de una gaussiana. scipy.signal.ricker fue eliminada de las versiones recientes de scipy, así que la definimos nosotros mismos.

<Figure size 640x480 with 1 Axes>

La ondícula de Ricker es suave y de banda limitada. Grafiquemos el valor absoluto de su espectro de amplitud de Fourier.

<Figure size 1000x500 with 2 Axes>

¿Cómo se ve la distribución de los datos del evento?

<Figure size 640x480 with 1 Axes>

3.2 Ruido gaussiano y relación señal-ruido

Creamos una señal pura. Ahora creamos una serie de tiempo de ruido para sumársela.

Una precaución: np.random.uniform(0, 1) sortea valores entre 0 y 1, cuya media es 0.5. Sumar eso a una señal introduce un desplazamiento de nivel medio (DC) — un pico espurio en la frecuencia cero del espectro. Usamos en cambio rng.standard_normal, que tiene media cero por construcción.

Una sola definición, usada en todas partes. Definimos la relación señal-ruido (SNR) como la amplitud absoluta pico de la señal dividida por la desviación estándar del ruido. Esta es la definición incorporada en mlgeo_synth y la que usan los experimentos con detectores del capítulo 4.3. Mantener una sola definición importa: una afirmación como «el detector falla por debajo de SNR 1» no significa nada si dos capítulos miden la SNR de maneras distintas. En el código, normalizamos la señal a amplitud pico unitaria y el ruido a desviación estándar unitaria, de modo que s + noise / SNR tenga exactamente la SNR declarada.

<Figure size 640x480 with 1 Axes>

Compare el espectro de amplitud de Fourier del ruido contra el espectro de la señal.

<Figure size 640x480 with 1 Axes>

Se ven muy distintos en el dominio espectral. El ruido reparte su energía en todas las frecuencias; la señal concentra la suya en una banda angosta.

Ahora sumamos el ruido a la señal, escalado por una relación señal-ruido (SNR). Definimos aquí la SNR como el cociente entre la amplitud absoluta máxima de la señal y la del ruido. Primero normalizamos ambas y luego dividimos el ruido por la SNR.

<Figure size 640x480 with 1 Axes>

3.3 Ruido con un color elegido: aleatorizar la fase

El ruido puede tener distinto contenido de frecuencia, o color. Una receta general construye una serie de tiempo de ruido a partir de un espectro de amplitud de Fourier elegido:

  1. Elija la amplitud en cada frecuencia no negativa (plana para el ruido blanco, un decaimiento 1/f1/f para el ruido coloreado, o el espectro de datos reales).
  2. Sortee una fase aleatoria en cada frecuencia, uniforme entre 0 y 2π2\pi.
  3. Invierta con np.fft.irfft.

Una serie de tiempo de valores reales requiere un espectro con simetría hermitiana: el valor en −f-f debe ser el conjugado complejo del valor en +f+f. Programar esa simetría a mano es propenso a errores. irfft toma solo las frecuencias no negativas y aplica la simetría por construcción, así que la usamos. Los intervalos (bins) de la frecuencia cero y de Nyquist deben ser reales, así que fijamos sus fases en cero.

<Figure size 640x480 with 1 Axes>

Sume el ruido nuevo a la señal (la ondícula de Ricker), esta vez con una SNR mucho menor, y grafique en los dominios del tiempo y de la frecuencia.

<Figure size 1000x500 with 2 Axes>

Compare las distribuciones de los datos de la señal pura, del ruido y de su suma.

<Figure size 1000x500 with 2 Axes>

3.4 Momentos estadísticos

Calcule los momentos estadísticos de la señal limpia y del ruido. ¿Qué momentos discriminan entre señal y ruido, y qué tan sensibles son al nivel de ruido?

The first moment of the signal is: 0.0
The second moment of the signal is: 0.013293333067181288
The third moment of the signal is: 0.006430940666682904
The fourth moment of the signal is: 0.007489819567378666
The first moment of the noise is: 0.0
The second moment of the noise is: 1.0
The third moment of the noise is: 0.005114776505715557
The fourth moment of the noise is: 3.134334832790316

El primer momento central es cero por definición. El segundo momento (la varianza) depende de la normalización, así que es un discriminante débil. El cuarto momento es el interesante: una serie de tiempo que es casi toda ceros con un solo transitorio corto tiene colas pesadas respecto de una gaussiana, así que su curtosis es grande. La curtosis es el caballo de batalla entre las características para detectar eventos impulsivos en el ruido, y volveremos a usarla en la lección de ingeniería de características (2.11).

4. Ruido sintético informado por la física a partir de datos reales

El ruido blanco y el coloreado son idealizaciones. El ruido sísmico real tiene estructura: picos microsísmicos, ruido antrópico durante el día, respuesta instrumental. Para hacer realista el ruido sintético, podemos tomar prestado el espectro de amplitud de ruido real y aleatorizar solo la fase.

Descargamos dos horas de ruido sísmico registrado en la estación UW.RATT antes de un sismo, con el cliente FDSN de obspy.

/Users/marinedenolle/Dropbox/CLASSES/ESS490/curriculum-book/.pixi/envs/default/lib/python3.12/site-packages/obspy/clients/fdsn/client.py:251: ObsPyDeprecationWarning: IRIS is now EarthScope, please consider changing the FDSN client short URL to 'EARTHSCOPE'.
  warnings.warn(msg, ObsPyDeprecationWarning)
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_25811/286277248.py:8: ObsPyDeprecationWarning: attach_response is deprecated and will be removed in a future release. Use remove_response() instead.
  N = fdsn_client.get_waveforms(network=network, station=station, location='--',
1 Trace(s) in Stream:
UW.RATT..HHZ | 2021-07-29T04:15:00.000000Z - 2021-07-29T06:14:59.990000Z | 100.0 Hz, 720000 samples

4.1 Ruido de espectro apareado

Ahora generamos ruido sintético cuyo espectro de amplitud coincide con el del ruido real de RATT. El paquete del curso mlgeo_synth ofrece spectrum_matched_noise, que conserva el espectro de amplitud de un registro de referencia y sortea fases aleatorias. Usa rfft/irfft, así que la simetría hermitiana necesaria para una salida de valores reales queda garantizada por construcción. Ejecute help(mlgeo_synth.spectrum_matched_noise) para ver la firma.

(720000,) (720000,)

Compare ambos en el dominio del tiempo y en el de la frecuencia. Las series de tiempo deberían verse distintas — las fases son diferentes, así que las formas de onda no se alinean en ningún punto. Los espectros de amplitud deberían superponerse exactamente, porque el espectro sintético es una copia del real.

<Figure size 1000x600 with 2 Axes>
<Figure size 1000x500 with 1 Axes>

Los dos espectros caen uno encima del otro: la información de amplitud es idéntica. La información de fase no lo es. Podemos comprobarlo directamente.

<Figure size 1000x400 with 2 Axes>

4.2 Aumentar un evento con ruido realista a SNR variable

Ahora combinamos las dos ideas: un evento sintético (la ondícula de Ricker) más ruido realista de espectro apareado, en un rango de niveles de SNR. Esta es una estrategia estándar de aumento de datos (data augmentation): una sola plantilla de evento limpia produce muchas muestras de entrenamiento con niveles de ruido controlados.

spectrum_matched_noise acepta una longitud de salida n, así que podemos sortear una realización de ruido de la misma longitud que la ventana de la señal.

<Figure size 1000x2000 with 10 Axes>

A SNR baja el evento desaparece dentro del ruido. En algún punto de esta escalera, cualquier detector — humano o algorítmico — empieza a fallar. Encontrar ese punto de falla es exactamente lo que hacemos a continuación, con el detector más antiguo de la caja de herramientas sismológica.

4.3 El generador de datos sintéticos del curso

El paquete mlgeo_synth usado arriba es el generador de datos sintéticos del curso. Además de spectrum_matched_noise, puede generar sismogramas sintéticos completos con una onda P, una onda S, una coda y ruido coloreado, además de conjuntos de datos etiquetados para ejercicios de clasificación. Nos apoyaremos en él en capítulos posteriores, cuando necesitemos datos con verdad de referencia conocida.

<Figure size 1000x400 with 1 Axes>
{'t_p': 10.0, 't_s': 13.571428571428571, 'peak_amplitude': np.float64(0.4841325757612952), 'snr': 5.0}

4.4 Ejemplo desarrollado: el piso de detección del STA/LTA

Antes de que cualquier red neuronal toque trazas como estas (el capítulo 4.3 entrena una), les debemos un modelo de referencia clásico. El detector STA/LTA (promedio de corto plazo sobre promedio de largo plazo; Allen, 1978) desliza dos ventanas sobre la traza: una ventana corta (aquí de 1 s) sigue el nivel instantáneo de la señal, una ventana larga (de 10 s) sigue el fondo, y su cociente se dispara cuando llega un transitorio. Ha corrido en las cadenas de disparo de los observatorios durante décadas, no tiene nada que entrenar, y cualquier detector aprendido debe superarlo para justificar su complejidad.

Un detector necesita un umbral, y el umbral debe salir del ruido solamente. Generamos 200 ventanas de puro ruido, registramos el pico de STA/LTA de cada una y fijamos el umbral en el percentil 99 de esa distribución: una tasa de falsas alarmas del 1 % por construcción, antes de mirar señal alguna. Luego corremos el detector sobre s + noise / SNR a lo largo de la escalera de SNR np.logspace(-1, 2, 20), con 60 realizaciones frescas de ruido por SNR, y contamos la fracción de ensayos en los que el pico de STA/LTA cruza el umbral. Cada punto son 60 ensayos de Bernoulli, así que le adjuntamos un intervalo binomial (de Wilson).

noise-only peak STA/LTA: median 4.34, 99th percentile 6.09
SNR   0.10: detection probability 0.02
SNR   0.14: detection probability 0.02
SNR   0.21: detection probability 0.02
SNR   0.30: detection probability 0.02
SNR   0.43: detection probability 0.05
SNR   0.62: detection probability 0.00
SNR   0.89: detection probability 0.02
SNR   1.27: detection probability 0.03
SNR   1.83: detection probability 0.00
SNR   2.64: detection probability 0.07
SNR   3.79: detection probability 0.22
SNR   5.46: detection probability 0.62
SNR   7.85: detection probability 0.85
SNR  11.29: detection probability 0.98
SNR  16.24: detection probability 1.00
SNR  23.36: detection probability 1.00
SNR  33.60: detection probability 1.00
SNR  48.33: detection probability 1.00
SNR  69.52: detection probability 1.00
SNR 100.00: detection probability 1.00
detection probability crosses 50% near SNR 4.97
<Figure size 1100x400 with 2 Axes>

El panel izquierdo es la lógica del umbral: el ruido puro produce una distribución de picos de STA/LTA (mediana cercana a 4.3 para este ruido rico en microsismos), y colocar el umbral en su percentil 99, 6.1, fija la tasa de falsas alarmas en 1 % antes de que cualquier señal entre al experimento. El panel derecho es el piso de detección: la curva se mantiene en el nivel de falsas alarmas de diseño, de 1-2 %, hasta una SNR de alrededor de 2.6, sube y cruza el 50 % cerca de SNR 5, y satura en 1 por encima de una SNR de alrededor de 16 — una transición de menos de una década de ancho. La moraleja: con su tasa de falsas alarmas clavada en 1 %, este detector clásico necesita una amplitud pico de aproximadamente cinco desviaciones estándar del ruido antes de dispararse con confiabilidad ante este evento. Ese número es el punto de referencia que cualquier detector aprendido tiene que superar; en el capítulo 4.3 superponemos esta curva sobre el piso de detección de una CNN entrenada, medido de la misma manera.

5. Ejercicios

Ejercicio 1: ruido de espectro apareado a partir de otra ventana de referencia

Descargue una ventana de ruido de dos horas distinta de UW.RATT (por ejemplo, 24 horas antes), genere ruido de espectro apareado a partir de ella y compárela con la ventana usada arriba:

  1. Compare la varianza y la curtosis de las dos ventanas reales y de sus dos versiones sintéticas.
  2. Superponga los cuatro espectros de amplitud. ¿Se mueven los picos microsísmicos de un día a otro?
  3. Comente: ¿basta una sola ventana de ruido sintético para representar el ruido de esta estación?

Ejercicio 2: mueva el piso de detección del STA/LTA

El piso de detección de la sección 4.4 se midió con un par de ventanas específico: 1 s la corta, 10 s la larga. Esas longitudes son las únicas perillas de ajuste del detector, y los observatorios las eligen según las señales que persiguen.

  1. Repita el barrido de la sección 4.4 para tres pares de ventanas: (0.5 s, 5 s), (1 s, 10 s) y (2 s, 20 s). Para cada par, recalcule primero el umbral a partir de ventanas de puro ruido — la distribución del ruido cambia con las ventanas — y luego mida la curva de detección y su cruce del 50 %.
  2. Grafique las tres curvas de detección sobre un mismo eje. ¿Qué par de ventanas detecta el evento de Ricker a la SNR más baja, y cómo se relaciona eso con la duración del evento, de aproximadamente 1 s?
  3. ¿Qué sale mal si conserva el umbral de 1 s / 10 s mientras cambia las ventanas?

6. ¿Cuándo son admisibles los datos sintéticos?

Los datos sintéticos son una herramienta, no un sustituto de la observación. Use esta lista de verificación antes de incorporar datos sintéticos a un proyecto.

Admisible:

  • Desarrollo de métodos. Prototipar un detector, un filtro o un modelo sobre datos con verdad de referencia conocida antes de tocar datos reales.
  • Evaluación comparativa (benchmarking). Medir la exactitud (accuracy), la exhaustividad (recall) o los puntos de falla (como en el ejercicio 2) exige conocer la respuesta verdadera. Los datos sintéticos la proporcionan.
  • Aumento de datos con procedencia declarada. Agregar realizaciones de ruido o eventos sintéticos a un conjunto de entrenamiento es práctica estándar, siempre que el artículo o el informe declare qué muestras son sintéticas y cómo se generaron.
  • Conjuntos de prueba ocultos. Los instructores y quienes mantienen los benchmarks usan datos sintéticos para construir conjuntos de prueba que los modelos no pueden haber memorizado.

No admisible:

  • Sustituir la validación real en afirmaciones científicas. Un modelo validado solo con datos sintéticos no ha sido validado. Una conclusión científica sobre la Tierra debe descansar en observaciones reales.

Declare siempre. Cada vez que los datos sintéticos entren en un análisis, dígalo: qué muestras, qué generador, qué semilla. El código del generador y las semillas pertenecen al repositorio junto con el resto del proyecto. En este curso, use el paquete mlgeo_synth como generador para que la procedencia quede a un import de distancia, y siga la guía de proyecto en las instrucciones del proyecto MLGeo.