Story Transcript
UNIVERSIDAD COMPLUTENSE DE MADRID FACULTAD DE CIENCIAS FÍSICAS DEPARTAMENTO DE FÍSICA DE LA TIERRA, ASTRONOMÍA Y ASTROFÍSICA I (GEOFÍSICA Y METEOROLOGÍA)
PROGRAMA OFICIAL DE POSTGRADO DE FÍSICA: MÁSTER DE GEOFÍSICA Y METEOROLOGÍA TRABAJO DE INVESTIGACIÓN
LOCALIZACIÓN DE TERREMOTOS A PARTIR DE FORMAS DE ONDA Almudena Gomis Moreno Directora: Dra. Elisa Buforn Peiró Junio 2007
AGRADECIMIENTOS En primer lugar quisiera agradecer al Departamento de Física de la Tierra, Astronomía y Astrofísica I de la Facultad de Físicas de la Universidad Complutense de Madrid el haberme facilitado los medios para poder llevar a cabo este trabajo de investigación del Máster de Geofísica y Meteorología. En particular, querría mostrar mi más profundo agradecimiento a la directora de este trabajo, la Dra. Elisa Buforn Peiró, sin cuya ayuda y consejos no se podría haber realizado. Este trabajo ha sido llevado a cabo gracias al proyecto ERSE (Escenarios Realistas de Riesgo Sísmico en España) (REN2003-5178-C03-01), financiado por el Ministerio de Ciencia y Tecnología (en la actualidad Ministerio de Educación y Ciencia) y al proyecto RISTE (Riesgo de Terremotos y Tsunamis en España) (CGL2006-10311-C03-01/BTE), financiado por el Ministerio de Educación y Ciencia. Gracias a Hélène Lyon-Caen, Francesco Pacchiani y Aurore Franco, de la Ecole Normal Supérieur de París, por dedicar su tiempo a enseñarme los conceptos fundamentales y los algoritmos, y facilitarme los programas para llevar a cabo este trabajo. Quisiera agradecer también al Instituto Geográfico Nacional su colaboración con la aportación de los registros e información instrumental necesarios para el desarrollo de este proyecto. Mis compañeros de estos dos años han sido una parte importantísima de este proceso. Gracias por todos los consejos, los ánimos, los fantásticos momentos que hacen que todo sea más fácil, y por todas las cosas que he aprendido de vosotros: Tatiana, Carmen, Juan Luis, Paco, Dani, Jacques, Esther, Ana, Javi. Todos habéis puesto vuestro granito de arena. Tatiana, muchísimas gracias por el apoyo en todos los sentidos, incluso logístico, hasta el último minuto.
i
Gracias a Salva y Lucía, que, en estos dos años y hasta el último momento, me han sacado de más de un apuro. Sin todas esas pequeñas ayudas del día a día no se podría acabar ningún trabajo de este tipo. De manera especial quiero expresarle mi agradecimiento a Joaquín por sus ratos de escucha, consejos e interés, cuando la geofísica no le toca ni de cerca, y por sus intentos de atraerme al lado oscuro. Tentaciones de ese tipo siempre son bienvenidas. ¡Qué bonita es la Física!, ¿verdad? Carmen, gracias, porque has hecho que sienta que esta etapa de mi vida no está aislada de todo lo demás, sino que sea continuidad de todo lo que pasamos juntas en el colegio. Llevo dos semanas intentando encontrar el agradecimiento adecuado para Juan, y todo resulta insuficiente. Solo sé que me siento muy afortunada de que seas tú el que tenga que aguantar mis enfados (bendita paciencia), y con el que pueda compartir todas mis alegrías. En último lugar quiero darles las gracias a mis padres, del los que he aprendido a querer aprender. Me habéis dado todo vuestro apoyo en los momentos más difíciles y me habéis ayudado a elegir el camino correcto en los momentos claves.
Y puede que este sea un buen principio. Principio de incertidumbre…*
*
“Principio de incertidumbre”. Ismael Serrano.
ii
LOCALIZACIÓN DE TERREMOTOS A PARTIR DE FORMAS DE ONDA
Trabajo de Investigación elaborado por Almudena Gomis Moreno, encuadrado en el Programa Oficial de Postgrado de Física: Máster en Geofísica y Meteorología, impartido por los Departamentos de Física de la Tierra, Astronomía y Astrofísica I y II de la Facultad de Ciencias Físicas de la Universidad Complutense de Madrid.
Madrid, 15 de junio de2007
Almudena Gomis Moreno
iii
ÍNDICE
1. Introducción ....................................................................................... 1 2. Metodología ........................................................................................ 3 2.1. Identificación de familias o multipletes....................................... 6 2.2. Correlación cruzada ................................................................... 8 2.3. Relocalización de los eventos.................................................... 10 3. Implementación del algoritmo........................................................... 12 3.1. Determinación de las familias de sismos .................................. 12 3.1.1. Cálculo de la transformada wavelet de las ondas P .......... 12 3.1.2. Normalización de la wavelet............................................. 13 3.1.3. Cálculo de las distancias de Levenstein ........................... 14 3.1.4. Creación de un diagrama de árbol ................................... 14 3.1.5. Extracción de los multipletes de sismos .......................... 15 3.2. Cálculo de la correlación cruzada............................................. 16 3.2.1. mkhdr ............................................................................. 16 3.2.2. mkdat ............................................................................. 16 3.2.3. mklist.............................................................................. 16 3.2.4. mktbl .............................................................................. 17 3.2.5. dbh ................................................................................. 17 3.3. Relocalización .......................................................................... 17 4. Aplicación a la serie sísmica de 2002 de Gérgal (Almería).................. 18 4.1. Introducción ............................................................................ 18 4.2. Procesamiento de datos............................................................ 20 4.3. Resultados ............................................................................... 23 5. Conclusiones .................................................................................... 34 6. Bibliografía ....................................................................................... 36 7. Anexo 1: sismograma y sonograma ................................................... 39
iv
ÍNDICE DE APÉNDICES
Se adjunta un CD-ROM en el que se tiene acceso a los siguientes apéndices: A. Programa en .bash para la automatización de los programas param, scale, dist-spares y cluster-sparse B. Programa para el cálculo de coordenadas UTM C. Tiempos de llegada manuales de la fase P D. Resultados de la relocalización con Hypo71 E. Diagramas de árbol de la distancia de Levenstein F. Extracción de multipletes G. Archivos con extensión .hdr H. Archivos con extensión .dat I. Otros I.1. Lista de sismos de Gérgal :gergal.lis I.2. Lista de sismos de la familia A: mult0001.lis I.3. Lista de sismos de la familia B: mult0002.lis I.4. Información de la familia A: mult0001.table I.5. Información de la familia B: mult0002.table I.6. Retardos calculados para la familia A I.7. Retardos calculados para la familia B
v
1. INTRODUCCIÓN La ocurrencia de terremotos es uno de los fenómenos naturales que mayores daños económicos y humanos puede causar. El estudio detallado de la ocurrencia de sismos es fundamental para poder prevenir y mitigar los efectos destructores de los mismos. Uno de los problemas fundamentales en el estudio de terremotos es su localización espacio–temporal, que permitirá identificar las fuentes sismogénicas de la región y asociarlas a las fallas activas. Los primeros métodos y más sencillos para la localización de un sismo (determinación hipocentral) fueron los métodos gráficos, que requerían un mínimo de tres estaciones para obtener la localización. Posteriormente se desarrollaron los métodos numéricos. Hoy en día, la mayoría de los métodos numéricos que se utilizan se basan en la minimización de los residuos entre los tiempos teóricos y los observados de llegadas de las diversas fases del sismograma: las ondas P, S, Lg, etc. Esta minimización realizada para diferentes modelos de Tierra permite obtener las coordenadas del epicentro, la profundidad y la hora origen del terremoto. Este es el fundamento de los métodos más utilizados, como el Hypo71 (Lee y Lahr, 1972), HYPOINVERSE (Klein, 1978), etc. De todos los parámetros a determinar, la profundidad es el que, con mayor frecuencia, queda peor determinado, con grandes errores, al no converger el sistema de ecuaciones y debiendo fijarse este parámetro. Este problema se presenta normalmente cuando existe una mala geometría de la red y se ve muy afectado por las discontinuidades en los modelos de Tierra, estructurados en capas plano-paralelas de velocidad, lo que le resta exactitud a la relocalización hipocentral. Recientemente se han desarrollado nuevas metodologías basadas en la utilización de formas de onda, que permiten determinar con mayor precisión la profundidad. Los métodos que utilizan formas de onda se basan en la comparación de la forma de onda de distintos sismos en cada una de las estaciones en las que se han registrado, agrupándolos en familias según la similitud de sus formas de onda (Poupinet et al., 1984; Ito, 1990). Una vez identificadas las diferentes familias, se procede a computar el retardo entre las diferentes señales, usando métodos de correlación cruzada (Waldhauser y Ellsworth, 2002; Hauksson y Shearer, 2005; Shearer et al., 2005). La hipótesis en que se basa el método consiste en que la diferencia de tiempos de llegada de una misma fase de dos sismos, con un mecanismo focal semejante y próximos en el espacio, es debida a la diferencia hipocentral de los
1
mismos. Conocido el retardo será posible calcular la posición relativa de los hipocentros de estos terremotos mediante distintos procesos, como por ejemplo la doble diferencia de Waldhauser y Ellsworth (2000). La aplicación de estos métodos permite una relocalización más exacta que la obtenida con los métodos tradicionales, minimizando la influencia del modelo de Tierra en la relocalización. Hay que tener en cuenta que, debido a que las zonas sísmicamente más activas de España son el sur de la Península Ibérica, el mar de Alborán, el Golfo de Cádiz y el norte de Marruecos, la cobertura azimutal de estaciones no es óptima, lo que produce en algunos casos grandes errores en la profundidad. Una metodología como la que aquí se plantea permitirá mejorar la relocalización hipocentral, lo que implica una mejor definición de las zonas sismogénicas y su posterior aplicación a estudios de sismotectónica o de riesgo sísmico. El objetivo de este trabajo es la relocalización por medio de formas de onda de terremotos ocurridos en la región Ibero-Magrebí (sur de la Península Ibérica, Golfo de Cádiz, mar de Alborán y norte de África). En concreto se centrará en la denominada serie de Gérgal (Almería), que comenzó el 4 de febrero de 2002 con la ocurrencia de un terremoto principal de magnitud mbLg = 5.1 que fue seguido por un total de 39 sismos que se extendieron hasta el día 26 de dicho mes. Para llevar a cabo este proceso se van a emplear registros digitales de banda ancha y de período corto. La serie de Gérgal ha sido determinada previamente por el Instituto Geográfico Nacional (IGN) utilizando métodos tradicionales de lectura de fases P y S. Los resultados obtenidos en este trabajo utilizan formas de onda, y serán comparados con los resultados obtenidos con la metodología tradicional utilizada por el IGN. El trabajo se va a desarrollar según el siguiente esquema: 1. identificación de posibles familias, 2. determinación hipocentral de parámetros focales por medio de formas de onda, 3. comparación de los resultados con los obtenidos por el IGN.
2
2. METODOLOGÍA La localización de un terremoto consiste en determinar las coordenadas epicentrales, profundidad y hora origen del sismo (xo, yo, zo, to), minimizando en la medida de lo posible sus errores. Para un grupo de terremotos, la localización puede ser: a) absoluta, en la que se localiza de manera absoluta el hipocentro de los sismos, o b) relativa, en la que los sismos se relocalizan en relación a otro sismo o al baricentro de un grupo de sismos. La precisión de la localización absoluta depende de muchos factores, tales como la geometría de la red, la precisión de la toma del tiempo de llegada de la fase y el conocimiento del medio (Pavlis, 1986). El método que se va a emplear en este trabajo se basa en la determinación del tiempo de retardo de una misma fase de dos sismos diferentes registrados en la misma estación, lo que se llama doble diferencia (Waldhauser y Ellsworth, 2000). Un esquema de esta situación se muestra en la figura 1. Puede decirse que la doble diferencia proporciona una medida de la distancia entre los dos hipocentros (Menke, 1999). El resultado es una relocalización relativa de los sismos. La ventaja de utilizar relocalizaciones relativas es que minimizan los errores introducidos por el medio (Poupinet et al., 1984; Frechet, 1985; Got et al., 1994).
Figura 1. Disposición de los hipocentros y la estación. DD indica doble diferencia.
3
Se define el retardo entre dos señales, o doble diferencia, como la diferencia del tiempo de recorrido desde el foco del sismo a una misma estación para una fase determinada en dos sismos diferentes. Este retardo se debe únicamente a la diferencia en las localizaciones hipocentrales de los sismos que las han generado, si se cumplen las siguientes condiciones: 1. que sus mecanismos focales sean similares, 2. que sus focos estén próximos, en comparación con la distancia foco-estación, de modo que las ondas hayan seguido aproximadamente la misma trayectoria y hayan atravesado el mismo medio. Un criterio a seguir es que dos sismos tendrán sismogramas incoherentes cuando sus hipocentros estén separados más de la cuarta parte de la longitud de la frecuencia predominante en el sismograma (Geller y Mueller, 1980). El conocimiento del medio siempre es necesario, lo que constituye una dificultad añadida a los procesos de relocalización, ya que el medio atravesado por las ondas rara vez es bien conocido. Pero en el caso de dos sismos con hipocentros muy próximos la influencia del medio es mucho menor, ya que la influencia de las anomalías de velocidad es igual para ambos y, por tanto, no tienen influencia directa en la relocalización. No obstante, en ausencia de errores proporcionados por el medio, en la relocalización siempre estaría presente el error de la determinación manual de fases mediante un programa de procesado sísmico, que en los mejores casos es de tan solo 0.1s (Lyon-Caen et al., 2004), lo que conlleva un error epicentral de 500 m y de alrededor de 1 Km en profundidad. Por tanto, la gran ventaja del empleo de la correlación cruzada para el cálculo de la doble diferencia es que las incertidumbres de las coordenadas hipocentrales disminuyen en más de un orden de magnitud en comparación con otras técnicas que no emplean la doble diferencia (Waldhauser y Ellsworth, 2000), ya que se elimina la incertidumbre que acarrea determinar de manera manual el tiempo de llegada de las distintas fases que van a intervenir en la relocalización. La precisión temporal de este método puede llegar a ser de 1 ms frente al rango de 10-30 ms en el caso de la determinación manual del tiempo de llegada, lo que proporciona una relocalización relativa entre terremotos con un error que oscila entre unos metros a unas decenas de metros. El error del que no nos vamos a poder desprender es el introducido por la geometría de la red. Una mayor cobertura azimutal proporciona unos mejores resultados.
4
En el proceso de relocalización, el cálculo del retardo se puede llevar a cabo básicamente de dos maneras distintas: 1. tomando el tiempo de llegada de la fase directamente sobre el sismograma y calculando la diferencia, 2. empleando el método de correlación cruzada entre dos señales que son similares. La ventaja de la correlación cruzada es que proporciona menor error que la lectura de fases de un sismograma. La correlación cruzada determina directamente el retardo entre dos sismos, lo que conlleva asociado un error más pequeño. Algunos autores (Pacchiani, 2006) denominan relocalización relativa a la relocalización llevada a cabo con la correlación cruzada y relocalización absoluta a la que emplea el tiempo de llegada de cada fase, si bien no hace referencia al proceso de relocalización en sí. Es decir, a si los sismos son relocalizados en función de otro sismo o de un baricentro, o si lo son de manera absoluta. El procedimiento a seguir en este trabajo es el que se muestra a continuación: 1. Identificación de familias o multipletes. Un multiplete, o familia, se define como un conjunto de terremotos cuyos sismogramas tienen un alto grado de similitud debido a que han tenido lugar muy próximos en el espacio y con un mecanismo focal muy parecido. Gracias a esta característica, se puede calcular la doble diferencia para una misma fase mediante la correlación cruzada. A dos sismos pertenecientes a una misma familia se les denomina doblete. 2. Correlación cruzada. Es la comparación de la forma de onda de los terremotos de un mismo multiplete. Se procede al cálculo del retardo de la misma fase en cada par de sismogramas. Esta diferencia temporal es fundamental para la relocalización hipocentral. Es la etapa más importante del proceso de relocalización, pues de su precisión dependerá el éxito de la relocalización. 3. Relocalización de los eventos. Mediante cualquier método que emplee el retardo en la llegada de una fase para cada par de estaciones. Si se relocalizan los sismos en función de un sismo principal el método se denomina “master event” (Ito, 1985; Deichmann y García-Fernández, 1992). Si por el contrario la relocalización se lleva a cabo empleando todos y cada uno de los sismos de la familia como evento maestro, el método se
denomina “master event
generalizado”. Asimismo se puede emplear el programa de relocalización HypoDD (Waldhauser y Ellsworth, 2000), en el que se puede llevar a cabo tanto
5
una relocalización a partir de retardos de correlación cruzada, como con retardos obtenidos a partir de medidas absolutas de los tiempos de llegada, como con una combinación de ambos tipos. La solución se halla mediante mínimos cuadrados, de manera iterativa. A continuación se describe en detalle cada una de estas etapas del proceso de relocalización. 2.1. Identificación de familias o multipletes Para la agrupación de los sismos en distintas familias es necesaria la obtención previa de sus sonogramas. Se define un sonograma como la representación del contenido en frecuencias de un sismograma frente al tiempo. Esta representación se obtiene al aplicar transformadas wavelet a los sismogramas correspondientes a los eventos sísmicos de interés. La transformada wavelet consiste en ir analizando la señal temporal (sismograma) usando unas determinadas funciones (wavelet madre) de distintos tamaños que se trasladan en el tiempo a lo largo de la señal original. Así se conocerá el contenido de las distintas frecuencias para cada instante de tiempo. La definición matemática de la transformada wavelet es la que sigue:
Ws (τ , s ) =
1 s
∫
+∞
−∞
⎛ t −τ ⎞ x(t )ψ * ⎜ ⎟dt ⎝ s ⎠
[1]
donde τ es el parámetro de traslación, s es el parámetro de escala, x(t) es la señal procesada y ψ(t) la wavelet madre. Hay muchos tipos diferentes de wavelets madre. En este trabajo se va a usar la wavelet de Morlet o gaussiana (Goupillaud, Grossman y Morlet,1984), definida como:
ψ (t ) = e
⎞ ⎛ −t 2 ⎟ 2⎟ ⎜⎜ t Δ o ⎠ ⎝
e i 2πf ot
[2]
que se encuentra centrada en to en el dominio del tiempo y en fo en el de frecuencias, con una duración Δt=s Δto en tiempo y Δf= Δfo/s en frecuencias.
6
A partir de la wavelet madre van a ir surgiendo otras funciones con distintos valores de τ y de s para llevar a cabo la transformada. El parámetro de traslación, τ, está relacionado con la posición de la función wavelet sobre la señal que queremos transformar: se refiere a la localización de la ventana de análisis. Corresponde a la información en el tiempo que va a proporcionar la transformada. Por su parte el parámetro de escala, s, está relacionado con el contenido en frecuencias: el escalamiento de una señal corresponde a una compresión o a una dilatación de la wavelet madre, por lo que, si la wavelet madre es estrecha, se estará obteniendo información de las altas frecuencias, y por el contrario, si es ancha, se obtiene información de las bajas frecuencias. En la figura 2 se muestra un ejemplo de una señal y de una wavelet madre para valores de s de 1 y 5 , y valores de to de 2, 40, 90 y 140 s, donde la zona sombreada en naranja es la wavelet que se desplaza a lo largo de la señal.
Figura 2. Aplicación de una wavelet (zona sombreada) sobre una señal, para diferentes valores de s y to. Las unidades del eje de ordenadas son cuentas, y en abscisas segundos. (Modificado de Polikar, 1994). La transformada wavelet se aplicada a cada uno de los sismogramas con los que vamos a trabajar. El resultado de la transformación wavelet es un sonograma. El siguiente paso es determinar qué sismos forman parte de un mismo multiplete, y para ello es necesario estudiar la similitud que hay entre los sonogramas que se han obtenido al aplicar la transformada wavelet a cada uno de los sismogramas. La similitud viene determinada por la distancia de Levenstein, que se define como el número de cambios, supresiones o inserciones que se han de llevar a cabo para transformar un patrón (en este caso los sonogramas) en otro (Levenstein, 1965). Para comprender mejor este proceso supongamos el siguiente ejemplo (figura 3): dos patrones formados por cuatro cuadrados de distinto color cada uno. Para
7
transformar uno en el otro bastará con intercambiar la posición del cuadro blanco y la del naranja. En este caso la distancia de Levenstein es 1.
Figura 3. Ejemplo del cálculo de la distancia de Levenstein. Por tanto, el siguiente paso a seguir es calcular la distancia de Levenstein para cada par de sonogramas de una misma estación. Pero hay que tener en cuenta lo siguiente: cada sismo ha sido registrado en varias estaciones, por lo que no es inmediato que dos sonogramas que sean similares en una estación impliquen que los sismos correspondientes sean similares. Para concluir que dos terremotos forman parte de un mismo multiplete será necesario un alto grado de similitud de los sonogramas en el mayor número de estaciones posible. Este número de estaciones ha de ser determinado dependiendo del problema que se esté tratando. 2.2. Correlación cruzada Una vez se han obtenido los multipletes o familias de sismos se procede a calcular el retardo de las parejas de señales pertenecientes a un multiplete mediante una correlación en el dominio de las frecuencias (Poupinet et al, 1984; Got et al., 1994), si bien también se puede llevar a cabo en el dominio temporal. La similitud entre dos señales viene determinada por el índice de correlación, que varía entre –1 y 1. Un valor de 1 indica la máxima similitud, 0 la mínima. Un valor negativo significa que la señal está invertida. Para comprender el proceso partamos de dos señales dependientes del tiempo, s1(t) y s2(t), entre las que se cumple la siguiente relación:
s2 (t ) = C te s1 (t + τ )
8
[3]
Para ver el grado de correlación entre dos señales, se comparan en el dominio de las frecuencias. Por tanto, el siguiente paso es expresar las señales en el dominio de la frecuencia, para lo que se aplica la transformada de Fourier: +∞
S1 ( f ) = ∫ s1 (t )e −2πift dt −∞
+∞
S 2 ( f ) = ∫ s 2 (t )e −∞
− 2πift
dt = C
te
∫
+∞
−∞
s1 (t )e
− 2πif (t +τ )
dt = C S1 ( f ) e te
− 2πifτ
[4]
El grado de correlación entre dos señales en el dominio de la frecuencia viene dado por el índice de coherencia:
γ(f )=
S12
2
[5]
A1 ⋅ A2
donde S12 es el interespectro, y A1 y A2 los autoespectros, que se definen como
S12 ( f ) = S1 ( f ) ⋅ S 2* ( f ) = C te S1 ( f ) e 2πifτ
[6]
A1 = S1 ⋅ S1*
[7]
A2 = S 2 ⋅ S 2*
[8]
2
Para el cálculo del retardo entre dos señales que son suficientemente similares se emplea la fase del interespectro, que se obtiene de la expresión [6]:
θ12 ( f ) = 2πτf = τω
[9]
El retardo entre las dos señales es la pendiente de la recta de la fase del interespectro frente a la frecuencia
ω . La pendiente está ponderada mediante unos
pesos w(f) que dependen del índice de coherencia, dando así mayor peso allí donde existe mayor semejanza entre las dos señales:
γ 2( f ) w( f ) = (1 − γ 2 ( f ))
9
[10]
Por tanto, la regresión en cada estación está limitada a un número determinado de frecuencias alrededor del máximo del interespectro. Para que una frecuencia se pueda considerar apta para tomar parte en la regresión, es necesario que cumpla dos condiciones: 1. que la energía asociada a esa frecuencia sea por lo menos el 75 % de la energía máxima del sismograma, 2. que la coherencia sea al menos de 0.8. 2.3. Relocalización de los eventos Una vez que se han obtenido los retardos para cada par de sismos ya se puede proceder a la relocalización. Para ello se puede usar cualquier método que emplee la doble diferencia. En general, el tiempo de recorrido de una onda a una estación se puede expresar en función de las coordenadas del sismo (xo, yo, zo), las de la estación i-ésima (xi, yi, zi), la hora origen (to) y la velocidad, v, como
t obs = t o +
( x i − x o )2 + ( y i − y o )2 + ( z i − z o )2 v
= to +
Δi = F (xo , y o , z o , t o ) v
[11]
En esta ecuación son cuatro las incógnitas que se han de determinar (xo, yo, zo, to). Se supone conocida la distribución de velocidades y que la lectura del tiempo de llegada de las fases está libre de error, pero esto no es cierto, por lo que es preferible tener un sistema sobredeterminado. Incluso si esto no pudiera ser, se puede llegar a encontrar la mejor solución posible por medio de una inversión. Lo más habitual es tratar de minimizar iterativamente el residuo mediante mínimos cuadrados. Se define el residuo como
Ri = t i
obs
− ti
cal
[12]
donde i hace referencia al sismo. El tiempo calculado de llegada de una determinada fase se puede expresar como
t cal = F ( xi , y i , z i , t i )
10
[13]
Como el problema directo no es lineal, para resolver el inverso se procede a desarrollar F(xi, yi, zi, ti) en serie de Taylor (Geiger, 1910) en torno a unos parámetros iniciales (xo+dxo, yo+dyo, zo+dzo, to+dto). Para facilitar la notación llamemos
r m =(m1, m2,
m3, m4)=(x0, y0, z0, t0); entonces el tiempo calculado, si nos quedamos en primer orden, se puede expresar como
(x − x ) (y − y ) (z − z ) r Rki = t obs − Fi m 0 = i o o Δx0 + i o o Δy 0 + i o o Δz 0 − Δt 0 vΔ i vΔ i vΔ i
( )
[14]
donde el subíndice k hace referencia a la estación. La expresión anterior se puede abreviar como
Rki =
∂t ki r r Δmi ∂m
[15]
Esta ecuación no es apropiada para emplear la doble diferencia, pero se puede definir una diferencia de residuos para dos sismos, i y j, en la misma estación k-ésima (Frechet, 1985):
(
dR = t − t ij k
i k
)
j obs k
(
− t −t i k
)
j cal k
∂t ki r ∂t kj r = r Δmi − r Δmi ∂m ∂m
[16]
En esta ocasión se pueden emplear tanto las dobles diferencias, como los tiempos de llegada absolutos obtenidos de la observación de sismogramas. Este
problema
se
resuelve
mediante
inversión,
bien
mediante
la
descomposición en valores singulares (SVD) para los grupos de sismos pequeños, bien por el método de los gradientes conjugados (LSQR), (Paige y Saunders, 1982) para sistemas con más ecuaciones.
11
3. IMPLEMENTACIÓN DEL ALGORITMO La metodología descrita en el apartado anterior se ha aplicado mediante diversos programas implementados por varios autores. Estos programas han sido adaptados para poder usarlos con los datos utilizados en este estudio, siendo necesario modificar la razón de muestreo de nuestras señales (Δt), el formato de los nombres de nuestros archivos de entrada, el rango de frecuencias de interés en nuestro trabajo, y la estructura de carpetas, es decir, la disposición de los directorios y archivos en el ordenador. A continuación se describen brevemente las principales características de los algoritmos utilizados de acuerdo con la metodología discutida en el apartado anterior. 3.1. Determinación de las familias de sismos. Se ha utilizado la metodología SPARS (Syntactic Pattern Recognition Scheme), (Zhizhin, 1992), de reconocimiento sintáctico de patrones. Para ello se ha implementado un programa en .bash (Apéndice A del CD-ROM que se adjunta) para la ejecución automática de los cuatro primeros programas que obtienen las familias de sismos. A continuación se describen brevemente los algoritmos utilizados y los parámetros de entrada, así como los resultados que proporciona cada uno de ellos. 3.1.1. Cálculo de la transformada wavelet de las ondas P. Se ha realizado mediante el algoritmo param (programa lenguaje en C; Zhizhin, 1992). Los parámetros de entrada de este programa son los siguientes: 1. n, el número de intervalos de frecuencia en los que se va a estudiar el contenido frecuencial (oscila entre 1 y 20). Se han tomado 20 intervalos para hacer un estudio más exhaustivo de los sismogramas, 2. u, límite superior del intervalo de frecuencias, expresado en Hz. No debe ser mayor que la frecuencia de Nyquist, que se define como fN=1/2Δt, donde Δt es el intervalo de muestreo, que en este caso es Δt = 0.02, y por tanto fN = 25 Hz. Se ha tomado u=10 Hz, ya que se va a proceder a filtrar los sismogramas entre 1 y 10 Hz, debido a que se ha observado que de este modo se obtiene una correlación cruzada óptima (Haukkson y Shearer, 2005),
12
3. l, frecuencia inferior del intervalo de frecuencias, también en Hz. No puede ser menor que la frecuencia de resolución, que se define como la inversa de la duración de la señal. Tomamos l=1 Hz, debido al filtrado paso-banda al que sometimos a los sismogramas, tal y como se ha explicado en el párrafo anterior, 4. d, duración de las ventanas temporales. Debe ser mayor que el ritmo de muestreo, pero menor que la duración total del intervalo. Tomamos un valor de 0.032 s, que es el más comúnmente empleado, ya que en nuestros sismogramas Δt=0.02 s. El resultado final de este programa es la transformada wavelet de la señal cortada. En nuestro caso, la fase utilizada es la onda P. Se definen los parámetros IOFFSET (señal antes de la llegada de la onda P) y WFLENGTH (longitud total de la señal). Los valores utilizados son 20 y 256 puntos, respectivamente, lo que asegura que la ventana utilizada contenga la onda P. 3.1.2. Normalización de la wavelet. El algoritmo empleado en esta etapa es scale (programa en lenguaje C; Zhizhin, 1992). Con él se normaliza la transformada wavelet y se convierten las unidades a dB. El resultado es el sonograma de la señal sísmica correspondiente. En la figura 4 se muestra un ejemplo del sonograma para la componente z del terremoto del 5 de febrero de 2002 en la estación ACLR. En el sismograma, el eje de ordenadas son cuentas, y en el sonograma, el intervalo de frecuencia empleado. En ambos casos el
Cuentas digitales
eje de abscisas son segundos.
t (s)
Figura 4. En la parte superior se muestra la ventana temporal utilizada para la obtención del sonograma –parte inferior-del terremoto del 5 de febrero de 2002 en la estación ACLR. Se puede consultar esta figura en mayores dimensiones en el Anexo 1 al final de este trabajo.
13
3.1.3. Cálculo de las distancias de Levenstein. Esta etapa se llevó a cabo mediante el algoritmo distpar (programa en lenguaje C; Zhizhin, 1992 modificado por Gaucher, 1996), cuya entrada son los archivos KnowBaseName y TestBaseName, que contienen la lista de eventos que se van a procesar, y LocBaseName, que contiene la lista de eventos con su localización preliminar. Asimismo, se implementó un programa en MATLAB (Apéndice B del CD-ROM) que transforma coordenadas geográficas en coordenadas UTM empleando el elipsoide de referencia del WGS84 (sistema de coordenadas mundiales que data de 1984), necesarias para el cálculo de la distancia euclidiana, que es la que es calculada por el programa. El resultado de este programa son las distancias de Levenstein entre todos los sismogramas, comparándolos de 2 en 2. 3.1.4. Creación de un diagrama de árbol. Este paso se lleva a cabo mediante el programa cluster-sparse (programa en lenguaje C; Zhizhin, 1992). Ordena las distancias de Levenstein obtenidas por el método anterior en un diagrama de árbol que muestra la similitud entre los sismos. Los sismos relacionados con las ramas más cortas del árbol son más semejantes que los relacionados mediante las ramas largas, entendiéndose que pueden pertenecer a una misma familia. En la figura 5 se muestra un ejemplo de diagrama de árbol para la estación EBAN. Las distancias de Levenstein aparecen a la derecha, oscilando entre 3.396×105 y 1.577×106.
14
Tree[root].Dist = 1577359.000000 2002.02.04-20.13.32--------------------------------> 8.682e+05+--------> 2002.02.04-21.37.53-----------------------> | | 6.240e+05+--------< | 2002.02.04-20.44.31-------------------> | | 5.010e+05+---< | 2002.02.05-05.28.20-------------------< | 1.106e+06+----------------> 2002.02.05-00.19.41-----------------------> | | 6.242e+05+---> | | 2002.02.14-03.16.10-----------------------< | | | 7.185e+05+> | | 2002.02.04-20.31.02-----------------------> || | | 6.184e+05+--- | | 5.963e+05+-----< | 2002.02.04-20.17.44----------------------> | 5.740e+05+ | 2002.02.05-13.51.22----------------------< | 1.577e+06+ 2002.02.05-12.41.16---------------------------------> | 8.840e+05+------------------------< 2002.02.04-21.43.21-----------------------> | 6.009e+05+----> | 2002.02.04-22.11.25-----------------------< | | 7.586e+05+----< 2002.02.04-22.00.03-------------> | 3.396e+05+--------------< 2002.02.05-21.42.00-------------<
Figura 5. Ejemplo de un diagrama de árbol para la estación EBAN. Los números a la izquierda identifican al sismo. Los restantes corresponden a las distancias de Levenstein entre unos sismos y otros. Los círculos rojos indican la mayor y la menor distancia. 3.1.5. Extracción de los multipletes de sismos. Esta fase se ha realizado mediante el programa mkclus (programa en lenguaje C; adaptado de Pacchiani, 2006). Obtiene las familias de sismos mediante la comparación de los diagramas de árbol de las distintas estaciones obtenidos en el apartado anterior. Determina qué sismos se parecen y en qué estaciones, a fin de, imponiendo un número mínimo de estaciones en las que los sismogramas deben parecerse, extraer los multipletes de manera objetiva. A continuación se describen los parámetros que controlan el proceso, así como sus valores óptimos: 1. número de estaciones en que se busca coincidencia: dos seísmos similares deben tener registros parecidos en el máximo número de estaciones posible. Restringir mucho este parámetro provoca familias poco pobladas, pero se asegura la similitud de los sismos que la constituyen. El inconveniente es que
15
se pueden discriminar sismos que pertenecen a esa misma familia, pero de los que no se tiene un registro suficientemente bueno. Por el contrario, un valor pequeño produce familias muy numerosas formadas por sismos que no tienen por qué ser similares realmente. Por tanto, es necesario llegar a un compromiso. En este trabajo se ha estudiado el rango comprendido entre 1 y 5 estaciones, 2.
distancia de Levenstein: indica el grado de similitud que se requiere para considerar que en una estación dos sismos son coincidentes. A menor umbral, más parecidos serán los sismos. Se ha trabajado con valores que oscilaban entre un valor de 4×105 y 10×105. Las familias propuestas por mkclus han sido confirmadas visualizando los
sismogramas correspondientes. 3.2. Cálculo de la correlación cruzada. Como paso previo a la correlación cruzada es necesario obtener una primera localización. Se ha utilizado el programa Hypo71, que proporciona un fichero de salida con extensión .vis, del que los programas utilizados en la etapa de correlación leerán la información de cada sismo, como es el nombre de las estaciones en las que se tiene registro, el azimut o el ángulo de incidencia. Para llevar a cabo esta relocalización es necesario tomar en los sismogramas la hora de llegada de las fases implicadas en el proceso de relocalización. En nuestro caso se trabajará con ondas P. Los algoritmos empleados para llevar cabo la correlación cruzada son: 3.2.1. mkhdr (programa en lenguaje Fortran; autor desconocido; modificado de Pacchiani, 2006). Escribe las cabeceras de los sismos, que contienen las coordenadas del hipocentro, hora origen, azimut, ángulo de incidencia y tiempo de llegada de la P, 3.2.2. mkdat (programa en lenguaje Fortran; autor desconocido; adaptado de Pacchiani, 2006). Corta los sismogramas con una ventana de una longitud determinada. En este trabajo se ha usado la zona de la onda P, empleando una ventana de 512 puntos de longitud, 100 puntos por delante de la llegada de la onda P. Se modificó el programa para poder visualizar los sismogramas cortados con cualquier programa de hoja de cálculo o similar, 3.2.3. mklist (programa en lenguaje Fortran; modificado de Pacchiani, 2006). Genera una lista con los sismos pertenecientes a un mismo multiplete,
16
3.2.4. mktbl (programa en lenguaje Fortran; Got, 1999; modificado por Pacchiani, 2006. Adaptado de Pacchiani). Crea dos archivos con información de las estaciones y de los rayos: distancia, azimut, ángulo de incidencia y tiempo de viaje de la onda P. Estos archivos son necesarios para el cálculo del retardo, 3.2.5. dbh (programa en lenguaje Fortran; modificado de Got, 1999; basado en el algoritmo de Frechet, 1985). Calcula el retardo de los sismos dos a dos, indicando la coherencia entre ellos. Fue modificado para que el formato de salida coincidiera con el de entrada del programa HypoDD, usado para la relocalización. 3.3. Relocalización. Este paso se llevó a cabo mediante el programa HypoDD (Waldhauser y Ellsworth, 2000), al que se puede acceder de manera libre en la red. Se trata de un algoritmo que permite realizar la relocalización relativa de sismos, tanto a partir de medidas de correlación cruzada, como a partir de tiempos absolutos del tiempo de llegada de la onda P, o mediante un procedimiento mixto. Tan sólo hubo que compilarlo y ponerlo al día a partir de la fe de errores proporcionada por los autores. Al ser un algoritmo suficientemente conocido, no se describe en detalle en esta memoria.
17
4. APLICACIÓN A LA SERIE DE 2002 DE GÉRGAL (ALMERÍA). El método descrito se va a aplicar a la serie de Gérgal con el objetivo de obtener una localización más exacta de los hipocentros, mejorando las profundidades con respecto a la solución calculada por el Instituto Geográfico Nacional (IGN) 4.1. Introducción El terremoto principal de la serie de Gérgal tuvo lugar el 4 de febrero de 2002 a las 20:09:30.48 (localización del IGN) (figura 6) con magnitud mbLg=5.1, y fue sentido con intensidad V en las poblaciones de Gérgal, Ohanes y Tabernas. Fue precedido por un terremoto precursor el mismo día a las 19:55:31.82 de magnitud mbLg=2.3 y seguido por un total de 39 réplicas (tabla 1), cuyos epicentros (IGN) se muestran en la figura 7. En esta localización se observa que los epicentros se distribuyen siguiendo una línea en dirección N-S a lo largo de una extensión de unos 160 Km de largo, con una anchura de unos 50 Km aproximadamente. Las profundidades determinadas por el IGN son superficiales (menores de 40 Km), si bien de 14 sismos no se ha podido calcular la profundidad.
Figura 6. Sismicidad de la Península Ibérica y norte de África. Se representan sismos con M≥3.0, en el período 1980-2004. En rojo los sismos superficiales (h ≤ 40 Km), en verde los de profundidad intermedia (40 Km ≤ h ≤ 150 Km), y en azul los profundos (600 Km ≤ h ≤ 700 Km). La estrella negra determina la localización del sismo principal de la serie de Gérgal (modificada de Goded, 2006)
18
Figura 7. Epicentros de los 39 sismos considerados de la serie de Gérgal según la localización inicial del IGN. El círculo azul corresponde al sismo principal. El tamaño de los puntos indica la magnitud del terremoto correspondiente. El triángulo verde indica la posición de la estación ACLR, que es la más cercana a la zona de estudio.
TABLA 1 SISMOS DE LA SERIE DE GÉRGAL (IGN) FECHA 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002 04/02/2002
HORA 19:55:31.82 20:09:30.48 20:13:33.20 20:15:16.51 20:17:45.12 20:18:23.92 20:20:09.25 20:24:59.86 20:31:02.64 20:34:08.35 20:35:24.40 20:40:14.88 20:44:32.54 21:09:35.20 21:37:54.77 21:43:21.71 22:00:04.88 22:11:26.89
LATITUD 37,0972º 37,0931º 37,0489º 37,0575º 37,1093º 37,0358º 37,0994º 37,0911º 37,0916º 37,0922º 37,0739º 37,0144º 37,0463º 37,0439º 36,9773º 37,0932º 37,1016º 37,0778º
19
LONGITUD -2,5347º -2,5379º -2,5001º -2,4993º -2,5366º -2,5318º -2,5251º -2,5211º -2,5852º -2,5355º -2,5230º -2,5282º -2,5326º -2,5201º -2,5444º -2,5462º -2,5457º -2,5470º
H (KM) 0,0f 0,8 11,3 0,0f 4,9 0,5 10,8 0,0f 1,3 0,0f 5,9 0,0f 2,7 0,0f 6 1,1 0 0,0f
mbLg 2,3 5,1 2,7 2,2 2,4 2,3 2,1 2,5 2,6 3,0 2,5 1,8 2,4 2,0 2,3 2,9 2,8 2,8
TABLA 1 CONTINUACIÓN FECHA 04/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 05/02/2002 06/02/2002 06/02/2002 10/02/2002 10/02/2002 14/02/2002 15/02/2002 26/02/2002
HORA 23:56:30.18 00:03:57.83 00:19:42.46 02:21:02.69 05:28:20.50 06:08:10.96 08:40:12.57 12:41:17.88 12:51:07.12 13:51:23.39 14:21:22.77 18:06:03.43 20:00:25.99 21:42:01.22 07:06:18.40 21:21:48.08 02:42:39.67 19:28:02.48 03:16:11.43 13:17:48.00 17:19:44.43
LATITUD 37,0955º 37,0970º 37,0749º 37,0194º 37,0622º 37,0386º 37,1054º 37,0622º 37,0703º 37,0395º 37,0494º 37,0155º 37,0485º 37,0723º 37,0587º 37,0554º 37,1118º 37,0940º 37,1079º 37,0967º 37,0924º
LONGITUD -2,5447º -2,5407º -2,5021º -2,5431º -2,5417º -2,5336º -2,5471º -2,5556º -2,5446º -2,5268º -2,5331º -2,5331º -2,5351º -2,5582º -2,5458º -2,5573º -2,5403º -2,5525º -2,5558º -2,5441º -2,5448º
H (KM) 2,7 1,4 10,6 11,4 10,6 16,1 6,4 10,5 3,8 4,6 0,0f 7,8 0,0f 1,2 11,2 0,0f 0,0f 7,3 0,0f 0,0f 0,0f
MBLG 1,6 2,2 2,2 1,7 2,4 2,2 2,4 3,1 2,3 2,2 2,1 2,1 1,9 3,0 2,6 2,0 1,7 2,1 1,8 2,0 2,2
4.2. Procesamiento de datos Para la relocalización de esta serie por medio de formas de onda, se han utilizado los datos correspondientes a las estaciones digitales de banda ancha y de período corto del Instituto Geográfico Nacional (IGN) y de la red conjunta del Real Observatorio de la Armada de San Fernando y de la Universidad Complutense de Madrid (ROA-UCM). En las tablas 2 y 3 se muestran las coordenadas y características de las estaciones empleadas en este trabajo. Su localización se muestra en la figura 8. TABLA 2 ESTACIONES IGN Estación EADA EALH EBAN EBER ECOG EALB EMAL EGUA ELOJ ELUQ ENIJ ERON EQES ETOB EVIA
Latitud 38.1673º 37.8582º 38.1710º 36.8978º 37.2772º 35.9398º 36.7620º 36.8337º 37.1480º 37.5592º 36.9715º 37.0180º 37.8028º 38.6447º 38.6387º
Longitud N N N N N N N N N N N N N N N
4.5772º 1.4197º 3.7900º 2.8897º 3.5663º 3.0343º 4.4280º 3.5653º 4.1530º 4.2678º 2.2070º 3.8050º 3.0712º 1.5478º 2.5025º
W W W W W W W W W W W W W W W
20
Muestreo (muestras/s) 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
Tipo banda ancha período corto período corto banda ancha período corto banda ancha banda ancha período corto período corto período corto período corto período corto banda ancha banda ancha período corto
TABLA 3 ESTACIONES ROA-UCM
-8º
Estación
Latitud
Longitud
ACBG ACLR ALB CART CEU CFS ELUQ ORGV
36.7672º 37.1885º 35.9300º 37.5868º 35.9160º 35.2300º 37.5605º 36.8605º
N N N N N N N N
-6º
2.1950º 2.5832º 3.0300º 1.0120º 5.3300º 2.4870º 4.2668º 3.4293º
Muestreo (muestras/s) 50 50 50 100 / 80 100 / 80 50 50 50
W W W W W W W W
-4º
EADA
0º ETOB
38º EALH
358û
EQES
CART
ELU y LUQ ACLR
ECOG ELOJ ERON EMAL
Gergal EBER
ENIJ
ORGV EGUA
ACBG
MA
0û
36º
CEU
-6º
-4º
R
M
ITE ED
R
N RA
EO
36º
ALB y EALB
MELI
-8º
ancha ancha ancha ancha ancha ancha ancha ancha
EBAN
PENINSULA IBERICA
LIJA
Banda Banda Banda Banda Banda Banda Banda Banda
-2º EVIA
38º
Tipo
AFRICA CFS
-2º
0º
Figura 8. Distribución de las estaciones sísmicas de la red IGN y ROA-UCM empleadas en este trabajo. El punto rojo indica el epicentro del sismo principal de la serie de Gérgal (IGN). De acuerdo con la metodología descrita, la primera etapa es la identificación de la fase P en los sismogramas, para lo que se han utilizado los registros en formato SAC (Seismic Analysis Code; Goldstein, 1987) (Apéndice C del CD-ROM). En algunos de los sismogramas se podía intuir la onda S, sin embargo, no se ha tenido en cuenta en este trabajo por dos motivos fundamentales:
al no observar la fase S claramente, el error introducido al tenerla en cuenta podía ser más grande que las ventajas derivadas de considerar una nueva fase,
para poder llevar a cabo la correlación cruzada es necesario observar la misma fase en bastantes sismogramas y, a ser posible, en las mismas estaciones para poder llevar a cabo una comparación satisfactoria. Con el número de fases S que se podía llegar a identificar, no se cumplía este requisito.
21
Para la identificación de la fase P se han seleccionado las estaciones óptimas considerando como tales las estaciones más cercanas, hasta un máximo de 214 Km, que son las que proporcionan un buen registro de la señal, con una relación señal/ruido baja. Especial interés tienen las estaciones ACLR, de la red ROA-UCM, situada en Calar Alto (a 10 Km del epicentro del sismo principal) y las estaciones EBER y ENIJ del IGN, que se encuentran a menos de 40 Km (figura 8). A partir del tiempo de llegada de las ondas P se ha realizado una primera localización con el programa Hypo71 (Apéndice D del CD-ROM). El siguiente paso consistió en remuestrear las señales sísmicas con el mismo Δt, ya que la razón de muestreo oscila entre 50 y 100 muestras/s según la estación (tablas 2 y 3). Es un requisito necesario para comparar unos sismogramas con otros, tanto en el paso para obtener las familias de sismos, como para calcular la doble diferencia. El muestreo elegido en este trabajo ha sido 50 muestras/s, es decir Δt=0.02s. Tanto la identificación de la onda P como el remuestreo se han realizado con SAC. Es importante comprobar que los espectros de los sismogramas remuestreados no varían respecto a los sismogramas originales, para asegurarse de que el contenido en frecuencias sigue siendo el mismo después del remuestreo. En las figuras 9a y 9b se puede ver un ejemplo, donde se puede apreciar cómo la resolución ha disminuido ligeramente al reducir el ritmo de muestreo, observándose menos detalles en la nueva señal sísmica. Sin embargo, los espectros de la señal original y remuestreada son idénticos, permaneciendo inalterados los valores de amplitud en su parte plana (2×101 cuentas) y la frecuencia esquina (7 Hz).
Figura 9 a. Sismograma original y remuestreado a 50 muestras/s del sismo del 5 de febrero de 2002 a las 05:28:20.50. Los círculos rojos indican dos instantes de tiempo para los cuales se aprecia, a modo de ejemplo, una significativa pérdida de resolución.
22
AMPLITUD
AMPLITUD
FRECUENCIA (Hz)
FRECUENCIA (Hz)
Figura 9 b. Espectros de amplitud del sismograma original y del remuestreado. En esta fase del tratamiento de datos es necesario llevar a cabo un filtrado de los registros, pues muchos de ellos llevan un ruido aleatorio asociado que puede interferir a la hora de proceder a la correlación cruzada. Asimismo, la correlación cruzada es más fiable si se lleva a cabo un filtrado entre 1 y 10 Hz (Hauksson y Shearer, 2005). De modo que se procedió a aplicar un filtro Butterworth paso-banda con 3 polos, entre 1 y 10 Hz. Con este paso queda finalizado el tratamiento de datos y se puede comenzar con el proceso de extracción de familias. 4.3. Resultados El siguiente paso consiste en obtener los sonogramas y normalizarlos, para finalmente obtener las familias de sismos. Empleando un valor umbral de 6×105 en la distancia de Levenstein se han obtenido dos familias, a partir de los 39 sismos iniciales de partida. La primera de ellas, la familia A, consta de 7 sismos, y la familia B de 6 (tabla 4). Con los sismos restantes no se ha podido llevar a cabo la relocalización, ya que han sido descartados debido a que, con la distancia de Levenstein elegida como umbral, no formaban parte de ninguna familia. Se pueden consultar los diagramas de árbol obtenidos en el Apéndice E del CD-ROM, y el estudio de la influencia del número de estaciones y de la distancia de Levenstein1 en la extracción de familias en el Apéndice F.
ID sismo 1 2 3 4 5 6 7
1
04 04 04 04 04 05 10
TABLA 4 FAMILIAS DE SISMOS Familia A Familia Fecha Hora Fecha / 02 /2002 20:13:32 04 / 02 /2002 / 02 /2002 20:18:22 04 / 02 /2002 / 02 /2002 20:34:07 05 / 02 /2002 / 02 /2002 21:43:21 05 / 02 /2002 / 02 /2002 22:00:03 05 / 02 /2002 / 02 /2002 21:42:00 05 / 02 /2002 / 02 /2002 02:42:39
Apartado 3.1.5., p. 15.
23
B Hora 20:44:31 21:37:53 05:28:20 06:08:10 08:40:11 12:51:06
En la figura 10 se muestran dos ejemplos de sismos de una misma familia, uno para cada caso, y en dos estaciones distintas, ACLR y EBER. Se puede apreciar que el contenido en frecuencias es prácticamente igual. La diferencia reside en la amplitud de las fases, debido a que se trata de sismos de diferente magnitud. Para la familia A se han seleccionado el sismo del 4 de febrero a las 20:18, que tuvo una magnitud mbLg=2.3, y el del 10 de febrero a las 02:42, con mbLg=1.7. La diferencia en magnitud es lo suficientemente grande como para apreciar diferencias en amplitud en los sismogramas. En la familia B se han seleccionado sismos del 5 de febrero, a las 05:28 y 06:08, con magnitud mbLg de 2.4 y 2.2 respectivamente. En este caso la diferencia de magnitud es menor, por lo que las amplitudes son más similares.
A
B
Figura 10. Ejemplos de dos sismos pertenecientes a la familia A (en la parte superior) y B (en la parte inferior). Una vez obtenidas las diferentes familias de terremotos, se procede al cálculo de la doble diferencia para todos los pares de sismos de una misma familia (resultados intermedios en los Apéndices G, H e I del CD-ROM). Los resultados obtenidos de la doble diferencia se muestran en las tablas 5 y 6, en las que se indica la identificación de los sismos que se van a correlacionar dentro del multiplete (ID1 e ID2), y el número de estaciones en las que se ha llevado a cabo el cálculo del retardo. A continuación se indica la estación, el retardo (×10-4 s) para esa estación y el peso, definido según la ecuación [10] que se le asigna a ese valor. Por ejemplo, en la primera línea de la tabla, ID1=1 e ID2=2 indican que se van a correlacionar los sismos 1 y 2 en seis estaciones diferentes. En el segundo caso, ID1=1 e ID2=3, señalan una correlación entre los sismos 1 y 3 en 10 estaciones distintas. Y así sucesivamente.
24
TABLA 5 RETARDOS DE LA FAMILIA A ID1 Sta 1 ACBG EGUA 1 ACBG EBER ENIJ EVIA 1 ACBG EBER ELUQ EVIA 1 ACBG EBER ENIJ EVIA 1 ACBG EBER ENIJ 1 ACLR EGUA 2 ACBG EGUA 2 ACBG EGUA 2 ACBG EGUA 2 ACBG EGUA 2 ACLR ENIJ 3 ACBG EBAN ELUQ EVIA 3 ACBG EBAN ENIJ EVIA 3 ACBG EBAN ENIJ
ID2 retardo ×10-4 s 2 1945.7 2112.6 3 657.2 2203.1 2003.6 4381.9 4 2177.3 3967.3 1980.9 975.0 5 5134.3 5482.7 5189.5 1581.8 6 -368.6 5988.3 4591.6 7 2478.3 2427.7 3 -1588.6 4198.5 4 184.0 3728.3 5 3201.1 3585.1 6 1822.7 2071.8 7 86.9 -83.6 4 1316.1 -3522.3 -714.3 -358.8 5 2438.0 2204.3 2232.7 -3343.4 6 2898.4 3853.4 2069.9
nº de estaciones en Peso Sta w 6 0.8291 ACLR 0.8689 ENIJ 10 0.7723 ACLR 0.7833 EGUA 0.5865 ORGV 0.647 10 0.6850 ACLR 0.8887 ECOG 0.6485 ENIJ 0.5335 10 0.8390 ACLR 0.9357 ECOG 0.9280 ORGV 0.7438 9 0.6382 ACLR 0.8693 ECOG 0.7980 ORGV 5 0.9105 EBER 0.9238 ENIJ 6 0.8618 ACLR 0.7796 ENIJ 6 0.8280 ACLR 0.8216 ENIJ 6 0.9206 ACLR 0.7770 ENIJ 6 0.8385 ACLR 0.7418 ENIJ 4 0.9687 EBER 0.9492 10 0.7965 ACLR 0.8584 EBER 0.6672 ENIJ 0.8081 10 0.8437 ACLR 0.9136 EBER 0.7821 ORGV 0.7257 8 0.9278 ACLR 0.7316 EBER 0.8156 ORGV
que se lleva a cabo la correlación cruzada. retardo Peso Sta retardo ×10-4 s ×10-4 s
Peso
2389.1 4325.1
0.9678 0.8847
EBER ORGV
1989.4 429.3
0.9599 0.7085
3110.4 2478.7 -3288.7
0.7195 0.8235 0.8138
EBAN ELUQ EQES
3341.9 -463.7 -289.2
0.8138 0.6420 0.7045
0 9344.8 4528.5
0 0.6602 0.6831
EBAN EGUA ORGV
9036.9 2511.9 -2356.9
0.8312 0.7970 0.8623
5496.1 7328.8 -1316.8
0.9038 0.7503 0.8189
EBAN EGUA EQES
5314.8 5694.0 3524.6
0.9389 0.9001 0.7068
3621.3 -291.6 349.2
0.8164 0.6461 0.7025
EBAN EGUA ERON
9703.8 412.2 8366.7
0.6126 0.7969 0.7988
255.7 2154.8
0.8800 0.9335
ECOG
835.4
0.6620
816.6 808.1
0.8239 0.6498
EBER ORGV
231.2 131.7
0.8520 0.7487
-729.6 709.9
0.8270 0.8318
EBER ORGV
1974.0 869.0
0.9328 0.7303
3146.2 2955.7
0.9517 0.9853
EBER ORGV
3521.8 3145.7
0.9808 0.8362
1216.6 2798.2
0.8183 0.8997
EBER ORGV
3970.7 3730.4
0.9364 0.8030
155.3
0.8883
EGUA
1863.9
0.7654
513.2 1914.7 2744.5
0.9017 0.8834 0.8665
EALH EGUA ORGV
-951.0 1354.2 1557.3
0.9720 0.9595 0.8107
0 3322.0 3046.8
0 0.8853 0.9252
EALH EGUA EQES
3205.2 3256.2 4139.9
0.9590 0.9451 0.8153
3638.8 3649.7 -152.4
0.9382 0.9638 0.7778
EALH EGUA
3701.3 3597.7
0.9859 0.9464
25
TABLA 5 (CONTINUACIÓN) 3 ACBG ENIJ 4 ACBG EALH CART ENIJ 4 ACBG EBAN ECOG ORGV 4 ACLR EGUA 5 ACBG EBAN ECOG ORGV 5 ACLR EGUA 6 ACLR EGUA
7 2356.5 -818.9 5 879.3 4380.4 985.4 958.9 6 3411.6 -199.1 76.5 360.4 7 -864.4 -194.0 6 485.3 2180.0 612.1 657.3 7 -3064.0 0 7 -1435.6 -557.5
4 0.8319 0.7735 12 0.9010 0.9899 0.7356 0.8624 10 0.8074 0.9053 0.8445 0.7865 5 0.7865 0.8244 10 0.8816 0.8308 0.9825 0.9432 5 0.9804 0 5 0.7699 0.7496
EBER
-10.2
0.7255
EGUA
-23.6
0.8832
ACLR EBAN ECOG ORGV
543.4 -4451.4 1498.0 46.0
0.7351 0.9330 0.8856 0.8449
EADA EBER EGUA EVIA
2897.6 1519.2 1909.3 539.3
0.7352 0.9615 0.9450 0.6681
ACLR EBER EGUA
-321.9 1928.1 0
0.7660 0.9688 0
EALH CART ENIJ
4622.8 4404.5 2124.7
0.9696 0.7306 0.8543
EBER ENIJ
-1826.7 -793.8
0.8356 0.8093
ECOG
-1662.9
0.6702
ACLR EBER EGUA
-1817.3 411.3 411.6
0.9048 0.9955 0.9639
EALH CART ENIJ
444.5 -1701.8 -176.7
0.9289 0.8780 0.9457
EBER ENIJ
-3365.5 -3045.0
0.8944 0.9557
ECOG
-3137.8
0.7649
EBER ENIJ
-6928.6 -2861.2
0.8126 0.9075
ECOG
-3771.5
0.7257
TABLA 6 RETARDOS DE LA FAMILIA B ID1 Sta 1 ACBG EBER ENIJ 1 ACBG EBER ENIJ 1 ACBG EBER ORGV 1 ACBG EGUA 1 ACLR EGUA 2 ACBG EBER ENIJ 2 ACBG EBER ORGV
ID2 retardo ×10-4 s 2 -3841.3 -2154.8 1179.9 3 7044.1 -1263.6 1610.7 4 -325.3 -1123.7 -1113.6 5 -1361.4 -327.0 6 -955.4 -281.3 3 -1000.8 886.6 298.8 4 907.2 1046.1 831.1
nº de estaciones en que se lleva Peso Sta retardo ×10-4 s 9 0.8192 ACLR -2927.4 0.9879 ECOG -1603.3 0.9079 ORGV 1015.5 8 0.7334 ACLR -4065.7 0.9813 ECOG 78.7 0.9173 ORGV -1351.7 7 0.6993 ACLR -3725.7 0.9568 EGUA -1208.1 0.8492 6 0.8050 ACLR 1124.9 0.7436 ENIJ 5017.9 6 0.8514 EBAN 619.7 0.8128 ENIJ 1628.0 8 0.7980 ACLR -1109.0 0.9854 ECOG -33.4 0.9694 ORGV 662.8 7 0.7694 ACLR -778.2 0.9660 EGUA 775.6 0.8489
26
a cabo la correlación cruzada. Peso Sta retardo ×10-4 s
Peso
0.8560 0.7955 0.7991
EBAN EGUA EQES
-165.3 -2352.0 1333.7
0.7934 0.7838 0.6909
0.8653 0.8859 0.8197
EBAN EGUA
-1038.3 -1511.7
0.7939 0.9146
0.8742 0.8773
EBAN ENIJ
0 1829.9
0 0.9153
0.8134 0.9231
EBER ORGV
1863.6 1899.5
0.9338 0.7911
0.7104 0.6862
EBER ORGV
1931.5 1903.8
0.9864 0.6903
0.9626 0.8521 0.7733
EBAN EGUA
-961.9 565.1
0.8985 0.8957
0.9712 0.8446
EBAN ENIJ
-664.1 568.0
0.8399 0.9856
TABLA 6 (CONTINUACIÓN) 2 ACBG EGUA 2 ACLR EGUA 3 ACBG EBER ORGV 3 ACBG EGUA 3 ACLR EGUA ERON 4 ACBG EGUA 4 ACLR EGUA EVIA 5 ACLR ENIJ
5 3875.8 3708.4 6 1934.1 1817.1 4 2026.4 160.9 1884.1 5 3056.9 3169.5 6 2928.7 0 9387.2 5 1255.7 3053.7 6 2604.2 2844.5 -1978.8 6 0 681.8
6 0.8729 0.9179 6 0.9884 0.7906 8 0.8771 0.9695 0.7802 6 0.7356 0.9637 8 0.9598 0 0.6020 6 0.6615 0.8441 7 0.9711 0.7854 0.7198 5 0 0.7405
ACLR ENIJ
4081.1 3799.7
0.8236 0.9643
EBER ORGV
4026.7 3868.4
0.9317 0.8631
EBAN ENIJ
5066.8 482.6
0.7274 0.7461
EBER ORGV
3945.9 -1116.9
0.9624 0.7839
ACLR EGUA EVIA
336.2 159.7 -4123.0
0.9964 0.8908 0.6388
EBAN ENIJ
268.8 270.3
0.9268 0.9729
ACLR ENIJ
3852.8 3457.6
0.9211 0.9652
EBER ORGV
3160.6 -372.1
0.9506 0.6873
EBAN ENIJ EVIA
2993.2 4024.5 -1194.1
0.7105 0.7062 0.6053
EBER ORGV
3116.0 -541.2
0.9713 0.7814
ACLR ENIJ
3500.5 3209.0
0.8994 0.9688
EBER ORGV
2979.3 3039.9
0.9472 0.8624
EBAN ENIJ
2608.8 -80.7
0.7320 0.7392
EBER ORGV
2974.1 -601.9
0.9553 0.7202
EBER ORGV
-13.0 -1846.6
0.9779 0.8392
EGUA
-34.3
0.7563
Un valor nulo del retardo en una estación indica que en esa estación los registros han sido lo suficientemente distintos como para no poder comparar la señal, y por tanto se les asocia un peso nulo. Se puede apreciar que el peso no depende únicamente de la distancia entre la estación y el hipocentro. Si así fuera, las estaciones más alejadas tendrían asociado un peso menor. Pero, tal y como se explicó en la metodología, el peso depende del índice de coherencia entre las dos señales que se correlacionan, dando mayor peso a señales más coherentes entre sí. Por tanto se observa que, en algunos casos, a retardos calculados en estaciones más lejanas se les ha asociado un peso mayor. No obstante, en términos generales, sí se observa que en las estaciones más cercanas se obtienen retardos con un mayor peso asociado, como es el caso de ACLR, ENIJ o EBER, y retardos con pesos menores en las más lejanas, particularmente en EVIA. La explicación reside en que, para estaciones más cercanas las señales que se registran generalmente son mejores, con una mejor relación señal/ruido, lo que conlleva que la correlación cruzada sea mejor, y por tanto, el peso también mayor. La excepción es la estación EALH, situada a unos 130 Km del epicentro del baricentro de la familia A, que proporciona unos pesos entre 0.9 y 1.0 (frente a pesos
27
entre 0.6 y 0.7 en estaciones más cercanas), si bien la correlación solo ha sido posible para 5 pares de sismos, por lo que no se considera significativo este resultado. La estación que proporciona en general los peores retardos es ELUQ, con pesos muy bajos, que varían entre 0.6 – 0.7. La estación que proporciona el peso menor de todos es EVIA, que es la estación más alejada del epicentro del sismo principal. Con estos resultados se procede a llevar a cabo la relocalización mediante el programa HypoDD para cada una de las familias por separado. El modelo de Tierra empleado fue obtenido por Dañobeitia et al. (1998) (figura 11). Es un modelo específico para el área de estudio, frente al modelo empleado por el IGN, de carácter más general. Se ha elegido este modelo de Tierra a partir de los resultados obtenidos por Rodríguez Abad (2004), que concluye que este modelo es más adecuado para la relocalización del sismo principal de esta serie.
Velocidad (Km/s) 5.70 5.75 5.90 5.95 6.70 6.75 7.90
Profundidad (Km) 0.00 8.00 12.00 18.00 23.00 29.00 35.00
Figura 11. Modelo de Tierra empleado en este trabajo. (Dañobeitia, 1998). Los resultados de la relocalización se muestran en las tablas 7, 8, 9 y 10, comparados con la localización inicial dada por el IGN.
TABLA 7 RELOCALIZACIÓN FINAL FAMILIA A. ID 1 2 3 4 5 6 7
Fecha 04 04 04 04 04 05 10
/ / / / / / /
02 02 02 02 02 02 02
/2002 /2002 /2002 /2002 /2002 /2002 /2002
Hora origen (hh:mm:ss.ss) 20:13:32.88 20:18:23.00 20:34:07.84 21:43:21.00 22:00:03.60 21:42:00.08 02:42:39.18
RMS
Latitud
Longitud
0.70 0.24 0.47 0.56 0.49 0.50 0.23
37.0903 37.0944 37.0951 37.0941 37.0944 37.0929 37.1009
-2.5656 -2.5580 -2.5622 -2.5661 -2.5611 -2.5619 -2.5600
28
Profundidad (Km) 9.1 ± 1.4 9.0 ± 4.5 10.0 ± 1.3 9.3 ± 1.3 8.8 ± 1.3 10.0 ± 1.3 10.4 ± 2.2
TABLA 8 LOCALIZACIÓN IGN FAMILIA A. ID 1 2 3 4 5 6 7
Fecha 04 04 04 04 04 05 10
/ / / / / / /
02 02 02 02 02 02 02
/2002 /2002 /2002 /2002 /2002 /2002 /2002
Hora origen (hh:mm:ss.ss) 20:13:33.20 20:18:23.92 20:34:08.35 21:42:21.71 22:00:04.88 21:42:01.22 02:42:39.67
RMS
Latitud
Longitud
0.65 0.73 0.75 0.97 0.68 0.70 0.40
37.0489 37.0358 37.0922 37.0932 37.1016 37.0723 37.1118
-2.5001 -2.5318 -2.5355 -2.5462 -2.5457 -2.5582 -2.5403
Profundidad (Km) 11.3 ± 0.5 0.5 ± 2.5 0.0f ± 0.0 1.1 ± 1.1 0.0 ± 1.3 1.2 ± 1.0 0.0f ± 3.7
TABLA 9 RELOCALIZACIÓN FINAL FAMILIA B. ID 1 2 3 4 5 6
Fecha 04 04 05 05 05 05
/ / / / / /
02 02 02 02 02 02
/2002 /2002 /2002 /2002 /2002 /2002
Hora origen (hh:mm:ss.ss) 20:44:32.02 21:37:54.00 05:28:20.33 06:08:10.76 08:40:11.86 12:51:06.21
RMS
Latitud
0.44 0.36 0.43 0.36 0.38 0.39
37.0742 37.0722 37.0659 37.0646 37.0660 37.0670
Longitud -2.5723 -2.5687 -2.5682 -2.5680 -2.5682 -2.5690
Profundidad (Km) 6.7 ± 2.1 7.6 ± 1.8 8.1 ± 1.8 7.7 ± 1.8 5.6 ± 7.0 7.3 ± 1.9
TABLA 10 LOCALIZACIÓN IGN FAMILIA B. ID 1 2 3 4 5 6
Fecha 04 04 05 05 05 05
/ / / / / /
02 02 02 02 02 02
/2002 /2002 /2002 /2002 /2002 /2002
Hora origen (hh:mm:ss.ss) 20:44:32.54 21:37:54.77 05:28:20.50 06:08:10.96 08:40:12.57 12:51:07.12
RMS
Latitud
Longitud
0.97 0.76 0.75 0.52 0.60 0.53
37.0463 36.9773 37.0622 37.0386 37.1054 37.0703
-2.5326 -2.5444 -2.5417 -2.5336 -2.5471 -2.5446
Profundidad (Km) 2.7 ± 2.5 6.0 ± 2.2 10.6 ± 0.3 16.1 ± 3.5 6.4 ± 2.7 3.8 ± 2.3
El número de correlaciones cruzadas de las que se dispone para cada sismo varía en cada caso (tabla 11).
TABLA 11 NÚMERO DE CORRELACIONES DE CADA SISMO Familia A Sismo Nº de correlaciones 1 49 2 34 3 48 4 53 5 53 6 48 7 28
Sismo 1 2 3 4 5 6
29
Familia B Nº de correlaciones 36 36 38 35 29 32 36
A partir de los resultados mostrados en las tablas 7, 8, 9, 10 y 11 se concluye que: 1. los valores en latitud y longitud son similares a los obtenidos por el IGN. No ocurre lo mismo con el valor de la profundidad. Para la familia A la profundidad de todos los terremotos está comprendida entre 8.8 y 10.4 Km, mientras que las profundidades del IGN tienen un mayor rango de variación (entre 1.1 y 11.3 Km) o no se ha podido obtener la profundidad (valores de 0.0f Km). En la familia B también se observa una clara mejora en la profundidad frente a las soluciones del IGN: todos los sismos de nuestra solución tienen profundidades comprendidas entre 5.6 y 8.1 Km, frente al rango de 2.7 a 16.1 Km del IGN, 2. los errores obtenidos en profundidad en la familia A son ligeramente mayores para el sismo 1 (1.4 Km frente a 0.5 Km del IGN) y similares en los sismos 4 y 6 (1.3 Km frente a 1.1 Km del IGN). En el resto de los sismos no es posible hacer esta comparación ya que para esos terremotos el IGN no determinó la profundidad. En la familia B los errores son menores salvo en los sismos 3 y 5, en que nuestros errores son sensiblemente mayores, 3. el valor cuadrático medio (RMS) de la hora origen en la familia A oscila entre (0.70 – 0.23) frente a (0.97 – 0.40) del IGN, lo que indica una clara disminución salvo para el sismo 1. En el caso de la familia B, se encuentra que el RMS oscila entre (0.36 - 0.44) frente a (0.52 - 0.97) de la solución del IGN, 4. los sismos para los que se dispone de un mayor número de retardos llevan asociado finalmente un error menor, de manera general. En la figura 12 se muestran los epicentros de la familia A (azul) obtenidos en este trabajo y los calculados por el IGN (verde). Se observa que, tras la relocalización, los epicentros se han concentrado agrupándose en una área reducida, siguiendo una línea de dirección de unos 52º con el norte geográfico y abarcando una extensión de unos 2 Km. La localización del IGN también mostraba una distribución NE-SW para la mayoría de los sismos pero en una área mucho más extensa. Las mayores diferencias corresponden a los sismos 1 y 2, muy alejados en la distribución del IGN, mientras que en la obtenida en este trabajo aparecen muy próximos a los restantes sismos. La figura 13 muestra la distribución de hipocentros en profundidad según un corte en dirección 52º N. Se observa que los focos se distribuyen entre los 8 y 10 Km de profundidad, siguiendo una línea con un buzamiento de 37º hacia el noreste.
30
Los epicentros de los sismos relocalizados de la familia B (azul) y los calculados por el IGN (verde) se muestran en la figura 14. Se observa una distribución de epicentros prácticamente N-S similar a la obtenida por el IGN, si bien los terremotos abarcan una extensión de unos 2 Km en este trabajo, frente a los cerca de 14 Km en el caso del IGN. En la figura 15 se muestra la distribución en profundidad de los focos de la familia B, para un corte N-S. Salvo el sismo 5 (con profundidad 5.6 Km), los demás se distribuyen a lo largo de una línea con un buzamiento de 56º hacia el sur. No obstante, este sismo es el que tiene un mayor error asociado, por lo que, si se tiene en cuenta este error, se puede considerar que el sismo sigue la tendencia general de un buzamiento de 56º hacia el sur. -2.60º
-2.55º
-2.50º
-2.45º
-2.40º
A Gergal 7
37.10º
5
7 5 2 3 4 6 1
A'
37.10º
Principal
4
3
6
37.05º
37.05º
1
2
37.00º
37.00º
36.95º -2.60º
-2.55º
-2.50º
-2.45º
36.95º -2.40º
Figura 12. Familia A. En verde se representan los epicentros de los sismos a partir de las coordenadas del IGN. En azul están representados los epicentros tras la relocalización con correlación cruzada. La estrella indica la localización del sismo principal, según la localización del IGN. Familia A
A'
A
0
Profundidad (km)
2 4 6 8 1 10 12 12
5 4 6 3
13
2 7 14
15
16
17
Figura 13. Corte en profundidad en dirección (37.07ºN, 2.06ºW) – (37.12º,2.52ºW) de los hipocentros de la familia A.
31
-2.60º
-2.55º
-2.50º
-2.45º
-2.40º
Gergal A
37.10º
5
37.10º Principal
1
3
2 6 5 4
6 3
37.05º
37.05º 1 4
37.00º
37.00º
2
36.95º -2.60º
-2.55º
-2.50º
36.95º -2.40º
-2.45º
Figura 14. Familia B. En verde se representan los epicentros de los sismos a partir de las coordenadas del IGN. En azul están representados los epicentros tras la relocalización con correlación cruzada. La estrella indica la localización del sismo principal, según la localización del IGN. Familia B
A
A'
0
Profundidad (km)
2 4 5
6
1 2
8
6 3
4
10 12
12
13
14
15
16
17
Figura 15. Corte en profundidad en dirección norte-sur de los hipocentros de la familia B. Cesca (2005) obtuvo para el sismo principal de Gérgal un mecanismo de falla normal con un plano de falla orientado en dirección NW-SE, con buzamiento hacia el NE y el foco a 9 Km de profundidad. El otro plano, orientado en dirección NE-SW, es prácticamente vertical. Por tanto no es posible relacionar los alineamientos de las familias A y B con el mecanismo del sismo principal. No obstante, hay que recordar que el sismo principal y sus réplicas no tienen por qué estar asociados a la misma falla, pudiendo ocurrir que el sismo principal desencadene el movimiento de fallas secundarias en las que se producen las réplicas.
32
Esta metodología ha sido aplicada también a la serie de réplicas del sismo de Alhucemas (Marruecos) del 24 de febrero de 2004 (mb=6.2). Sin embargo, para estos sismos, pese a que las réplicas son de mayor magnitud (4.0 – 5.0), no se han obtenido resultados satisfactorios ya que la estación más cercana se encontraba a 40 Km de distancia, estando las restantes a distancias entre 100 y 500 Km. Ya se ha visto que el método es muy sensible a poseer registros a distancias epicentrales pequeñas, lo que permite disponer de buenas formas de onda y en consecuencia poder realizar una buena correlación.
33
5. CONCLUSIONES La relocalización de terremotos por medio de formas de onda presenta claras ventajas frente a los métodos tradicionales, ya que a) permite disminuir los errores tanto de la profundidad como en la hora origen, al no utilizar los tiempos de llegada absolutos de la onda P, b) la menor influencia del modelo de Tierra permite una mejor determinación de la profundidad. Debido a la cercanía entre sí de los hipocentros, los rayos recorrerán trayectos similares hasta la estación, por lo que el modelo de Tierra no tiene una influencia directa sobre la relocalización. Los resultados de la aplicación de esta metodología a la serie de Gérgal (Almería) de 2002, formada por 39 sismos, ha permitido concluir que 1. es de gran importancia que los registros de cada sismo se hayan producido en las mismas estaciones para poder llevar a cabo el mayor número de correlaciones
posible,
con
el
fin
de
lograr
una
mejor
relocalización
especialmente en profundidad, 2. se han obtenido dos familias de terremotos formadas por 7 sismos (familia A) y 6 sismos (familia B) de entre los 39 sismos de la serie sísmica, 3. la familia A presenta una distribución de epicentros alineados según una dirección N52ºE, con un extensión de unos 2 Km. Los focos se distribuyen en profundidad entre 8 – 10 Km, siguiendo una línea de buzamiento de 37º hacia el noreste, 4. por su parte, los epicentros de la familia B se distribuyen según una línea con orientación N-S con una extensión de aproximadamente 1.5 Km. En profundidad, los hipocentros se sitúan entre los 5.6 Km y los 8 Km según una distribución con un buzamiento de 56º hacia al sur. 5. las soluciones de ambas familias presentan menores errores en profundidad y RMS frente a los valores calculados por el IGN. En particular, se ha logrado determinar la profundidad de varios sismos que el IGN no había determinado,
34
6. ninguna de las familias puede ser asociada al mecanismo focal propuesto por Cesca (2005) para el sismo principal de la serie de Gérgal, 7. esta metodología permite obtener detalles de la ocurrencia de series de réplicas, posibilitando el estudio más detallado del proceso de generación y ocurrencia de sismos, y 8. se ha empleado la transformada wavelet, frente a la transformada de Fourier, como viene siendo habitual, lo que permite un análisis simultáneo en el dominio del tiempo y de las frecuencias. Como futura línea de trabajo se propone: 1. la utilización de la onda S en la correlación cruzada, para obtener una relocalización más precisa al aumentar el número de correlaciones disponibles para cada par de sismos, y 2. una nueva línea de investigación, abierta por este trabajo, que va a permitir estudiar con mayor profundidad la relocalización tanto de series más numerosas de sismos, como enjambres de terremotos, con el fin de conocer las estructuras sismogénicas de diversas zonas de estudio, como es el caso particular del mar de Alborán.
35
6. BIBLIOGRAFÍA
Cesca, S. (2005). Inversión del tensor momento sísmico de terremotos superficiales a distancias regionales. Tesis doctoral. Universidad Complutense de Madrid. Madrid. 289 pp.
Dañobeitia, J.J., V. Sallares y J. Gallart. (1998). “Local earthquakes seismic tomography in the Betic Cordillera (southern Spain)”. Tectonophysics, 160, 225-239.
Deichmann, N. y M. García-Fernández. (1992). “Rupture geometry from highprecision relative hypocentre locations of microearthquake clusters”. Geophys. J. Int., 110 (3), 501-517.
Frechet, J. (1985). Sismogenèse et doublets sismiques. Tesis doctoral. Université Scientifique et Médicale de Grenoble. 206 pp.
Gaucher, E. (1998). Comportement Hydroméchanique d'un Massif Fracturé: Apport de la Microsismicité Induite. PhD thesis, Université de Paris VII.
Geiger, L.: Herdbestimmung bei Erdbeben aus den Ankunftszeiten. Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen 4 (1910), 331-349.
Geller, R. J. and S. Mueller. (1980). “Four similar earthquakes in central California”. Geophys. Res. Lett., 7(10), 821-824.
Goded Millán, T. (2006). Reevaluación de daños de los terremotos de Málaga de 1494 y 1680. Trabajo de investigación para la obtención del Diploma de Estudios Avanzados. Universidad Complutense de Madrid. 152 pp.
Goldstein, P. (1987). Seismic Analysis Code (SAC). Universidad de California.
Got, J.-L., J. Frechet, and F. Klein. (1994). “Deep fault plane geometry inferred from multiplet relative relocation beneath the south flank of Kilauea”. J. Geophys. Res., 99,15.375-15.386.
Goupillaud, P., A. Grossman, y J. Morlet. (1984). “Cycle-Octave and Related Transforms in Seismic Signal Analysis”. Geoexploration, 23:85-102.
36
Haukkson, E. y P. Shearer (2005). “Southern California hypocenter relocation with waveform cross-correlation, part 1: results using the double-difference method”. Bull. Seism. Soc. Am. 95, 3, 896–903.
Ito, A. (1985). “High resolution relative hypocenters of similar earthquakes by cross-spectral analysis method”. J. Phys. Earth, 33, 279–294.
Ito, A. (1990). “Earthquake swarm activity revealed from high-resolution relative hypocenters – clustering of microearthquakes”. Tectonophysics, 175, 47-66.
Klein,
F.W.
(1978).
Hypocenter
location
program
HYPOINVERSE,
U.S.
Geological Survey. Open-File Report 78-694, 113 pp.
Lee, W. H. K. y J. C. Lahr (1972). HYP071: A computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes, Open File Report, U. S. Geological Survey, 100 pp.
Levenstein, V. I. (1965) “Binary codes with use of deletions, insertions and substitutions of symbol”s. Dokl. A. N. SSSR, 163,845-848.
Lyon-Caen, H., P. Papadimitriou, A. Deschamps, P. Bernard, K. Makropoulos, F. Pacchiani y G. Patau (2004). “First results of the CRLN seismic array in the western Corinth Gulf: Evidence for old fault reactivation”. C. R. Geoscience, 336:343-351.
Menke, W. (1999). “Using waveform similarity to constrain earthquake locations”. Bull. Seism. Soc. Am. 89, 1143-1146.
Pacchiani, Francesco. (2006). Étude sismologique des failles normales actives du rift de Corinthe. Thèse de doctorat. Université Paris XI. 383 pp.
Paige, C.C. y M.A. Saunders (1982). « LSQR: Sparse linear equations and least squares problems”. ACM Transactions on mathematical software 8/2, 195-209.
Pavlis, G.L. (1986) “Appraising earthquake hypocenter location errors: a complete, practical approach for single-event locations”. Bull. Seism. Soc. Am., 76, 1699-1717.
37
Polikar, Robi. (1994). The wavelet tutorial. http://www-ccrma.stanford.edu/~unjung/mylec/WTpart1.html http://www-ccrma.stanford.edu/~unjung/mylec/WTpart2.html http://www-ccrma.stanford.edu/~unjung/mylec/WTpart3.html
Poupinet, G., W. L. Ellsworth, y J. Frechet (1984). “Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California”. J. Geophys. Res., 89(B7), 5719-5731.
Rodríguez Abad, I. (2004) Estudio de la determinación hipocentral e influencia de los errores. Trabajo de investigación para la obtención del Diploma de Estudios Avanzados. Universidad Complutense de Madrid. 133 pp.
Shearer, P., E. Hauksson y G. Lin (2005). “Southern California Hypocenter Relocation with Waveform Cross-Correlation, Part 2: Results Using SourceSpecific Station Terms and Cluster Analysis”. Bull. Seism. Soc. Am., 95, 904915.
Waldhauser, F. y W. L. Ellsworth. (2000). “A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault”. Bull. Seismol. Soc. Am., 90(6): 1353-1368.
Waldhauser, F. y W. L. Ellsworth. (2002). “Fault structure and mechanics of the Hayward Fault, California, from double-difference earthquake locations”. J. Geophys. Res., 107(B3), doi: 10.1029/2000JB000084.
Zhizhin, M., A. Gvishiani, S. Bottard, B. Mohammadioun y J. Bonnin (1992). “Classification of strong motion waveforms from different geological regions using syntactic pattern recognition scheme”. Cah. C. Europ. Geodyn. Sismol., 6, 33-42.
38
7. ANEXO 1
Cuentas digitales
Ampliación de la figura 4.
t (s)
39