Explorador de internet GMT Herrmann's Code

Como hacer una tomografía de Ruido Ambiental, paso a paso! 0.- Programas utilizados

Explorador de internet

rdseed

SAC

merge_21_csh (Youkai)

ccNoise (Fenglin) alternativa SAC cor Matlab xcorr

SWT (Fenglin)

GMT Herrmann's Code 1.- Datos 1.1 Búsqueda de datos:

Lo primero que hay que reconocer es qué estaciones existen en el área de estudio. Esto implica que primero hay que verificar qué canales BHZ (banda ancha verticales) han estado activos, en qué tiempo, lo han estado y si el IRIS dispone de los datos. Esto se hace en esta página

http://www.iris.edu/SeismiQuery/channel.htm

Es importante llenar el canal de CHANNEL que es donde se pone BHZ y las coordenadas. Se ve en cuales tiempo las estaciones están todas activas y se seleccionan la ventana de tiempo y las estaciones a utilizar.

1.2 Descarga de datos:

Para descargar los datos se va a esta pagina para hacer la solicitud de las estaciones.

http://www.iris.edu/SeismiQuery/breq_fast.phtml

La solicitud se llena de la siguiente forma

Network: TAL

Station: FULANITA

Channel: BHZ

Data start time: DIAo

Data end time: DIAf

El resto de la info puede ser dejada en blanco

PON UN LABLE!!!!!!!!

Luego

Media: Electronic (ftp)

Query over: DMC archived waveforms

Presionar: START QUERY

Una nueva ventana aparece donde ser revisa la solicitud y se selecciona el dato a descargar, seleccionar full SEED

y presionar SUBMIT

Un email llega al correo con la dirección para bajar los datos.

1.3 Abrir los datos

Para abrir los datos hace falta este programa

http://www.iris.edu/forms/rdseed_request.htm

Con el terminal vamos a la localización del .SEED que hemos descargado

Escribimos: reseed (enter)

Aparece unas preguntas que hay que contestar. Se llena de la siguiente forma:

< IRIS SEED Reader, Release 5.0

Input File (/dev/nrst0) or 'Quit' to Exit: DR.SDD-9-2011-9-2013.536556.seed

Output File (stdout) :

Volume # [(1)-N] :

Options [acCsSpRtde] : d

Summary file (None) :

Station List (ALL) :

Channel List (ALL) : *Z

Network List (ALL) :

Loc Ids (ALL ["--" for spaces]) :

Output Format [(1=SAC), 2=AH, 3=CSS, 4=mini seed, 5=seed, 6=sac ascii, 7=SEGY] : 1

Output file names include endtime? [Y/(N)]Y

Output poles & zeroes ? [Y/(N)]Y

Check Reversal [(0=No), 1=Dip.Azimuth, 2=Gain, 3=Both]: 1

