3 SAC - Espectro de Ruido

3.1 Ruido ambiental

Podemos ver el espectro de las vibraciones ambientales que mide un sismómetro para obtener información sobre las fuentes de ruido cerca (o lejos!) de un instrumento y su contribución a la señal.

Dado una señal sísmica, podemos hacer una transformada de Fourier (o, bueno, una transformada rápida de Fourier) para obtener el espectro de ruido en la señal. Vamos a trabajar con datos reales, un registro de un sismómetro que muestra las vibraciones ambientales en la componente vertical: las que vienen de lejos (del mar, por ejemplo), o las que están cerca el sitio (el viento, actividad humana etc.). El archivo de ruido es ruido_TRANB.HHZ, y el archivo de los polos y los ceros en esta instrumento es 3T.pz.

En sismología, siempre el espectro de las señales esta medido en decibelios. Se puede ver las primeras dos hojas de este manual del USGS, por ejemplo, para ver las convenciones usadas: sacpsd.pdf

Abre el programa SAC, y lee el archivo de ruido. Use "lh" y "ppk" para investigar las propiedades de la señal (sugiero "qdp off"): su largo en segundos, sus frecuencias dominantes, y su muestreo (delta).

3.2 Limpieza de señales

Antes de aplicar la transformada rápida de Fourier, debemos limpiarlo. Siempre en la análisis de datos sismológicos uno se elimina el promedio y la tendencia de la serie de tiempo y se aplica un taper (siempre es de Hanning: para que la señal empieza y termina en cero). Los comandos en SAC para hacer eso son "rmean", "rtrend" y "taper". Observe los cambios después de la aplicación de cada comando:

SAC> r ruido_TRANB.HHZ 
SAC> qdp off
SAC> p
SAC> rmean
SAC> p
SAC> rtrend
SAC> p
SAC> taper
SAC> p

¿Por qué es necesario hacer estos cambios?

Si uno quiere, podemos aplicar un filtro también entre ciertas frecuencias, pero vamos a intentar trabajar con datos no-filtrados en esta clase. Finalmente, debemos eliminar la respuesta del instrumento también, usando su archivo de polos y ceros y el comando "transfer" (la última clase). Lo que queremos es la aceleración del suelo en ms$ ^{-2}$, debido a la definición de un espectro de potencia y de los decibelios.

SAC> transfer from polezero subtype 3T.pz to acc freqlimits 
0.001 0.002 500 1000

Note que aquí, fijamos algunos límites de frecuencia al comando de transferencia para evitar problemas con: (i) un sismómetro tiene cero respuesta a cero frecuencia; y (ii) a altas frecuencias el muestreo puede causar problemas a veces. Esta modificación de la señal, en el dominio de frecuencia, es unidad (ningún cambio) entre los valores centrales, y cero en los extremos. En el caso de arriba, hemos elegido un rango de frecuencias mas amplio que el rango de operación del instrumento. Se puede leer mas acerca de este comando en el manual.

3.3 El espectro de potencia

Para una serie de tiempo discreto, tomamos el espectro de frecuencia usando la transformada rápida de Fourier con el comando "fft". Por el espectro de potencia estamos interesados en la información de la amplitud, no la fase, entonces usamos "keepam".

SAC> fft
SAC> keepam

Si aplicamos la transformada de Fourier a la señal medida en ms$ ^{-2}$ (con la respuesta del instrumento ya eliminado), estamos listos para convertir el espectro a un espectro de potencia en unidades de decibelios. Este requiere las siguientes operaciones:

  1. Hay que elevar el espectro al cuadrado para obtener potencia.
  2. Multiplicamos el espectro por un factor de 2 para tomar en cuenta la contribución de frecuencias negativas (convención USGS).
  3. Hay que dividir el espectro por el largo de la serie de tiempo en segundos, para normalizarlo (nuestro archivo es 1000) segundos en largo.
  4. Hay que tomar el logaritmo (base 10) del espectro, y multiplicarlo por 10, para obtener el resultado en decibelios [ $ 10\log_{10}$m$ ^2/$s$ ^4/Hz$] [dB].

En SAC, la serie de comandos necesarios es:

SAC> sqr
SAC> mul 2
SAC> div 1000
SAC> log10
SAC> mul 10

El espectro de potencia es creado en este punto. Vamos a verlo; es una buena idea definir los limites del eje-x en el gráfico, y que sea una escala logarítmica, con "xlim" y "xlog":

