ARTÍCULOS ORIGINALES

Análisis comparativo de diferentes métodos de picado automático de fases en terremotos registrados en la Estación Sismológica de La Plata (LPA)

 

Juan I. Sabbione^1,2, María Laura Rosa¹, Danilo R. Velis^1,2 y Nora C. Sabbione¹

¹ Facultad de Ciencias Astronómicas y Geofísica- Universidad Nacional de La Plata.
² CONICET
jsabbione@fcaglp.unlp.edu.ar*, mlrosa@fcaglp.unlp.edu.ar, velis@fcaglp.unlp.edu.ar, nora@fcaglp.unlp.edu.ar

 

---

RESUMEN

Analizamos ocho métodos de detección y picado automático de fases internas de terremotos. Entre ellos, tres son métodos convencionales, dos son métodos nuevos desarrollados a partir de mejoras a los métodos convencionales, y los otros tres son adaptaciones a métodos originalmente concebidos para picado de primeros arribos en datos de sísmica de exploración. Para analizar los métodos, seleccionamos once registros de una hora de duración obtenidos en la Estación Sismológica de la Plata (LPA) correspondientes a terremotos regionales y telesísmicos con magnitudes mayores a 6 (Mw) y diversas distancias epicentrales y profundidades de foco. Los registros elegidos fueron previamente picados manualmente por un analista. Se buscó picar de manera automática un total de 30 fases, de las cuales 25 habían sido originalmente observadas por el analista. Las otras 5 fases, fueron detectadas a partir de la aplicación de los métodos automáticos, y luego corroboradas mediante el modelo IASP91. Se estudió la eficiencia de los métodos a partir del porcentaje de detecciones obtenidas y del porcentaje de picados falsos generados. Se analizaron también los errores en los picados desde un enfoque estadístico para estimar las exactitudes y precisiones de los métodos propuestos. Este análisis permitió realizar una valoración relativa de los ocho métodos de picado automático propuestos. Los resultados obtenidos revelan que algunos de los métodos exhiben ventajas significativas por sobre otros. Estos métodos que muestran mejor desempeño serían apropiados para realizar sistemáticamente el picado automático, tanto para estudios de sismología global como para estudios de sismicidad, especialmente cuando se cuenta con un gran volumen de datos."

Palabras clave: Detección de eventos; Fases de terremotos; Picado automático; Falsas alarmas.

ABSTRACT

We analyzed eight methods of detection and automatic body phase picking. Three of them are conventional methods, two are new methods developed from improvements to conventional methods, and the other three are methods originally designed for first-break picking of seismic exploration data. To analyze the methods, we selected eleven one-hour records obtained at La Plata Seismological Station (LPA) corresponding to regional and teleseismic earthquakes with magnitudes greater than 6 (Mw) and diverse epicentral distances and focus depths. The records were previously manually picked by an analyst. We attempted to pick a total of 30 phases automatically, 25 of which had originally been observed by the analyst. The other 5 phases, were identified from the application of the automatic methods, and then corroborated by using the IASP91 model. The efficiency of the methods was studied based on the detection percentage obtained and the false phase percentage generated. The picking errors were also analyzed from a statistical point of view in order to estimate the accuracy and precision of the proposed methods. This analysis allowed us to make a relative assessment of the eight automatic picking methods proposed. The results show that some of the methods exhibit significant advantages over the others. These methods that show better performance would be appropriate to systematically carry out the automatic picking, both for global seismology studies and seismicity studies, especially when a large volume of data is available.

Keywords: Event detection; Earthquake phases; Automatic picking; False alarms.

---

 

INTRODUCCIÓN