Select Data Type [(E=Everything), D=Data of Undetermined State, M=Merged data, R=Raw waveform Data, Q=QC'd data] :E

Start Time(s) YYYY,DDD,HH:MM:SS.FFFF :

End Time(s) YYYY,DDD,HH:MM:SS.FFFF :

Sample Buffer Length [20000000]:

Extract Responses [Y/(N)] : Y

Se generan un chorreo de archivos .SAC que contienen el dato continuo que se necesita.

2. Preprocesar el dato y calcular la Función de Cross Correlación

El preprocesamiento del dato consiste en: poner todo el dato en un intervalo de muestreo pequeño (5 muestras por segundo o una muestra por segundo), dividir y reorganizar el dato del desastre en el que baja a archivos que contienen un día de información (86400 o 86401 s), remover la respuesta instrumental y finalmente calcular las funciones de crosscorrelación que son las que necesitamos para estimar cosas después.

2.1 Cambiar el intervalo de muestro

Esto se hace con sac. Para saber el intervalo de un archivo, se usa:

r sac

lh delta

para cambiarlo se hace

decimate 4

w over

en el caso de que haya que interpolar porque nos pasamos se usa

interpolate delta 0.2

w over

para muchos archivos se puede usar csh

#!/bin/csh

foreach i (`ls *.SAC`)

echo "…decimating $i..."

sac<<END >/dev/null

r $i

decimate 2

w over

q

END

end

También se puede hacer sin son menos de 1000 archivos

r *.sac

decimate 4

w over

2.2 Dividir en secciones de un día:

Como el dato es continuo, hay que separarlo en pequeños pedazos de un día, para poder manejarlo de una forma más sencilla. Para esto se usa el código de Youkai llamado merge_sac_21. Para usarlo existe una lista que se llama event.lst donde se guardan todos los días del año que se quieren estudiar

Año Mes D Hora Cero Fin del día

2011 01 01 000000.00 86400

2011 01 02 000000.00 86400

2011 01 03 000000.00 86400

2011 01 04 000000.00 86400

2011 01 05 000000.00 86400

Para cordero para una sóla estación hacemos lo siguiente. De aquí saldrán los BHZ de un día de longitud. El código hace lo que puede, rechaza dato muy malo automáticamente y rellena los gaps con la media de los datos

foreach station ( NA.SABA )

if ( ! -e $station ) mkdir $station

cd $station

cp ../event.lst .

foreach comp ( BHZ )

/home/tangyc/bin/merge_sac_21 $station $comp event.lst ../../NewData/$station

end

cd ..

end

2.3 Remover la respuesta del instrumental y remover otras cosas

Una vez que tengamos los datos separados en un día hay que remover una tendencia lineal y la media para cada día (para evitar problemas en el dominio de la frecuencia), usar el filtro taper para hacer cero los extremos y remover la respuesta del instrumento. OJO esta se puede dejar si todos los instrumentos de las estaciones usadas tienen el mismo instrumental. Para removerla hace falta el archivo de polos y ceros que se genera a partir de leer el SEED (e.g. SAC_PZs…). Un ejemplo para hacerlo en una carpeta con muchos BHZ de la misma estación sería

foreach i (`ls *.BHZ`)

echo "…nRESPONSEing $i..."

sac<<END

SETBB pzfile "SAC_PZs"

r $i

rmean

rtrend

taper

bp cor 0.01 0.5 n 2 p 2

transfer from polezero s %pzfile to none

bp cor 0.01 0.5 n 2 p 2

w over

q

END

end

2.4 Calcular la función de Cross-Correlación

Esta función es el primer objetivo a alcanzar y es la que define parcialmente la función de Green (la función de Green describe a dos receptores cómo si ambos fueran o la fuente o el receptor de una onda que pasa entre ellos). Para esto podemos usar el código de Fenglin, alternativamente se puede usar SAC, o Matlab pero no se si hay alguna diferencia entre estos.

Para usar el código de Fenglin se hace

ccNoise -m sta1.SAC -c sta2.SAC -T4096 -N 256 -o output.cc -0 -d -g -p -f0.005/1

T es el tamaño de la ventana N es el número de corridas que hace para el cálculo.

Existe la posibilidad que la amplitud de las señales sea muy pequeña y que el código las desprecie. en ese caso se multiplican para hacerlas grandes y que no haya problema. Por ejemplo

mul 1.e6

un código automático que haga todo para muchas estaciones es:

#!/bin/csh

foreach i (CM.URI CU.ANWB CU.GRGR CU.SDDR DR.SC01 NA.SABA PR.PCDR TR.MCLT TR.SVB TR.TOSP)

foreach j (CM.URI CU.ANWB CU.GRGR CU.SDDR DR.SC01 NA.SABA PR.PCDR TR.MCLT TR.SVB TR.TOSP)

if($i == $j) then

echo "…Misma Estación $i $j…"

else

echo "…Creando directorio…"

mkdir CC.$i.$j

cd CC.$i.$j

set path1 = ../$i

set path2 = ../$j

echo "…Creando path $path1 $path2…"

ls $path1 | awk -F_ '{print $1}' > ls.ls

ls $path2 | awk -F_ '{print $1}' >> ls.ls

sort ls.ls | uniq -d > jday.ls

echo "…Iniciando CC…"

foreach day ( "`cat jday.ls`" )

echo ..CCing jday $day..

ccNoise -m $path1/$day* -c $path2/$day* -T4096 -N 256 -o $i.$j.$day.cc -0 -d -g -p -f0.005/1.

end

cd ../

endif

end

end

Para correrlo en sac es muy sencillo

r sta1.SAC sta2.SAC

cor

w autoSta1.cc CCsta1sta2.cc

Sac calcula automáticamente la autocrosscorrelación para la primera estación así que hay que tener cuidado de no usar este para los estudios. SAC tmb hace unas CC que son del largo de dos días 84600 x2 s, las señales por lo general están entre -800 y +800 s lo cual implica que hay que cortar la CC a ver si de verdad hay señal, mejor no asustarse por el desastre que a veces SAC produce. Matlab tmb las hace con xcorr, pero tienen el mismo detalle.

Para verlas, un filtro bonito es: bp cor 0.03 0.06 n 2 p 2

uno correcto es bp cor 0.02 0.2 n 2 p 2

Para hacer la Pirámide de ruido se utiliza una subrutina de SAC que se llama sis

SAC> read *AZ*BHZ*

SAC> sss

Signal

...