Geofisica Colombian a pp. 55-60 Bogota, D.C.N°.7 diciembre de 2003 ISSN - 0121 - 2974 Detecci6n en linea de onda-P en sismogramas 3-D utilizando redes neuronales CARLOS ALlATE GERMAN CASTELLANOS Grupo de Control y Procesamiento Digital de Seriales. Universidad Nacional de Colombia Sede Manizales, canalza@epm.net.co RESUMEN Se desarrolla un detector en linea del tiempo de arribo de ondas P sobre registros electronicos de tres componentes de eventos tectonicos. Se lleva a cabo la deteccion en linea de la onda-p empleando redes neuronales (RN), utilizando dos diferentes for- mas: en la primera se sintetiza el detector de LRV, mientras en la segunda se desarrolla un clasificador estadistico con las res- pectivas clases (seiial con onda-p presente y sefial con solo ruido presente). La deteccion se realiza empleando redes neuronales de tipo perceptron multi-capa, en las cuales se aprovechan las siguientes cualidades: su capacidad de mapeo no li- neal de entrada-salida, buena generalizacion, bajo costo computacional, entre otras. Las anteriores cualidades hacen que las RN sean apropiadas para el proceso de sefiales en tiempo real. Las entradas a la RN corresponden a datos normalizados del re- gistro sismico y las caracteristicas que miden el grado de polarizacion de la onda. Los registros procesados se obtuvieron de la base de datos del Instituto Iris y de Ingeominas, y su analisis da como resultado un mejor desernpefio del sistema en la detec- cion del tiempo de arribo, aunque no presenta una exactitud aceptable en la estimacion del tiempo de arribo de la onda P. PALABRAS CLAVE: REDES NEURONALES, PERCEPTRON MULn-CAPA. DETECCION ONDA P ABSTRACT We have developed a P-wave arrival time detector of tectonic events over digital three-component records online. This pro- cess uses two different methods: the first one works with a LVR detector, and the second one uses a statistical classifier of clas- ses (signal with P-wave and signal with seismic noise). Detection uses multi-layer perceptron neural network where is taken advantage of non-lineal input-output, relation, well-generalization, low computational demands, etc. These appropriate quali- ties for real time signal processing. The input data for the system becomes a normalized vector of seismic signal and its polari- zation arguments. Testes were done with data from IRIS and INGEOMINAS database. Although the system has a good perforance, its accuracy is not acceptable. KEYWORDS: NEURAL NETWORK, PERCEPTION, MULllCAPA, DETECTING P-WARE INTRODUCCION La detecci6n de la onda-p es referida a la tarea de deteccion de cambios abruptos, la cual consiste en identificar el instante del tiempo en el que las propiedades del proceso cambian internpesti- vamente, considerando que antes y despues del momento del cambio, las propiedades del proceso SOninvariantes en el tiempo. Tal planteamiento permite que la derivacion formal de algoritmos de detecci6n de cambios abruptos se realice dentro de la c1ase de metodos estadisticos de analisis de procesos aleatorios estacio- narios, para los cuales se consideran dos modelos de procesos Manuscrito recibido para evauaclon el 30 de octubre de 2002. Articulo aceptaco para publicaciOn par el comte Editorial el 30 de octubre de 2003. GEOFislCA COLOMBIANA. 7, DICIEMBRE DE 2003 aleatorios: completamente independientes y los basados en mo- delos regresivos (AR). La prueba de hip6tesis mas empleada corresponde allogaritmo de relacion de verosimilitud (LRV), de la cual se derivan los algoritmos analizados: la suma acumulativa (CUSUM), el LRV generalizado (GLRV), X2 CUSUM YCUSUM bilateral, entre otros. Sin embargo, en caso de tener ruido no gaus- siano a la entrada del detector, la funci6n LRV tiene un caracter no lineal con la entrada y, practicarnente, es muy dificil hallar Unaex- presion explicita con el vector de entrada, 10 que dificulta la sintesis practica del detector. Asi, el detector basado en el LRV para ruido no gaussiano debe considerarse como un dispositivo COnmapeo no lineal. Por cuanto las redes neuronales pueden aproximar Unagran c1asede mapeo no lineal, su aplicacion se ha considerado en la de- tecci6n de sefiales (Lippmann et al., 1989). 55 En el presente articulo se lleva a cabo la deteccion en linea de la onda-p empleando RN, utilizando dos diferentes formas: en la primera se sintetiza el detector de LRY, mientras en la segunda se desarrolla un clasificador estadistico con las respectivas clases (sefial con onda-p presente y sefial con solo ruido presente). La deteccion en linea se realiza deslizando una ventana de N muestras a traves de todo el registro sismico de tres componentes. En cada instante de tiempo se obtienen las funciones de caracterizacion que corresponden a las entradas de la red neuronal, se realiza la propagacion y se decide si en la ventana actual esta el tiempo de arribo de la onda-p. DETECCION EN LINEA DE ONDA-P UTILIZANDO REDES NEURONALES Detector de LR V. La deteccion de sefiales implica decidir si una sefiaI objetivo esta presente 0 no, a partir de datos observados. Existen dos criterios comunes para el disefio e implementacion de un detector de sefiales: el criterio de Bayes y el criterio de Ney- man-Pearson. Ambos criterios dan como resultado el detector de LRY (Luo et aI., 1997), el cual esta compuesto de dos partes: la primera determina la LRY de la funcion densidad de probabilidad (fdp) en presencia de la sefial objetivo y de la fdp en ausencia en la sefial objetivo (ruido). La segunda parte determina la deteccion (salida del detector igual a uno) si la razon de verosimilitud es mayor que un umbra I 0 , Figura 1 La deteccion de un cambio abrupto de la onda-p, se puede sim- plificar al caso de analisis cuando se tiene un solo patron dado s(k), 11m, que representa el cambio abrupto mas ruido blanco gaus- siano (k), 1:-::;k s; n, con media E[(ll)]=O y varianza E[1l2(k)]=cr2EI algoritmo de deteccion puede ser simplificado hasta el correlador que calcule el producto escalar de y con s [4]: Z(n)=(y,s) = !s(k)x(k) > Ao, k=1 Acordes con la estructura del correlador de la figura I, se su- giere una red de tres capas con n neuronas en la capa de entrada, m neuronas en la capa oculta y funcion de transferencia fi (-)y una neurona como capa de salida con funcion de transferencia Sr-), En principio, la estructura del detector secuencial descrito puede ser 56 ALZATE ET AL. implementado con una red neuronal del tipo retropropagacion (figura 2). Los pesos w./l, l:-::;i:-::;ndenotan la conexion entre la neurona i en la capa de entrada y la neuronaj en la capa oculta; mientras yJCk) yhJCk) , 1:-::;j:-::;mcorresponden a su salida y funcion de activacion, respectivamente; WjPJ denota la conexion entre la neurona j de la capa oculta y la neurona en la capa de salida. Se define respectivamente z(k),h(k) y .y como la salida, la funcion de activacion, y el umbral para la neurona en la capa de salida. Sin embargo, en caso de tener ruido no gaussiano a la entrada del de- tector, la funcion de verosimilitud tiene un caracter no lineal con la entrada y, practicamente, es muy dificil hallar una expresion ex- plicita con el vector de entrada, 10que dificulta la sintesis del de- tector. EI detector basado en la relacion de verosimilitud para ruido no gaussiano debe considerarse por tanto, como un dispo- sitivo con mapeo no lineal. Por cuanto las redes neuronales del tipo MLP pueden aproximar una gran clase de mapeo no lineal, su aplicacion se ha considerado en la deteccion de sefiales. En Luo et aI., (1997) se sugiere la realizacion de la funcion de verosimilitud para sefiales de entrada con diferente fdp a la gaussiana, entre- nando adecuadamente los pesos de la RN para obtener el nivel de no linealidad requerido por el detector. Detector basado en la clasificacion simple de eventos: En este caso, se desarrolla un clasificador estadistico con las respectivas clases (sefial con onda-p presente y serial con solo ruido presente). La RN identifica la clase a la cual el patron pertenece debido a su capacidad d~ extraccion de caracteristicas importantes a partir del conjunto de entrenamiento. La clasificacion es una transformacion que mapea los puntos intermedios obtenidos de la extraccion de caracteristicas a puntos correspondientes a cada clase en un es- pacio de decision n-dimensional, donde n es el numero de clases por reconocer. En este esquema los perceptrones multicapa rea- lizan reconocimiento de patrones debido a que poseen un sistema no supervisado para extraccion de caracteristicas (capas ocultas) y un sistema supervisado para clasificacion (capa de salida) .• ENTRENAMIENTO DEL SISTEMA Funciones caracterfsticas Con el objeto de aumentar la efectividad del estimador durante la extraccion de parametres, se emplean las funciones caracteristicas. Entre las mas empleadas estan las siguientes: Razon de energia. En las sefiales sismicas de tres componentes, se puede construir una funcion para detectar movimiento sismico coherente en un plano bidimensional. La razon entre la energia en el plano horizontal y la energia total en el instante de tiempo i esta dado: ER i+N-lI (xl(n)2 +x2(n)2) 1/=; (1)i+-N II(Xl (n)2 +x2(n)2 +XJ (nf) t1=1 GEOFislCA COLOMBIANA, 7, DICIEMBR DE 2003 DETECCI6N EN LiNEA DE ONDA.P EN SISMOGRAMAS 3·D UTILIZANDO REDES NEURONALES donde X" X2 YX3 son las tres componentes de la sefial sismica, N es el tamafio de la ventana. ER tiende a I si la mayoria de la energia sismica esta concentrada en el plano xy, y tiende a 0 si esta concentrada en la componente vertical. Short Term Average/Long Term Average (STA/LTA). Se calcula la razon entre el valor medio de la amplitud absoluta de dos ven- tanas de tiempo consecutivas, una corta (STA) y una larga (LTA), y se compara con un umbral predeterminado que declara la pre- sencia del evento sismico. La determinacion del tamafio optimo de las ventanas depende del contenido frecuencial del sismo (Oonicx, 1998). Se considera una sefial x con N muestras, con la que se construye una funcion caracteristica dada por: para i ~ I, con C(i) como parametro dependiente del tiempo. Se define a 0 (0) =X(0)2, Y en forma recursiva los promedios sobre las ventanas cortas y largas se establecen como: STA(i) = C\ 0 (i) + (1 - CI)STA(i - 1) LTA(i) = C2 0 (i) + (1 - C2)STA(i - 1) para i ~ 1, Y valores iniciales ST A(O) = LTA(O) = 0 (0). Los va- lores de las constantes CI y C2 estan entre 0 y I, valores que se deb en sintonizar en la practica segun la aplicacion, La determi- nacion del cambio es declarada como: Ip = inf {i EN I STA(i) I LTA(i) > C3} siendo C3 el umbral de decision. Funcion rectilineal F. Esta mide el grado de polarizacion lineal de una sefial sismica de tres componentes, la cual aumenta si la energia sismica se concentra en una direccion principal. La onda-p esta generalmente polarizada de forma lineal, de esta forma se espera el maximo valor de F en el tiempo de arribo de la onda-p (An ant et aI., 1995). F=I- A2 A, donde AI es el mayor valor propio y 2 es el segundo mayor valor propio de la matriz de covarianza: [ Var(xl) CoV(XI,X2) Cov(x1 ,XJ)] M= CoV(XI,X2) Var(x2) COV(X2,xJ) COV(XI,XJ) COV(X2,xJ) Var(xJ) donde x], X2 y X3, Y son series de tiempo correspondientes a los tres componentes de una sefial sismica. Var(xj) es la varianza de y COV(XI ,X2) es la covarianza entre XI y X2: GEOFislCA COLOMBIANA, 7, DICIEMBRE DE 2D03 donde 111 Y 112 son los val ores medios de XI y X2 respectiva- mente y N, es el numero de muestras. La funcion recti lineal F es igual a I cuando una onda esta total mente polarizada de forma lineal (arribo de la onda-p) y a 0 cuando no existe polarizacion. Funcion grado de polarizacion DO? Mide el grado de polari- zacion de una sefial de tres componentes: DOP= (A,- A2i + (A2 - AJ)2 + (AJ - A,)2 2(11., + 11.2 + AS (9) (2) donde AI, A2, A3 son los valores propios de la matriz de cova- rianza M. Esta funcion caracteristica puede ser usada para detectar ondas-p y ondas-s debido a que estas poseen alto grado de polari- zacion (Oonincx, 1998). Funcion valor absoluto. Definida de la siguiente forma: Fabs =abs {x (i)},'If i= 1,2, ... , N , (10) (3) (4) don de N es el mimero de muestras y X es una sefial sismica. Esta funcion tiene alta fidelidad y velocidad de proceso (Dai et al., 1997). Funcion envolvente analitica compleja. Usada ampliamente en sistemas de deteccion con el algoritmo de STAlL T A Y esta de- finida como (Kuang et al., 1998): (II) (5) donde xU) es la transformada Hilbert de la sefial sismica xU). Funcion de planaridad. Registra el grado de polarizaci6n de una sefial de tres componentes en un plano, y esta definida por (Oonincx, 1998): P= A, + 11.2 - 2A.J A, + A2 (12) (6) donde 11.\, 11.2, 11.3 son los valores propios de la matriz de cova- rianza M. La onda-s esta polarizada de tal forma que el movi- miento de sus particulas esta confinado a un plano; por tanto se espera el maximo valor de P en el tiempo de arribo de esta onda. (7) Topofogfa de fa RN La funcion valor absoluto tiene alta fidelidad y velocidad de proceso (Dai et. al., 1996), mientras la funcion rectilineal (6) y el grado de polarizacion (DO?) (9) son medidas discriminantes debido a que miden el grado de polarizacion de una serie de tiempo, y la onda-p es unaonda compresional que general mente esta pola- rizada de forma lineal. De esta forma, el espacio de entrada 10 con- forman la funcion valor absoluto, la funcion rectilineal y la funcion DO? Numero de neuronas de entrada. Esta determinado por el mimero de muestras de la ventana correspondiente a la funcion valor absoluto, la funcion rectilineal y la funcion DO? El numero (8) 57 de muestras de la ventana es determinado a traves de la funcion recti lineal y la norma varimax (Anant et. al., 1995): 1.50\-~-t,;-~-;;,.--....-~-,I;o~-,i;,,,,-~,,,,~--,\;;~--.i M (numero de muestras en Ia ...entaroa) Figura 3. Norma varimax promedio. Lx(i)4 V = ---;-,'c-' ~~,- c (~X(i)2) donde xes la serie de tiempo correspondiente a la funcion recti- Ftmcion Reclilineal (N = , 50) Figura 4. Ejernplo de funcion rectilineal. ALZATE IT AL. En las figuras 4 y 5 se muestra la funcion rectilineal y la funcion DOP, respectivamente, de la sefial de tres componentes de la figura 6 para una ventana de 150 muestras con 5 de traslape. EI maximo de las funciones corresponde a la ventana que contiene el tiempo de arribo de la onda-p. Numero de capas ocultas. Cualquier espacio de decision arbi- trario puede ser aproximado por una red neuronal de perceptrones multicapa con una capa oculta de neuronas no lineales (Haykin, 1999). De esta forma es necesaria solo una capa oculta para se- parar el espacio de entrada que corresponde al conjunto de entre- namiento, pero esto no implica que sea el valor optimo en el sentido de generalizacion. Funci6n DOP (N ;< 150) (13) Figura 5. Ejemplo de funcion DOP. Numero de neuronas ocultas. Este mirnero determina la can- tidad de hiperplanos que van a formar el espacio de decision. En general, la cantidad de neuronas ocultas debe ser la minima po- sible para evitar sobreentrenamiento y mala generalizacion. Se realizaron pruebas empezando con dos neuronas ocultas y entre- nando con el algoritmo de Levenberg-Marquardt y el criterio de X 10G corrccoeote n , il ~"w, __ ,___", j .. , 0 500 1000 1500 2000 2500 3000 3500 .000 'f corrccoeree eE~~"-'-j ·4 0 500 1000 1500 2000 2500 3000 3500 '000 If Componente z:.iF----. 1 0 500 1000 1500 2000 2500 3000 3500 .000 Mu.,~s{ll) lineal para ventanas deslizantes de M muestras, M = 2,3,4, ...MMAX' La idea es obtener el valor optirno de M para que la funcion recti- lineal tenga solo un pico prominente que corresponde a la ventana donde esta la onda-p, es decir, minimizar la entropia de la funcion recti linea!. La norma varimax mide la entropia de una serie de tiempo; mientras mas alta es, menor numero de picos tiene la sefial, Se determine que el valor optimo promedio de M es de 150 muestras en cada ventana como se muestra en la figura 3; por tanto, el nurnero de neuronas de entradas es 150 para la funcion valor absoluto, una para la funcion rectilineal y otra para la funcion DO? En total, el numero de neuron as de entrada es 152. Figura 6. Sefial de prueba de tres componentes. 58 GEOFislCA COLOMBIANA, 7, DICIEMBAE DE 2003 DETECCI6N EN LiNEA DE ONDA-P EN SISMOGRAMAS 3-D UTILIZANDO REDES NEURONALES Early Stopping [14]. Los resultados mostraron que el nurnero minimo de neuronas para obtener una buena generalizacion fue de 5 neuronas en la capa oculta. Conjunto de entrenamiento. Inicialmente se utilizaron con- juntos pequefios de entrenamiento correspondientes a 10 datos de onda-p y 10 datos de ruido. EI desernpefio del sistema fue pobre, con una elevada tasa de falsas alarmas y falso rechazo. Posterior- mente, se aumento de forma significativa el nurnero de patrones de ruido. Finalmente, el sistema fue entrenado con 465 patrones de ruido y 13 patrones de onda-p. Algoritmo de entrenamiento. Este corresponde al de Leven- berg-Marquardt, el cual es una modificacion del metodo de mini- mizacion de Gauss-Jordan. Este algoritmo tiene menor tiempo de convergencia que la mayoria de los metodos clasicos como la Regia Delta Generalizada. Se utilizo el criterio de Early Stopping para detener el entrenamiento cuando el error de validacion al- canzo el minimo. Este metodo garantiza que el sistema entrenado tiene buena generalizacion para un conjunto predeterminado de pa- trones de validacion. El error de entrenamiento fue del orden de 10-5• EL esquema del sistema neuronal se ilustra en la figura 7: En el caso del detector neuronal de razon de verosimilitud, la sintesis se realize con los siguientes parametres: Tamaiio de la ventana (N). Debe ser determinado a partir del tipo y del contenido frecuencial de la sefial sismica (Oonincx, 1998). Se determino el valor de N a partir del contenido frecuencial de las sefiales sismicas de tipo regional y local, 10 que resulta en un valor de N=80 muestras, que corresponde a 4s para sefiales mues- treadas a 20 Hz. Numero de neuronas ocultas. Se determine empezando con dos neuronas y probando el desempefio de la red, para un resultado final de 25 neuronas en la capa oculta. En cuanto al conjunto de entrenamiento, este fue determinado de la misma forma que el clasificador antes analizado, empezando Figura 7. Diagrama del detector basado en el clasificador. GEOFislCA COlOMBIANA. 7. OICIEMBRE OE 2003 con conjuntos pequefios y luego incrementando dependiendo del desernpefio del sistema. Al final de las pruebas el conjunto de en- trenamiento estuvo compuesto de 2300 patrones de ruido y 48 pa- trones de onda-p. El algoritmo de entrenamiento fue una modificacion de la Regia Delta Generalizada descrito en Luo et al. (1977). Resultados y conclusiones Los sistemas fueron probados utilizando 77 sefiales sisrnicas de tres componentes muestradas a 20 Hz, que corresponden a sismos regionales y locales de la zona de California y Japon. Adicional- mente se agrego ruido blanco gaussiano a las sefiales para deter- minar el desempefio de los algoritmos frente a condiciones de ruido. La relacion sefial a ruido se calculo hallando la varianza del ruido antes de la onda-p, y hallando la varianza de la senal-ruido despues de la onda-p en una ventana estacionaria: SNR= 1Olog,o(crlr -IJ Los resultados del tiempo de proceso se representan en la tabla 1, los cuales se midieron para los algoritmos en la estacion SUN Enterprise 3500, como tambien en Matlab y en Fortran 95 en un PC sobre Pentium 4. Tabla 1 i i - ! ~-l tiempo prOCCSil Tiempo proceso !, Tiempo proceso Detector (SUN) (PC F95) I (PC Matlab) f-------- ------ -------------"-- -- -------- ClasIficador I 212 ms 950 ms I 4 s LRV I 59 ms I - 400 ms lis En el caso del detector neuronal de razon de verosimilitud, el desempefio dio como resultado un error promedio de 9,1% para una relacion serial a ruido minima promedio de 15dB. Mientras, para el clasificador, el desernpefio del detector resulto en un error promedio de 1,3 % y la relacion sefial a ruido minima promedio fue de 9,67dB. Como se observa el trabajo del clasificador es su- perior al detector basado en la LRY. Se presentan ejemplos de deteccion con el clasificador neu- ronal sobre sefiales sismicas de tres componentes con alta SNR (Figura 8) y baja SNR (Figura 9). El recuadro indica la ventana donde se encuentra el tiempo de arribo de la onda-p. BIBLIOGRAFIA KUANG, C. and MENDEL, J. (1994). First break refraction event picking using fuzzy logic systems, IEEE Transactions on Fuzzy Systems, Vol 2, NO.4. November 1994. 59 LUO, F. and UNBEHAUEN, R. (1997). Applied neural networks for signal processing, University of Erlangen-Nuremberg, 1997. DAr, H. and MACBETH, C. (1997). The application of back-propagation neural network to automatic picking seismic arrivals from single-com- ponent recordings, British Geological Survey, Edinburgh, Scotland. ANANT, K. S. and DOWLA, F. A. (1995). Wavelet transform methods for phase identification in three-component seismograms, University of California, Davis and Lawrence Livermore National Laboratory. 60 ALZATEET AL. OONINCX, P. 1. (1998). Automatic Phase Detection in Seismic Data using the Discrete Wavelet Transform, Centrum voor Wiskunde en Informatica. HA YKIN, S. (1999). Neural Networks: A comprehensivefoundation, Prenti- ceHal!. LIpPMANN, R. and BEKMAN P. 1989. "Adaptive Neural Ney Preprocessing for signal detection in non gaussian noise", Adv. Neural Information Processing System, Vol I, pp. 124-132. GEOFlslCA COLOMBIANA. 7, DICIEMBRE DE 2003