En sismología global, la detección y picado de fases de eventos representa un paso fundamental para la localización y análisis de datos de terremotos en cualquier estación sismológica. Tradicionalmente, esta tarea era realizada por inspección visual mediante la lectura de las bandas de los registros de las estaciones.
Actualmente, la gran mayoría de los instrumentos utilizados en las estaciones sismológicas digitalizan y almacenan los datos. Al mismo tiempo, el monitoreo sísmico se ha multiplicado significativamente en todo el planeta, lo que genera la necesidad de analizar un inmenso volumen de datos. Esta abundancia de datos produjo que la inspección visual de los registros se vuelva una tarea muy demandante en cuanto a tiempo y recurso humano. Por esta razón, en las últimas décadas se han desarrollado diversos algoritmos para realizar la detección y picado de fases de terremotos de manera automática o semiautomática. Muchos de estos algoritmos están actualmente incluidos en paquetes de procesamiento utilizados por diferentes organizaciones o instituciones.
Por otro lado, cabe destacar que puede trazarse una analogía entre el problema que representa la detección de fases de terremotos en sismología y el picado automático de primeros arribos en sísmica de exploración, proceso fundamental en la industria petrolera. Claramente, la semejanza entre estos dos problemas se debe principalmente a que los fundamentos físicos que subyacen en ambas teorías son muy similares. Esta semejanza queda también evidenciada en la literatura por la presencia de métodos que usan algoritmos muy similares para realizar una tarea o la otra.
Entre los artículos presentados en revistas de sismología son muy populares los métodos que comparan el valor de algún atributo de la señal en dos ventanas de tiempo diferentes. Este atributo es calculado en una primera ventana de corta duración en tiempo que intenta extraer las características actuales de la señal, y en otra de mayor duración que tiene por objeto dar la referencia del ruido de fondo en el que está enmascarado el evento a detectar. Allen (1982) introduce un nuevo concepto definiendo dicho atributo como "characteristic function" (CF). Debido a que estos métodos suelen comparar los promedios de dicha CF entre las dos ventanas de tiempo, suelen ser denominados métodos de tipo STA/LTA ("short term average / long term average"). El evento se suele declarar cuando este cociente supera un valor umbral THR ("threshold"). El método que fue presentado originalmente en Allen (1978) ha sido ampliamente utilizado por el U.S. Geological Service (USGS), en particular en el proyecto "USGS Earthworm system" (Johnson et al., 1995).
Además del método de Allen, existen otros métodos del tipo STA/LTA muy difundidos. Uno de los más citados es aquel presentado por Baer y Kradolfer (1987), que es actualmente utilizado como el motor del software de picado "Manneken Pix" desarrollado por Aldersons (2004). La rutina AUTOPIC contenida en el Seismic Analysis Code (SAC) utiliza el método presentado por Ruud y Husebye (1992), también de tipo STA/LTA. Por su parte, Earle y Shearer (1994) definen la envolvente analítica de la señal como la CF para calcular el STA/LTA en su método, el cual utilizan para hacer un estudio de sismología global a partir de la construcción de tablas de tiempos de recorrido. Todos estos métodos del tipo STA/LTA, que podríamos denominar "convencionales" por su amplia difusión y aceptación a lo largo del tiempo, guardan ciertas similitudes con el tradicional método para picado de primeros arribos sísmicos desarrollado por Coppens (1985), en donde se introduce un criterio de comparación de energías entre dos ventanas anidadas de diferentes tamaños.
Se han desarrollado también diversos algoritmos de picado no convencionales que utilizan modelos autorregresivos (AR) junto con el Criterio de Información de Akaike (AIC) presentado en Akaike (1971) y Akaike (1974). Los trabajos más difundidos con este enfoque son Sleeman y van Eck (1998), Leonard y Kennett (1999), Bai y Kennett (2000), y Zhang et al. (2003), entre otros.
La teoría de los fractales fue también utilizada tanto en sismología para picar fases de terremotos (Tosi et al., 1999), como en sísmica de exploración para picar primeros arribos (Boschetti et al., 1996; Jiao y Moon, 2000; Sabbione y Velis, 2010). Recientemente, se han realizado también algunos intentos para abordar el mismo problema a partir del uso de estadística de órdenes superiores (Saragiotis et al., 2004; Küperkoch et al., 2010), en los que nuevamente se puede encontrar una correlación con métodos de exploración sísmica (Yung e Ikelle, 1997).
No obstante, se debe remarcar que existe una diferencia substancial entre el picado de primeros arribos sísmicos y el picado de fases de terremotos. Mientras que en el problema relativo a la exploración sísmica se busca detectar sólo un evento (el primer arribo) por cada registro (la traza sísmica), en los registros de estaciones sismológicas arriban diferentes fases de terremotos y por ello los métodos intentan picar varios tiempos. Así, mientras que en los métodos símicos se busca un único rasgo distintivo a lo largo de toda la traza, los algoritmos utilizados en sismología tienden a detectar varias fases cuando se supera un umbral o tolerancia dada. Por esta razón, si se pretende utilizar algoritmos originalmente desarrollados para picar primeros arribos sísmicos, éstos deben ser previamente adaptados para identificar fases de terremotos en registros de estaciones sismológicas.
Uno de los mayores problemas de los métodos de picado automático de fases de terremotos es el de los "picados falsos". En este trabajo vamos a diferenciar dos clases de picados falsos. Definimos como "falsas alarmas" a aquellas detecciones de algún evento que no están relacionadas a un sismo (e.g., detecciones que ocurren antes de la llegada de las ondas primarias) o mucho después de la coda, y llamaremos "falsas fases" a aquellas fases picadas luego de la llegada de las fases primarias de un sismo, pero que no logramos asociar a ninguna fase teórica. Esta distinción cobra importancia para aquellas estaciones donde se realizan estudios de sismicidad.
En dichos estudios se intenta acotar el número de las falsas alarmas, procurando paralelamente no perder sensibilidad en la detección de eventos reales. Para definir entonces los picados falsos, dada una fase detectada mediante algún método automático que no se corresponda con el tiempo de alguna de las fases obtenidas mediante la lectura de bandas, se debe corroborar primero que el tiempo detectado no coincida (aproximadamente, con cierta tolerancia) con el tiempo de arribo de alguna fase teórica no distinguida por el analista al leer la banda. Si esto efectivamente no ocurre, dicha detección automática representa un picado falso.
La mayoría de los métodos nuevos que aparecen en la literatura suelen ser validados mediante la comparación contra métodos convencionales, pero no contra métodos menos difundidos o más recientes. Se han publicado también diversos trabajos en donde, además de presentar eventualmente nuevas estrategias, se realizan algunas recopilaciones y comparaciones de los métodos convencionales preexistentes (Allen, 1982: Withers et al., 1998; Leonard, 2000). Es claro que la performance de un dado método está ligada no sólo a la calidad de los datos (generalmente caracterizada por la relación señal - ruido), sino también a las diferentes características de los mismos. Desafortunadamente, resulta dificultoso hallar trabajos en los que se comparen y analicen diferentes métodos en función de las características del dato inspeccionado de manera de servir de guía a un sismólogo o analista a la hora de discernir qué método de picado automático utilizar para procesar los datos con los que cuenta.
Con la motivación de lo recién expuesto, hemos seleccionado once registros de terremotos adquiridos en la Estación Sismológica de La Plata (LPA). Los registros contienen terremotos regionales y telesísmicos con distancias epicentrales mayores a 1000 Kms y magnitudes mayores a 6 (Mw). Hemos utilizado este conjunto de datos para analizar y comparar ocho métodos de picado de fases de terremotos, incluyendo métodos convencionales, no convencionales, y variaciones de métodos preexistentes propuestas en este trabajo. Entre ellos, testeamos tres métodos convencionales: el método presentado en Allen (1978), aquel introducido por Baer y Kradolfer (1987), y el método de Earle y Shearer, (1994); dos métodos basados en modificaciones a estos mismos métodos (que llamaremos método de Allen modificado y método de Baer y Kradolfer modificado); y por último tres adaptaciones a métodos que han resultado ser muy útiles en el picado de primeros arribos en exploración sísmica: el método de Coppens modificado, el método de la entropía, y el método de la dimensión fractal, todos ellos presentados en Sabbione y Velis (2010).
El presente trabajo tiene también por objetivo brindar herramientas a la Estación Sismológica de La Plata que permitan la detección de fases de terremotos de manera automática. De esta manera, se pueden aprovechar la consistencia y la rapidez que se obtienen al realizar esta tarea de forma automática. Paralelamente, el conjunto de métodos de picado de fases generado podría utilizarse para realizar estudios de sismicidad en regiones de interés, que suelen involucrar un gran volumen de datos, y para los que la detección automática de eventos resulta primordial. En particular, este sería el caso del sistema de Fallas Magallanes-Fagnano en la provincia de Tierra del Fuego, donde el Departamento de Sismología de La Plata cuenta con una serie de estaciones permanentes y personal que realiza monitoreo continuo de los registros de los mismas (Sabbione et al., 2007; Buffoni et al., 2009).
Los resultados obtenidos con los ocho métodos de picado automático son analizados a partir de diferentes enfoques. Se estudia primero la capacidad de los métodos de detectar las fases leídas en las bandas. Uno de los resultados más interesantes que arroja la aplicación de los métodos es la detección de algunas fases reales originalmente omitidas por el analista. Además, para complementar el análisis de las detecciones, se calcula también el número de picados falsos obtenidos. Este análisis, junto con el estudio estadístico de los errores en los picados, nos permite recomendar cuáles de los métodos aquí presentados serían más aptos para realizar el picado automático de fases de terremotos de forma sistemática en estudios sismológicos de gran escala, o bien para realizar estudios continuos de sismicidad, que requieren el procesamiento de un gran volumen de datos.