SAC> xlim 0.01 60
SAC> xlog
SAC> p

¿Cuál es su reacción inicial en ver el espectro? Cabe mencionar en este punto que el espectro va a tener muchas fluctuaciones sistemáticas en su amplitud debido al hecho que la señal original no es continua ni infinita. Podemos suavizar el espectro con una media móvil: para cada punto en el espectro, su valor es el promedio del punto mas los 80 puntos alrededor de él. El comando en SAC ... "smooth":

SAC> smooth mean halfwidth 40

Ahora esta listo el espectro. Debe aparecer como lo en la figura 3.1. Sugiero guardar el espectro con

w espectro_ruido_TRANB.sac

¿Por qué hay un gran cambio en el espectro alrededor de 50 Hz?

En la sección 3.1, notaron las frecuencias dominantes en la señal. ¿El espectro que ustedes tomaron esta de acuerdo con su respuesta en la sección 3.1?

3.4 Los modelos de ruido alto y ruido bajo del USGS

De un estudio del USGS sobre su red mundial, definieron límites para ruido alto y ruido bajo en estaciones sísmicas. El informe se puede bajar del siguiente enlace, los gráficos en las páginas 14, 31 y 35 son los mas relevantes para esta clase.

Enlace: peterson_usgs_seismic_noise_ofr93-322.pdf

Se puede bajar los dos archivos que definen espectros de ruido "alto" y "bajo" en estaciones sísmicas aquí:

New Low Noise Model: nlnm_freq.acc
New High Noise Model: nhnm_freq.acc

Los archivos son de texto que se puede leer con el comando "readtable". Los siguientes comandos en SAC generan el imagen en la figura 3.1:

SAC> readtable content p nlnm_freq.acc
SAC> w nlnm_freq.sac 
SAC> readtable content p nhnm_freq.acc
SAC> w nhnm_freq.sac 
SAC> r espectro_ruido_TRANB.sac nlnm_freq.sac nhnm_freq.sac
SAC> p2

\includegraphics[width=12 cm]{espectro_ruido_TRANB.eps}
Fig 3.1: El espectro del ruido del archivo ruido_TRANB.HHZ, comparado con los modelos de ruido alto y bajo del USGS.

Use su entendimiento de la clase teórica para explicar la forma del espectro de ruido en esta estación.

3.5 Espectro de potencia - función de probabilidad

Para una estación sísmica, se pueden calcular miles de espectros de ruido dentro de un periodo de tiempo de varias días y amontonarles. La figura 3.2 muestra la probabilidad que el espectro de ruido calculado por una hora de datos cae dentro de un cierto píxel en el gráfico. Con eso, se puede ver el promedio del ruido y una estimación de la variación en el ruido con el tiempo.

Para mas información se puede leer el informe del USGS:

McNamaraBuland_usgs_psdpdf_3.0.pdf


\includegraphics[width=12 cm]{1.ps}
Fig 3.2: Espectro de ruido para la estación TRANB (la función de densidad de probabilidad para 30 días de datos continuos).

  1. ¿Qué significan los valores de la probabilidad en este gráfico?
  2. Sismos locales, microsismos del mar, ondas de superficie. ¿Cuáles influencian la variación en el espectro a un periodo de (i) 6 s; (ii) 20 s; (iii) 0.5 s.

3.6 Trabajo adicional

El archivo 2004361005853_XJ_TRANB.HHZ es la señal registrada en la estación TRANB (Puerto Río Tranquilo, Aysén) de un famoso terremoto histórico.

Este archivo tiene un largo de 7000 segundos y muestra, mayormente, una onda Rayleigh. El espectro de potencia de un terremoto va a estar muy distinto del espectro de ruido ambiental. En SAC, es posible aplicar un filtro (tipo pasa banda) usando el comando "bandpass":
SAC> bandpass co f1 f2 np 6 p 1
Este es un filtro con frecuencias de esquina de f1, f2 [Hz], con 6 polos y 1 paso (si quieren saber mas de los polos y pasos de este filtro pregunta en las clases, o alternativamente se puede investigar su efecto sobre una función impulso, en el dominio de frecuencia, en el estilo de la última parte de la clase 2).

About this document ...

This document was generated using the LaTeX2HTML translator Version 2008 (1.71)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 -white -no_navigation sac_03_ruido.tex

The translation was initiated by matt on 2016-12-21


matt 2016-12-21