MÉTODOS

En esta sección realizamos una muy breve descripción de los métodos utilizados en el presente trabajo. Para un análisis más profundo de los mismos, se citan las referencias de aquellos artículos en donde los métodos fueron presentados. En algunos casos, describimos las mejoras propuestas en forma detallada.
Cualquier método de picado automático requiere que el usuario seleccione y fije un determinado número de parámetros para su aplicación. Esta selección representa un punto crucial en la detección de fases de eventos mediante algoritmos computacionales. Si los parámetros no son elegidos de manera adecuada, los métodos tenderán a fallar. Asimismo, si los resultados que se obtienen fueran muy dependientes de la selección de los parámetros, se deberá tener especial cuidado en su selección. Por ello, cabe resaltar que los métodos más robustos serán aquellos para los que se precise seleccionar el menor número de parámetros y cuyos resultados sean poco dependientes de dicha selección.
Para fijar los parámetros de los métodos se suelen considerar características de amplitud y tiempo de la señal, su contenido espectral, relación señal - ruido, etc. Por su parte, estas características están relacionadas con la distancia epicentral y profundidad del foco del terremoto, la magnitud, la geología local, el ruido de la estación, las características del instrumento, etc. Dado que La Plata es una región con sismicidad local moderada a baja, los datos seleccionados corresponden a terremotos regionales y telesísmicos de diferentes características. Uno de los mayores desafíos en este trabajo se debe justamente a que al ser los datos tan disímiles, la selección de los parámetros para cada método resulta muy delicada. Al final de esta sección se incluye una tabla (Tabla 1) donde se resumen los valores para los distintos parámetros de cada método descrito a continuación.

Tabla 1. Parámetros utilizados para el picado automático de fases.

Antes de comenzar con la descripción particular de cada método, debemos mencionar que el primer paso de todos los algoritmos consiste en remover la componente continua de la señal y normalizar la amplitud de los registros a . De esta manera homogeneizamos la escala vertical del dato de entrada en todos los métodos y facilitamos la elección de los parámetros.

Método de Allen (RAM)

El método presentado en Allen (1978), que llamaremos RAM (Rex Allen method), representó un gran avance en el picado de fases de terremotos mediante algoritmos computacionales. La función característica CF utilizada combina atributos de amplitud y fase de la señal s_i, donde representa el número de muestras:

CF_i = s_i² + C_i • (s_i - s_i-1)²           (1)

donde C_i es un factor de peso que vincula en la función característica los dos términos que la definen, uno relacionado con la amplitud de la señal, y el otro relacionado con su frecuencia. Tanto Baer y Kradolfer (1987) como Küperkoch et al. (2010) sugieren elegir:

   (2)

de manera de darle peso similar al término que depende de la amplitud y al que depende de la frecuencia.

A continuación, se calcula el "short term average" (STA) y el "long term average" (LTA) mediante:

(3)

donde C₃ y C₄ son constantes que controlan el tamaño de las ventanas de los promedios STA y LTA. El algoritmo declara un evento cuando:

STA_i >THR • LTA_i,       (4)

donde THR es un umbral ("threshold") que determina la sensibilidad del método. Finalmente, se define un "criterio de continuación" para determinar cuándo el evento ha finalizado (Allen, 1978).

Método de Baer y Kradolfer (BKM)

Baer y Kradolfer (1987) introducen la función E²_i, que es similar a la función característica utilizada en el RAM:

       (5)

Además, exponen motivos para modificar el criterio STA/LTA y en su lugar realizan el picado a partir de la variable estadística normalizada de E_i⁴:

       (6)

Cabe mencionar aquí que en la expresión presentada por Baer y Kradolfer (1987) se utiliza la varianza de E_i⁴ en el denominador de la ec. (6) en lugar de su dispersión. Hemos introducido esta pequeña modificación en la expresión original pues obtuvimos mejores resultados con la inclusión de este cambio. Finalmente, en el método de Baer y Kradolfer (BKM), los eventos se declaran si BK_i supera el umbral THR y sólo son aceptados si esta situación se mantiene por un tiempo mayor a T_up y no se interrumpe por un tiempo mayor a T_down, caso contrario el evento se rechaza. El objetivo de T_up y T_down es evitar tanto las falsas alarmas por cortos incrementos en BK_i como rechazar eventos reales por breves disminuciones abruptas, respectivamente.

Método de Earle y Shearer (ESM)

En el último método "convencional" que presentamos, Earle y Shearer (1994) utilizan la envolvente de la señal como función característica:

        (7)

donde es la transformada de Hilbert de s_i.
Luego, se definen dos ventanas de tiempo consecutivas de tamaño T_STA y T_LTA para el cálculo de STA y LTA y se calcula su cociente muestra a muestra, asignándolo a la primera muestra de la ventana más adelantada en tiempo, que corresponde a STA. Para suavizar este cociente y así prevenir las rápidas fluctuaciones, en lugar de utilizar los parámetros T_up y T_down como en el BKM, Earle y Shearer (1994) optaron por aplicar un filtro pasabajos de Hanning de longitud T_Han. Esta idea resulta muy útil y la utilizaremos más adelante en los nuevos métodos.
En el método de Earle y Shearer (ESM) los eventos se declaran cuando el cociente suavizado supera un dado umbral THR. La fase es picada en el punto de inflexión inmediatamente anterior el máximo de este cociente suavizado.

Método de Coppens modificado (MCM)

En el método de Coppens modificado (MCM), presentado en Sabbione y Velis (2010), la función característica utilizada es la energía de la señal. Si bien no es estrictamente un método del tipo STA/LTA porque no se calculan promedios, el esquema es muy similar. Se definen dos ventanas móviles anidadas, una de longitud T_C1 que intenta captar el carácter presente de la señal, y una segunda ventana, que a diferencia del trabajo original de Coppens (1985) hemos c onsiderado de longitud fija T_C2, y que pretende caracterizar el ruido de fondo del dato. Dicho cociente de energías, ER_i, es entonces calculado muestra a muestra mediante:

     (8)

donde N_C1 y N_C2 son el número de muestras que corresponden a las ventanas T_C1 y T_C2. Claramente, ante el arribo de una señal impulsiva, este cociente de energías presentará un rápido aumento que se acercará al valor 1. Vale aclarar que el cociente de energías es asignado a la última muestra de las ventanas.
De la misma manera que en Sabbione y Velis (2010), aplicamos a ER_i un filtro EPS ("edge preserving smoothing") (Luo et al., 2002) de longitud en tiempo T_EPS para resaltar los cambios bruscos. Finalmente, las fases son picadas en aquellas muestras donde el valor de la derivada de la resultante del filtro supera el umbral THR.

Método de la Entropía (EM)

En el método de la Entropía (EM) se introduce la entropía de la serie de tiempo (Denis and Crémoux, 2002) como un atributo sísmico. Esta entropía es calculada mediante:

          (9)

donde L_i es la "longitud" de la curva al tiempo t_i. Dado que la entropía se calcula en una ventana de longitud temporal T_H (que equivale, digamos, a N_H muestras) que se mueve muestra a muestra a lo largo de la traza, la expresión de la entropía resulta:

            (10)

     Este atributo tiende a aumentar rápidamente ante el advenimiento de la señal en un contexto ruidoso. Aquí hemos introducido una leve variación respecto del método publicado en Sabbione y Velis (2010). Obtuvimos mejores resultados filtrando la entropía H_i mediante la convolución con una ventana de Hanning de longitud T_Han de la misma manera que en el ESM. Luego, al igual que en el MCM, aplicamos un filtro EPS, y tomamos su derivada para picar las fases en aquellas muestras que superen un dado umbral THR.

Método de la dimensión fractal (FDM)

El método de la dimensión fractal (FDM) es otro método que resultó ser muy útil para picar primeros arribos sísmicos. Los detalles de la versión aquí usada del algoritmo pueden encontrarse en Sabbione y Velis (2010). El método se basa en la idea de considerar la serie de tiempo del dato como una curva fractal en el plano (Mandelbrot, 1983). Como tal, se puede extraer su dimensión fractal mediante diferentes métodos (Turcotte, 1997). La dimensión fractal teórica de una serie de tiempo correspondiente a ruido blanco tiende al valor 2, mientras que se aproxima a 1 en la presencia de señal. Aquí se calcula la dimensión fractal de la curva a partir del "método del variograma" (Turcotte, 1997). El método se resume a continuación.
El primer paso consiste en sumar ruido blanco de muy baja amplitud, definido por su relación señal-ruido (SNR) a todo el dato. El objetivo de esto es descorrelacionar el ruido de baja amplitud presente en la estación para así diferenciar mejor su dimensión fractal de aquella que corresponde a la señal de interés. La selección de la relación señal - ruido SNR sumada al dato es fundamental para que el método funcione correctamente. Este es un aspecto delicado del algoritmo cuyo efecto ya ha sido observado por Jiao y Moon (2000) y que se encuentra explicado en Sabbione y Velis (2010).
Una vez sumado el ruido, el variograma es calculado dentro de una ventana móvil de longitud T_F y N_F muestras mediante:

       (11)

para cuatro retardos diferentes h=1, 2, 3, y 4.
La teoría de los fractales (Korvin, 1992) supone que existe una relación exponencial entre h y V_i(h) de la forma:

      (12)

donde D_i es la dimensión fractal de la curva (en la ventana considerada). Entonces, tomando logaritmo a ambos lado en la ec. (12), la relación se vuelve lineal y los pares (log[h], log[V_i(h)]) obtenidos para los 4 retardos h se alinean (aproximadamente) en una recta. Ajustando dicha recta por mínimos cuadrados se obtiene su pendiente, que conforme la ec. (12) será igual a , y así se despeja finalmente la dimensión fractal D_i de la ventana móvil considerada. Desplazando la ventana muestra a muestra, se obtiene D_i punto a punto.
Luego, al igual que en el MCM y el EM, se aplica el filtro EPS (en este caso a la dimensión fractal). Por último, teniendo en cuenta que se espera un decrecimiento en D_i ante el arribo de la señal, se pican las fases en aquellas muestras en donde la derivada cambiada de signo de la resultante del filtro supere un dado umbral THR.

Método de Allen modificado (MAM)

A partir de haber hecho algunas pruebas con los distintos métodos aquí presentados hemos incorporado algunos cambios al algoritmo propuesto en Allen (1978) y hemos desarrollado un método similar que denominamos método de Allen modificado (MAM). La función característica elegida es la misma que en el método original, y está dada por las ecs. (1) y (2). A diferencia del método original, aquí usamos ventanas móviles de longitud T_STA y T_LTA para calcular STA y LTA, respectivamente. Dado un intervalo de muestreo constante estas longitudes equivalen a, digamos, N_STA y N_LTA muestras. Así, la función cociente RF entre STA y LTA se escribe, para la muestra i:

         (13)

Este esquema de ventanas es el mismo al utilizado por Earle y Shearer (1994) y se ilustra en la Figura 1 (b). Diferenciándonos nuevamente del RAM, donde se utiliza un criterio de continuación para decir cuándo un determinado evento ha finalizado, en el MAM suavizamos la función cociente mediante una ventana de Hanning (Figura 1(c)). Luego, el evento es declarado cuando dicho cociente suavizado (SRF) supera un dado umbral THR, y se lo considera finalizado cuando esto deja de ocurrir. Las fases se pican en los máximos locales de SRF donde el umbral es superado (Figura 1).

Figura 1. Método de Allen modificado (MAM). (a) Registro de la componente vertical de LPA. El picado automático obtenido mediante el MAM se muestra con líneas verticales grises. (b) Función característica CF y esquema de ventanas para STA y LTA del MAM. (c) STA/LTA suavizado por una ventana de Hanning y umbral THR utilizado para detectar las fases picadas (líneas verticales grises).

Método de Baer y Kradolfer modificado (MBKM)

El último método que presentamos aquí es el método de Baer y Kradolfer modificado (MBKM). Este método es esencialmente el mismo que el MBK, pero con una diferencia fundamental. En lugar de utilizar T_up y T_down para evitar las rápidas fluctuaciones y así reducir las falsas alarmas y los falsos rechazos de eventos, suavizamos nuevamente con una ventana de Hanning de longitud T_Han (Figura 2). La gran ventaja de este cambio aparentemente menor, pero clave, es la reducción del número de parámetros a fijar para correr el método. De esta manera, sólo dos parámetros deben ser seleccionados para aplicar el MBKM (ver Tabla 1). Aclaramos por último que en el MBKM marcamos el arribo de las fases en la primera muestra para la cual se supera el umbral THR, como se ve en la Figura 2 (c).

Figura 2. Método de Baer y Kradolfer modificado (MBKM). (a) Registro de la componente vertical de LPA. El picado automático obtenido mediante el MBKM se muestra con líneas verticales grises. (b) Función característica CF del MBKM. (c) CF suavizada por una ventana de Hanning y umbral THR utilizado para detectar las fases picadas (líneas verticales grises).

Tabla 2. Características de los terremotos seleccionados y las 30 fases utilizadas para analizar los métodos. Las 5 fases que se muestran en negrita cursiva son aquellas sólo detectadas por algunos de los métodos, y cuyos arribos teóricos fueron calculados mediante el modelo IASP91.

Selección de parámetros

Para completar la descripción de los métodos, enumeramos y describimos todos los parámetros utilizados en la Tabla 1. Los parámetros están enumerados en el orden en que son utilizados en el flujo de los algoritmos que implementan los distintos métodos. Notar que todos los algoritmos precisan de la elección del umbral THR, parámetro que controlará la sensibilidad del método. Del mismo modo, muchos algoritmos comparten otros parámetros, pero sus valores no son necesariamente coincidentes. El único parámetro que no fue fijado igual para todos los datos analizados es SNR en el FDM, dado que este método requiere que dicho parámetro sea específica y cuidadosamente elegido para cada dato analizado.

DATOS

En la Estación Sismológica de La Plata (LPA) se obtienen registros analógicos con un instrumento de largo período (Sprengnether) y registros digitales con un instrumento de banda intermedia (Guralp CMG 3ESPC). En la actualidad, un analista realiza la lectura de las bandas analógicas obtenidas con el instrumento Sprengnether para obtener los tiempos de arribos de las diferentes fases ante un dado evento. La información se transmite en forma codificada al ISC (International Seismological Center, Gran Bretaña) y al NEIC-USGS (National Earthquake International Center - United States Geological Service, Estados Unidos).
La Estación Sismológica de La Plata se caracteriza por realizar una exhaustiva lectura de fases secundarias. Las fases más importantes de estas lecturas correspondientes a once registros de terremotos fueron las utilizadas para poder evaluar los resultados del picado de los métodos aquí presentados. Los eventos utilizados para analizar los métodos corresponden a terremotos regionales y telesísmicos de magnitud mayor a 6 (Mw) y con focos ubicados a diferentes distancias epicentrales y profundidades. Seleccionamos la componente vertical de los registros en donde a priori deberían observarse fases más claras y con mejor relación señal/ruido. Para poder contar con los registros digitales correspondientes a estas bandas, se utilizaron los datos registrados en la componente vertical del sensor Guralp (CMG-3ESPC) y se los pre-procesó removiendo su respuesta de instrumento y agregando la respuesta del Sprengnether de manera de simular las bandas leídas. Este pre-proceso fue realizado mediante el Seismic Analysis Code (SAC). Teniendo en cuenta que el Sprengnether es un instrumento de largo período, el intervalo de muestro de los datos digitales generados fue fijado en 0,5s.
Como se mencionó anteriormente, la amplitud de cada registro fue normalizada a (-1,1) para preparar los datos de entrada para los métodos de picado automático. Hemos aislado registros horarios para cada sismo. Los once registros utilizados para el análisis de los métodos se observan en la Figura 3. En la Figura 3 y en la Tabla 2 que se presenta a continuación, puede apreciarse que los terremotos fueron seleccionados con diferente relación señal/ruido, como así también con distintas profundidades, distancias epicentrales, magnitudes, etc., como mencionamos anteriormente.
En la Tabla 2 se resumen las características de los terremotos y las fases a picar por los métodos. El conjunto de las 30 fases que se muestran en la tabla, está compuesto por 25 fases picadas por el analista en la lectura de bandas y por otras 5 fases (en negrita cursiva en la Tabla 2), picadas por algunos de los métodos, pero que no fueron originalmente detectadas por el analista. Para comprobar que estas últimas fueran reales, se realizó el cálculo de los tiempos de arribo teórico utilizando el modelo de Tierra IASP91 (Kennett, 1991). Luego, una vez realizada esta comprobación, agregamos dichas fases originalmente no distinguidas al conjunto de fases utilizadas para analizar los métodos. Buscamos entonces picar de manera automática las 30 fases que se muestran en la Tabla 2. Notar que la gran mayoría de estas fases corresponden a ondas que arriban a la estación como ondas compresionales (ondas P), que son aquellas que esperamos detectar con más claridad en la componente vertical de la estación, en especial para eventos regionales y telesísmicos como los estudiados.

Figura 3. Registros en la componente vertical de LPA de los 11 terremotos utilizados para analizar los métodos. En las abscisas se muestra el tiempo (UTC) y en las ordenadas la amplitud normalizada a .

RESULTADOS

Para evaluar los resultados de los métodos se consideraron las fases de la Tabla 2. No se realizó ningún tipo de restricción sobre el número de fases a picar. Esta última es una hipótesis bastante fuerte, ya que produce un mayor riesgo de picados falsos. Dado que los algoritmos no son complejos, los métodos corren muy rápido y una vez fijados los parámetros (Tabla 1) son muy simples de aplicar.
El primer resultado que se debe resaltar se refiere a las falsas alarmas (detección errónea de fases antes de la llegada de las ondas primarias, o posterior a la finalización de la coda). Se obtuvieron solamente un total de cuatro falsas alarmas con el método de Baer y Kradolfer (BKM) en dos registros diferentes, una falsa alarma con el método de la entropía (EM), y ninguna con los otros seis métodos analizados. Este bajo número de falsas alarmas muestra que los parámetros correspondientes fueron seleccionados de manera de no hacer los métodos demasiado sensibles, lo que a su vez brinda mayor confiabilidad respecto de los picados obtenidos.
Por otro lado, como puede apreciarse en la Figura 3, las fases de los dos terremotos ocurridos en Argentina (los dos más cercanos), es decir, las fases P y S del terremoto de Jujuy del 18 de noviembre de 2007, y las fases P y sP del terremoto de Santiago del 3 de septiembre de 2008, son muy claras y fácilmente detectables. Para estos dos eventos, como se esperaba, todos los métodos detectaron sus dos fases, sin ningún picado falso, con una sola excepción dada por el BKM, que en el terremoto de 2007 generó un único picado falso. Los resultados obtenidos para los registros de estos dos sismos muestran que los parámetros fueron seleccionados de manera tal de asegurar que se detecten correctamente las fases claras.
Para evaluar la eficiencia de los métodos calculamos el porcentaje de detecciones acertadas frente al total de las 30 fases a picar. Junto con este valor, calculamos también el porcentaje de picados falsos generados frente al total de picados realizados por cada método. Si bien debemos mencionar que el número total de fases analizadas no es muy grande, la muestra es suficiente para realizar el presente análisis cualitativo de los métodos. Estos resultados, resumidos en la Figura 4, nos indican la eficiencia de los métodos, puesto que brindan información sobre qué tan adecuados resultan para detectar fases reales sin generar picados falsos.

Figura 4. Eficiencia de los métodos. Las barras negras muestran el porcentaje de detecciones correctas y las barras grises el porcentaje de picados falsos para cada método.

En primer lugar, vemos que a excepción del BKM, los métodos presentan un porcentaje de picados falsos que no supera el 32%. Esta cota responde a que los parámetros fueron ajustados con el mismo criterio para todos los métodos: obtener un porcentaje alto de detecciones, pero intentando al mismo tiempo acotar los picados falsos.
Podemos ver también en la Figura 4 que el RAM no realiza ningún picado falso. Esto se debe principalmente a que la CF que utiliza el método, dada por la ecs. (1) y (2), es muy apropiada para detectar la llegada de señal, incluso en entornos ruidosos. El RAM resulta entonces un método muy fiable cuando detecta algún evento. No obstante, el porcentaje total de detecciones obtenido con este método es uno de los más bajos (menos del 50% de las fases a picar). Esto ocurre a causa del complicado "criterio de continuación" del método, ya que una vez detectado un evento, al algoritmo suele interpretar que continúa en la coda y no se pican fases posteriores a la P. El criterio de continuación del RAM debería ser eventualmente adaptado en diferentes estaciones sismológicas. Esta es la motivación por la cual aquí presentamos el método de Allen modificado (MAM).
Contrariamente al RAM, el BKM, que es el segundo método convencional aquí analizado, es el método que mayor porcentaje de picados falsos genera. Sin embargo, cabe destacar que presenta un porcentaje de detecciones relativamente bueno (cerca del 67%). El elevado porcentaje de picados falsos se debe al criterio adoptado en Baer y Kradolfer (1987) para decidir cuándo validar un evento y cuándo descartarlo (a partir de los parámetros T_up y T_down). Por esta razón es que intentamos en este trabajo mejorar dicho método e introducir el Método de Baer y Kradolfer modificado (MBKM), como ya fue mencionado.
Se destaca también en la Figura 4 el alto porcentaje de fases detectadas por el ESM y por el MAM, cercano al 90% (en promedio pican 27 de las 30 fases analizadas). El ESM y el MAM resultan ser los más eficientes a la hora de detectar fases de terremotos. Además, los dos métodos generan un porcentaje relativamente pequeño de picados falsos (15%, aproximadamente). Contrariamente, y a pesar del intento de mejora introducido, la performance del MBKM no se destaca por sobre la media de los otros métodos analizados.
Por último, y respecto de los tres métodos desarrollados originalmente para detectar primeros arribos en procesamiento de exploración sísmica (MCM, EM, y FDM), observamos que el MCM y el FDM funcionan ligeramente mejor que el EM. Ambos pican algo más que el 60% de las fases, y además generan alrededor de 10% de picados falsos solamente. Sin embargo, esta performance es inferior a la del ESM y a la del MAM, especialmente en cuanto al porcentaje de detecciones correctas. Interpretamos estos resultados como relativos a la dificultad de adaptar métodos de picado de primeros arribos sísmicos para detectar fases de terremotos.
Como ya hemos mencionado, algunos de los métodos arrojaron cinco tiempos de fases que no habían sido observadas por el analista, pero que mediante el uso del modelo IASP91 corroboramos que se debían efectivamente a fases reales, y no a falsas alarmas. De este interesante resultado, debemos destacar que el ESM y el MAM son los métodos que reconocieron la mayoría de estas fases, ya que ambos detectaron cuatro de estos cinco casos.
Analicemos ahora los resultados desde un punto de vista estadístico a partir de los errores en los picados. Para estimar entonces la exactitud y la precisión de los métodos, se calcularon los errores en el picado de las 25 fases que fueron observadas por el analista tomando como referencia los tiempos que se encuentran en la Tabla 2. En este análisis no se incluyeron las 5 fases (cuyos errores se podrían haber calculado a partir de los tiempos de arribo obtenidos mediante el modelo IASP91) ya que el error sistemático del modelo en la Estación Sismológica de La Plata habría sesgado los resultados. Como es sabido, la media m de los errores es un indicador de la exactitud de los métodos, y la dispersión s es un indicador de su precisión. En la Figura 5 se muestran las medias de los errores y el intervalo (m - s, m + s) que ilustran el comportamiento estadístico de los errores.
Es notable la similitud de las dispersiones de los errores de todos los métodos, ya que para los ocho métodos analizados se encuentran dentro del intervalo . Si bien estos valores son bastante altos, hay que tener en cuenta que estamos trabajando con dos muestras por segundo, y que además estamos considerando absolutamente todos los errores de las fases picadas en los cálculos, algunos de ellos bastante grandes, y que provocan un aumento significativo en la dispersión.

Figura 5. Errores de los diferentes métodos. Los cuadrados negros representan las medias de los errores: los valores positivos son retrasos en el picado y los negativos adelantos. Las barras representan las magnitudes de las respectivas dispersiones.

Respecto de la exactitud de los picados (dada por las medias de los errores), observamos que el EM y el FDM tienden a atrasarse en los picados cerca de 4s (8 muestras), que es un valor demasiado alto. Como se ha propuesto en otros trabajos, esto es fácilmente solucionable incorporando un adelanto adecuado en los algoritmos para los tiempos picados. Sin embargo, para confiar en este procedimiento, resultaría necesario contar con un mayor volumen de datos.
El RAM y el BKM también tienden a atrasarse un poco en sus picados (algo más de 1,5s), que no obstante, es un error por demás aceptable. Por otro lado, el resto de los métodos (ESM, MAM, MBKM, y MCM), presentan una media de errores de picado muy pequeño, en general menor a 1s en valor absoluto. Si tenemos en cuenta que esto equivale solamente a dos muestras, podemos afirmar que estos resultados son muy alentadores. En particular, el ESM y el MAM, además de ser métodos bastante eficientes a la hora de detectar eventos, son métodos bastante exactos para picar las fases.

CONCLUSIONES

Se presentaron y analizaron ocho métodos para realizar picado automático de fases internas de terremotos. Entre ellos, tres son convencionales y ya han sido ampliamente usados y estudiados: el método de Rex Allen (RAM), el método de Baer y Kradolfer (BKM), y el método de Earle y Shearer (ESM). Dos métodos fueron generados para este trabajo a partir de los métodos convencionales: el método de Rex Allen modificado (MAM), y el método de Baer y Kradolfer modificado (MBKM). Y por último, tres métodos que fueron originalmente desarrollados para detección de primeros arribos sísmicos, fueron adaptados en este trabajo para realizar picado de fases de terremotos: el método de Coppens modificado (MCM), el método de la entropía (EM), y el método de la dimensión fractal (FDM). Para analizar y comparar los métodos, se utilizaron once registros de una hora de duración adquiridos en la Estación Sismológica de La Plata (LPA). Estos registros corresponden a eventos de distancias epicentrales superiores a 1000 Kms, magnitud mayor a 6 (Mw), y diversas profundidades de foco.
La aplicación de los métodos es muy simple, y dado que son relativamente sencillos, los tiempos de procesamiento son muy pequeños. La única dificultad para implementar los métodos reside en seleccionar un número determinado de parámetros de forma adecuada, que dependiendo del método, oscila entre dos y cuatro solamente. Los resultados obtenidos con todos los métodos son los esperados en relación a las características de los datos picados, si bien algunos presentan una performance superior a otros.
Los picados falsos obtenidos se debieron en su gran mayoría a falsas fases en la coda de los eventos, y no así a falsas alarmas antes de que llegue la señal o luego de la coda, lo que indica que los métodos presentados son apropiados para realizar futuros estudios de sismicidad. Desde el punto de vista estadístico, si bien todos los métodos exhibieron precisiones muy similares y razonables teniendo en cuenta los datos utilizados, las dispersiones de los picados son algo mayores a lo deseado. Por ello, sería conveniente agregar un último paso a los algoritmos de picado para refinarlos y lograr así resultados más precisos. Este proceso podría incluir, por ejemplo, la búsqueda del punto del inflexión, o del máximo (o mínimo) de la señal más próximo al tiempo arrojado por el método de picado utilizado. Esto por supuesto dependerá del criterio que se adopte para definir el arribo de la fase.
Algunos de los métodos aquí presentados, en especial el método de Earle y Shearer (ESM) y el método de Allen modificado (MAM), detectaron fases internas teóricas que no habían sido observadas por el analista que originalmente realizó la lectura de las bandas. Más allá de las sabidas ventajas que ofrece el uso de métodos efectivos de picado automático en cuanto al ahorro de tiempo y recurso humano, este interesante resultado representa otro argumento a favor de su utilización en la estación sismológica LPA, donde se realiza una exhaustiva lectura de fases secundarias.
El MCM, el EM y el FDM, originalmente desarrollados para detectar un único evento (el primer arribo de los datos sísmicos), fueron correctamente adaptados para su aplicación en problemas de sismología. Sin embargo, esta adaptación resulta compleja para picar de manera exacta y precisa distintas fases y su performance no se destaca por sobre la media del resto de los métodos. En particular, el EM y el FDM tienden a realizar picados de fase demasiado retardados en tiempo.
El método de Allen (RAM) es un método cuyas detecciones son muy confiables puesto que no hemos observado que realice ningún picado falso. Por su parte, el ESM y el MAM resultan ser los más eficientes pues detectan con muy buena exactitud casi todas las fases de interés, y lo hacen generando un número muy pequeño de picados falsos. En base a estos resultados, recomendaríamos su uso en primer término para realizar el picado automático de un dado conjunto de datos. En el mismo sentido, el ESM y el MAM serían dos métodos apropiados para realizar un estudio de sismología global en la Estación Sismológica de La Plata utilizando una gran cantidad de registros.

REFERENCIAS

1. Akaike, H., 1971. Autoregressive model fitting for control. Annals of the Institute of Statistical Mathematics, 23:163-80.         [ Links ]

2. Akaike, H., 1974. Markovian representation of stochastic process and its application to the analysis of autoregressive moving average processes. Annals of the Institute of Statistical Mathematics, 26:363-387.         [ Links ]

3. Aldersons, F., 2004. Toward three-dimensional crustal structure of the Dead Sea region from local earthquake tomography. PhD thesis. Tel Aviv University, Israel.         [ Links ]

4. Allen, R. V., 1978. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68:1521-1532.         [ Links ]

5. Allen, R. V., 1982. Automatic phase pickers: their present use and future prospects. Bulletin of the Seismological Society of America, 72:S225-S242.         [ Links ]

6. Bai, C. Y. and B. L. N. Kennett, 2000. Automatic phase-detection and identification by full use of single three-component broadband seismogram. Bulletin of the Seismological Society of America, 90:187-198.         [ Links ]

7. Baer M., and U. Kradolfer, 1987. An automatic phase picker for local and teleseismic events. Bulletin of the Seismological Society of America, 77:1437-1445.         [ Links ]

8. Boschetti, F., M. Dentith, and R. List, 1996. A fractal-based algorithm for detecting first arrivals on seismic traces. Geophysics, 61:1095-1102.         [ Links ]

9. Buffoni, C., N. C., Sabbione, G. Connon, J. L. Hormaechea, 2009. Localización de hipocentros y determinación de su magnitud en Tierra del Fuego y zonas aledañas. Geoacta 34:75-85. ISSN 0326-7237.         [ Links ]

10. Coppens, F., 1985. First arrivals picking on common-offset trace collections for automatic estimation of static corrections. Geophysical Prospecting, 33:1212-1231.         [ Links ]

11. Denis, A., and F. Crémoux, 2002. Using the entropy of curves to segment a time or spatial series. Mathematical Geology, 34:899-914.         [ Links ]

12. Earle, P. S., and M. Shearer, 1994. Characterization of global seismograms using an automatic-picking algorithm. Bulletin of the Seismological Society of America, 84:366-376.         [ Links ]

13. Jiao, L., and W. M. Moon, 2000. Detection of seismic refraction signals using a variance fractal dimension technique. Geophysics, 65:286-292.         [ Links ]

14. Johnson, C. E., A. Bittenbinder, B. Bogaert, L. Dietz, and W. Kohler, 1995. Earthworm: a flexible approach to seismic network processing. Incorporated Research Institute for Seismology Newsletters, 14(2):1-4.         [ Links ]

15. Kennett, B. L. N. 1991. IASPEI 1991 Seismological tables. Research School of Earth Sciences. Australian National University.         [ Links ]

16. Korvin, G., 1992. Fractal models in the earth sciences. Elsevier.         [ Links ]

17. Küperkoch L., T. Meier, J. Lee, W. Friederich, and EGELADOS Working Group, 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics. Geophysical Journal International, 181:1159-1170.         [ Links ]

18. Leonard, M., 2000. Comparison of manual and automatic onset time picking. Bulletin of the Seismological Society of America, 90:1384-1390.         [ Links ]

19. Leonard M. and B. L. N. Kennett, 1999. Multi-component autoregressive techniques for the analysis of seismograms. Physics of the Earth and Planetary Interiors, 113:247-263.         [ Links ]

20. Luo, Y., M. Marhoon, S. A. Dossary, and M. Alfaraj, 2002. Edge-preserving smoothing and applications. The Leading Edge, 21:136-158.         [ Links ]

21. Mandelbrot, B. B., 1983. The fractal geometry of nature. W. H. Freeman and Company.         [ Links ]

22. Ruud B. O., and E. S. Husebye, 1992. A new three-component detector and automatic single-station bulletin production. Bulletin of the Seismological Society of America, 82:221:237.         [ Links ]

23. Sabbione, J. I., and D. Velis, 2010. Automatic first-breaks picking: New strategies and algorithms. Geophysics, 75:V67-V76.         [ Links ]

24. Sabbione, N., G. Connon, J. L. Hormaechea, M. L. Rosa, 2007. Estudio de sismicidad en la Provincia de Tierra del Fuego, Argentina. Geoacta, 32: 41-50. ISSN 0326-7237.         [ Links ]

25. Saragiotis C. D., L. J. Hadjileontiadis, I. T. Rekanos, and S. M. Panas, 2004. Automatic P-phase picking using maximum kurtosis and κ-statistics criteria. IEEE Geoscience and Remote Sensing Letters, 1:147:151.         [ Links ]

26. Sleeman, R., and T. van Eck, 1998. Robust automatic P-phase picking: an online implementation in the analysis of broadband seismogram recordings. Physics of the Earth and Planetary Interiors, 113:265-275.         [ Links ]

27. Tosi, P., S. Barba, V. D. Rubeis, and F. D. Luccio, 1999. Seismic signal detection by fractal dimension analysis. Bulletin of the Seismological Society of America, 89:970-977.         [ Links ]

28. Turcotte, D. L., 1997. Fractals and chaos in geology and geophysics, 2nd ed. Cambridge University Press.         [ Links ]

29. Whiters M., R. A. Aster, C. Young, J. Beiriger, M. Harris, S.Moore, and J. Trujillo, 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection. Bulletin of the Seismological Society of America, 88:95-106.         [ Links ]

30. Yung, S. K., and L. T. Ikelle, 1997. An example of seismic time picking by third-order bicoherence. Geophysics, 62:1947-1951.         [ Links ]

31. Zhang, H., C. Thurber, and C. Rowe, 2003. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings. Bulletin of the Seismological Society of America, 93:1904-1912.         [ Links ]

Recibido: 15 de septiembre de 2011.
Aceptado: 8 de noviembre de 2011