Estimación mejorada de modelos AR multivariados en el análisis de señales EEG

Improved Estimation of Multivariate AR Model for EEG Signal Processing

Jorge Iván Padilla Buriticá^(1)*, Luis David Avendaño Valencia^(2)*, Germán Castellanos Domínguez^(3)*

(1)MSc. Ingeniería-Automatización.jipadilla@unal.edu.co

(2)Ph.D.(c) en Ingeniería línea de automatización. ldavendanov@unal.edu.co

(3)Ph.D. en Ingeniería. Profesor titular, Universidad Nacional de Colombia, sede Manizales. cgcastellanosd@unal.edu.co

(*)Signal Processing and Recognition Group, Universidad Nacional de Colombia, sede Manizales.

Recibido octubre 8 de 2011, aprobado 17 de abril de 2013.

PALABRAS CLAVES

Electroencefalografía, filtro Kalman, modelos autorregresivos multivariados, factor de olvido, recocido simulado.

RESUMEN

Se propone una metodología para mejorar la estimación de las matrices de parámetros de los modelos autorregresivos multivariados, utilizando representación en espacio de estados y el filtro de Kalman, para mejorar la precisión de los parámetros estimados, con un bajo costo computacional. Se consideran dos métodos de adaptación de las matrices de covarianza del filtro de Kalman, para mejorar la velocidad de convergencia conservando la precisión del estimador. Se hacen pruebas sobre datos simulados, sobre una base de datos de electroencefalogramas y también se hace pruebas de la efectividad de la metodología. La contribución está dada en términos de la precisión de los parámetros estimados y el tiempo de estimación, que se reduce hasta en un 40% con un error cuadrático medio de 3%.

KEY WORDS

Electroencephalography, multivariate autoregressive models, Kalman filter, forgetting factor, simulated annealing.

ABSTRACT

This paper proposes a methodology to enhance the estimation of parameter matrices in multivariate autoregressive models. It uses the Kalman filter and state space representation to improve the precision of parameter estimation, while maintaining a reduced computational burden. Two methods of covariance matrix adaptation are considered within the Kalman filter framework to improve the estimator convergence rate while preserving precision. The methodology is tested on simulated data, electroencephalographic recordings and performance. As a result, a reduction of the computational burden of up to 40% is achieved, but reconstruction error reaches as much as 3%.

INTRODUCCIÓN

Los registros de electroencefalograma (EEG), que brindan información espacio-temporal del comportamiento cerebral, son ampliamente utilizados debido a su adquisición no invasiva, al bajo costo para los pacientes, a la alta sensibilidad en el diagnóstico de enfermedades neurológicas (epilepsia, enfermedades neurodegenerativas, tumores cerebrales, entre otras), a la determinación de conexiones fisiológicas relevantes en el cerebro, como afirma Hytti (2006) y al establecimiento de relaciones de causalidad entre las regiones del cerebro, por ejemplo en tareas cognitivas (Gómez-Herrero, 2008). El desarrollo de métodos automatizados de diagnóstico a partir del análisis de EEG implica la solución del problema inverso, que puede realizarse mediante modelos para-métricos o no paramétricos. Aunque se utilizan los últimos modelos, los primeros son preferibles debido a su precisión y capacidad de seguimiento de la dinámica de variables en el tiempo (Poulimenos,2006). En particular, se emplean los modelos autorregresivos multivariado como restricciones espacio-temporales, proporcionando una estructura específica para la distribución de corriente en el cerebro. Así, se han propuesto los modelos autorregresivos multivariados (MVAR), que consisten en una suma lineal ponderada, con la que se describe la señal multivariada, teniendo en cuenta los valores previos de la misma y su relación con los demás canales (Harrison, 2003).

No obstante, en el modelado MVAR se deben tener en cuenta dos problemas: la selección del orden adecuado del modelo y la estimación óptima de sus parámetros (Papas,2006). En cuanto al último problema, la estimación con métodos clásicos (Yule-Walker, Burg) puede generar inconvenientes en la precisión del estimador, sin embargo, el manejo de matrices de alto orden conlleva a un costo computacional excesivo (Haykin, 2001). Para lidiar con esta restricción, los filtros de Kalman se usan debido a que permiten minimizar el funcional de mínimos cuadrados, mediante la estimación recursiva de un modelo en espacio de estados, que toma ventaja en el álgebra matricial para mejorar el procedimiento de estimación de parámetros (Grewal, 2008).

En este artículo, se propone una metodología para mejorar la estimación de las matrices de parámetros de los modelos MVAR, utilizando representación en espacio de estados. Para mejorar la velocidad de convergencia del algoritmo, pero sin alterar la precisión del estimador, se consideran dos métodos de adaptación de las matrices de covarianza del filtro de Kalman: factor de olvido y recocido simulado. La efectividad de la metodología se prueba durante la discriminación de eventos epilépticos sobre datos simulados y en una base de datos de registros EEG. La estructura del artículo es como sigue: en la primera sección, se describen la metodología, los modelos MVAR y la estimación de parámetros mediante el filtro de Kalman. Luego se explica el marco experimental, en el que se describen las bases de datos y las medidas de desempeño empleadas para validar el comportamiento del modelo MVAR. Finalmente, de los resultados obtenidos sobre señales reales y sintéticas, se presentan las conclusiones.

METODOLOGÍA

Modelos Autorregresivos Multivariados (MVAR). Sea una serie de tiempo multivariada esta-cionaria y[k] = [y₁[k] y₂ [k]... y_n [k]]^T, y ∈ ℝⁿ, donde y_j[k] es la j-ésima serie de tiempo, siendo n el número de canales que conforma la señal multivariada de EEG. El modelo MVAR para la serie de tiempo y[k], se define como [8]:

donde Α_i ∈ ℝ^{n x n}, i = 1, ... , p es la i-ésima matriz de parámetros del modelo MVAR, que se define como, Α_i ={a_uv⁽ⁱ⁾ : u, v = 1, ... , n} mientras η[k] es un vector aleatorio no correlacionado con media cero y matriz de covarianza C_n ∈ ℝ^{n x n}, cuya distribución se asume Gaussiana. En el modelo MVAR, se puede definir un vector de parámetros x, tal que

tal que el modelo MVAR descrito en la Ec. (1) puede ser reescrito en forma alternativa como:

donde W [k] es una matriz de tamaño n² x pn², en la que se encuentran los valores previos de la serie de tiempo y[k], siendo [W₁[k] = ... W_p [k]] , con elementos:

La Ec. (2) considera el modelo MVAR como una relación lineal simple entre el vector de parámetros x y la serie de tiempo multivariada y[k], dada a través de la matriz de medición W [k] que contiene la información de los valores pasados de y [k]. Existen dos tareas usualmente relacionadas con el modelo MVAR, definido en las Ec. (1) ó (2): el primero, denominado problema directo, en el cual se conoce tanto la realización del ruido del proceso η[k], el valor inicial de la serie de tiempo y [0], como las matrices parámetros A₁, ... , A_p y se desea obtener una realización del proceso aleatorio multivariado y [k]; y el segundo, denominado problema inverso, que incluye la realización del proceso aleatorio multivariado y [k], cuando se desea conocer el orden del modelo MVAR(p), las matrices de parámetros A₁, ... , A_p y, eventualmente, las propiedades del ruido del proceso η[k].

El problema de estimación directo puede resolverse aplicando de forma iterativa cualquiera de las Ecs. (1) ó (2) bajo la condición inicial dada y [0]. Mientras, el problema inverso se resuelve utilizando métodos de estimación basados en mínimos cuadrados, como el descrito por Neumaier (2001). No obstante, la minimización directa de la función de costo de mínimos cuadrados generalmente hace necesaria la operación sobre matrices de tamaño muy elevado, que en el caso particular de registros EEG multivariados, puede hacer imposible la estimación de los parámetros del modelo. Una aproximación alternativa consiste en la estimación iterativa del vector de parámetros x, con métodos como el de mínimos cuadrados recursivos o el filtro de Kalman (Milde,2010). En este caso, el modelado debe considerar que el vector de parámetros evolucione de forma restringida a lo largo del tiempo hasta que converja en un tiempo finito. Así, se asume que la tasa de cambio del vector de parámetros entre dos instantes de tiempo es inferior a un valor dado σ_v, de acuerdo con ∆ x [k] = v[k], donde v [k] es un vector aleatorio de ruido blanco Gaussiano no correlacionado con media cero y matriz de covarianza C _v , definida como C _v = σ _v I . En particular, el operador de diferencia ∆, se puede aproximar en el caso más simple (conocido como de suavidad estocástica o smoothnesspriors (SP))en la forma [10]:

Las ecuaciones (2) y (3) definen una representación en espacio de estados, a partir de la cual es posible realizar la estimación del vector de parámetros x utilizando el filtro de Kalman.

Estimación de Parámetroscon el filtro de Kalman (KF): En la representación mediante espacio de estados dada por las Ecs. (2) y (3) se pueden hacer ciertas suposiciones sobre la forma de cambio temporal de la matriz de covarianza del filtro de Kalman C _v, a fin de mejorar la velocidad de convergencia. De hecho, en términos generales, se puede decir que un valor grande de esta covarianza hace que los datos antiguos sean descartados. En particular, se consideran dos aproximaciones descritas por Haykin (2001): factor de olvido (forgetting factor - FF) y recocido simulado (simulatedannealing - SA). En el método FF, la matriz de covarianza C _v se define como

tal que la actualización de la covarianza P [k] = λ -1 P [k - 1] en el filtro de Kalmanes P^-[k], donde λ, 0 « λ ‹ 1 o factor de olvido del estimador que define la velocidad con que se descartan los valores antiguos en el estimador. Una vez elegido C _v [k] de esta forma, se genera una ventana exponencial en el tiempo que sopesa los datos estimados, dándole mayor ponderación a las estimaciones de estado x^ [k] más recientes, haciendo de esta forma que los datos antiguos tengan menor influencia en los parámetros. La introducción de este factor de olvido flexibiliza el filtro, haciéndolo que la estimación sea más sensible a los datos nuevos.

En forma alternativa, el método SA define la matriz de covarianza C _v = αI , donde α ∈ ℝ ⁺ es un escalar positivo cuyo valor evidencia un decremento a medida que continúa el en-trenamiento. De esta forma, el valor de la covarianza dis-minuye a lo largo del tiempo, haciendo que a medida que avanza la estimación se va disminuyendo la variabilidad del vector de parámetros.

Estimación reducida de parámetros mediante el filtro de Kalman: Con el fin de reducir el número de parámetros para estimar en el caso de las señales de EEG, se pueden hacer las siguientes consideraciones acerca de la estructura de las matrices de parámetros del modelo MVAR. Por cuanto, los registros son adquiridos mediante el estándar EEG 10-20 (cuando se garantiza simetría en la ubicación de los electrodos en el cuero cabelludo del paciente) se tiene n = 19.

Como se detalla en la sección Marco experimental, se halla el orden del modelo MVAR (4), para el cual se hace la reducción de 19² * 4 = 1444 parámetros a 19 * 20 * 2 = 760 parámetros. La segunda consideración consiste en suponer que la interacción entre dos canales lejanos es nula. De esta forma, se puede suponer que los canales sólo se relacionan con los vecinos más cercanos, haciendo que la matriz de parámetros sea dispersa. La cantidad de parámetros para estimar dependerá de la forma de la matriz de parámetros y de la disposición espacial de cada canal. En ambos casos, se debe considerar que hay que modificar la estructura de la matriz W [k].

MARCO EXPERIMENTAL

Base de datos: La base de datos empleada pertenece al Ins-tituto de Epilepsia y Parkinson del Eje Cafetero y está conformada por 340 registros de 20 canales, correspondientes al estándar EEG 10-20. La base de datos consta de 120 registros de epilepsia generalizada, 120 registros de epilepsia focalizada y 100 registros normales, los cuales fueron adquiridos con una frecuencia de muestreo de 256 Hz, con resolución de 12-bit y 2 minutos de duración.

Los registros cuentan con el diagnóstico de epilepsia confirmado el análisis clínico de neurólogos. De igual manera, los modelos MVAR se prueban sobre un conjunto de datos consistente de registros sintéticos, que consisten en una señal multivariada de dos canales. La señal es generada por un modelo MVAR con matrices de parámetros conocidas y con diferentes realizaciones de η[k], donde la media es cero y la matriz de covarianza es conocida. Las matrices de parámetros y la matriz de covarianza del ruido son

Ambas matrices de parámetros definen un modelo MVAR con dos pares de polos complejos conjugados definidos en π/3 rad y π/6 rad, con magnitudes de 0.9 y 0.94 respectivamente.

Ajuste de los modelos MVAR. El ajuste de los modelos MVAR consiste en determinar el orden del modelo p así como determinar los valores de las matrices de parámetros A₁, ..., A_p Con el fin de determinar el orden del modelo, se puede utilizar el criterio de información Bayesiano (BIC), que se define como: BIC( p ) = -2 ln |C_n| + 2 p ln N , donde |C_n| denota el determinante de la matriz C_n, N es la cantidad de muestras en el modelo y p el orden. El orden óptimo se selecciona como aquel que maximiza el criterio BIC. La calidad del ajuste de los modelos MVAR se puede medir mediante dos criterios, el error cuadrático medio (MSE) y el índice de correlación. Ambas medidas se definen como:

donde y_i, es el i-ésimo canal de la señal EEG, ŷ_i es el i-ésimo canal estimado de la señal EEG, σ_yŷ² es la covarianza cruzada entre y_i eŷ_i; σ_y yσ_ŷ son las desviaciones estándar de y_i eŷ_i respectivamente.

RESULTADOS

Señales sintéticas: Se generan señales sintéticas que garantizan una estructura aleatoria de un sistema estable, en concordancia con el apartado anterior para diferentes valores iniciales y realizaciones de ruido. Mediante el algoritmo propuesto se estiman los parámetros del modelo MVAR (2) y se calculan las medidas de desempeño. El mismo procedimiento se repite 50 veces, y luego se calculan los valores medios. El resultado obtenido de MSE es de 0.0669 y es de -6.7254.

La Fig.1 muestra el patrón de convergencia del primer parámetro de una señal simulada con respecto al número de iteraciones del algoritmo. Se puede observar que al hacer la adaptación de la covarianza con FF, la estimación de los pa-rámetros converge más rápido que cuando se tienen los otros métodos utilizados.

Señales Reales. El primer paso en las pruebas de estimación de parámetros del modelo MVAR para señales EEG, consiste en calcular el orden óptimo del modelo, para lo cual se uti-lizó el criterio BIC, variando p desde 1 hasta 8. El resultado de esta prueba dice que un valor dep=3 es adecuado para estas señales, como se ve en la Figura 2.Este es un orden bajo comparado con los modelos univariados AR, que necesitan de un valor más alto para capturar de forma precisa la dinámica de las señales EEG (Chen, 1996)

Estimación completa de parámetros de señales reales: En la Fig. 3 se muestran las matrices 20x20 de parámetros del modelo que corresponden al número de electrodos del EEG. En la figura, para la interpretación de la escala de grises, los colores claros representan los electrodos que tienen menor interacción en el registro EEG. Se muestran los resultados de la prueba realizada sobre una señal de epilepsia generalizada de la base de datos de señales reales. En esta prueba los parámetros estimados con KF para el modelo MVAR(3), se agrupan en conjuntos de 4x4. El patrón de comportamiento se mantiene para todas las señales de la base de datos y para todos los estimadores mencionados en la metodología. El patrón de comportamiento encontrado es considerado para hacer la estimación reducida de parámetros.

En las notaciones “Focalizada 1" o “Generalizada 1" los índices numéricos corresponden al orden de numeración de cada registro en las bases de datos de epilepsia focalizada y generalizada, respectivamente, y sobre los cuales se hicieron las estimaciones. En la Tabla 1 puede observarse que KF muestra el MSE más bajo y mayor índice de correlación entre la señal real y la reconstruida, con un valor muy cercano a 1. También puede observarse que la estimación de parámetros hecha con las diferentes formas de adaptación de covarianza permite llegar a resultados de estimación cercanos a los obtenidos mediante la metodología simple. En la Figura 4, El error indicado corresponde a la diferencia absoluta de reconstrucción, que se calcula mediante la diferencia entre la señal original y su reconstrucción.

Tabla 1. Desempeño para la estimación de parámetros con estimación completa y con estimación reducida

La reconstrucción de la señal EEG de epilepsia focalizada 1, mostrada en la Fig. 4, obtenida con los parámetros del modelo MVAR (3), conserva las características temporales de la señal original. Este resultado es coherente con los resultados presentados en la Tabla 1 para la misma señal, en la cual el peor resultado se obtiene con SA.

Otro aspecto considerado, es el costo computacional en la estimación de estos parámetros. Se necesitan en promedio 3900 segundos para KF, un 40% menos para FF y un 4% más para SA, como se muestra en la Tabla 2.

Estimación reducida de parámetros de señales reales. Teniendo en cuenta el patrón de comportamiento que presentan los parámetros obtenidos con la estimación completa se reduce el número de parámetros para estimar. Se plantea una forma de estimación en la cual solo se tienen en cuenta los parámetros contenidos en grupos de . De esta forma, el costo computacional del algoritmo es reducido en un 96% para el caso de KF con señales generalizadas como puede verse en las Tablas 2 y 3 respectivamente.

La reconstrucción de las señales con el modelo MVAR (3), presenta un error pequeño como se ve en la Fig 6. Además, existe mayor relación e interacción entre los elementos de las diagonales de las matrices de parámetros; estos parámetros corresponden a la relación de cada electrodo consigo mismo como se muestra en la Fig 5.

Realmente, el algoritmo convencional de Kalman requiere el cálculo de la matriz inversa, cuya complejidad computacional se incrementa significativamente con el aumento en la dimensión de la matriz/vector espacio estado, con lo cual su complejidad computacional es alrededor de O(n²). Para manejar este problema, se introduce una modificación de los modelos canónicos de espacio estado, particularmente, se lleva a cabo la estimación mediante la estimación de parámetros con la adaptación de covarianza de factor de olvido en el Filtro de Kalman. Como resultado para la variante de For-getting Factor, la dimensión de n se reduce hasta en un 60%, así que una primera estimación de la complejidad obtenida es alrededor de O (0.6n²)

CONCLUSIONES

En este artículo se ha propuesto una metodología para mejorar la estimación de las matrices de parámetros de los modelos MVAR, utilizando representación en espacio de estados.

Esta representación permitió organizar las matrices de parámetros de forma tal que se pudieran aprovechar propiedades de simetría y dispersión de las matrices, para estimar un conjunto reducido de parámetros. Bajo estas condiciones es posible utilizar el filtro de Kalman para reducir el error de reconstrucción y el costo computacional en la estimación.

En general, los resultados de precisión del algoritmo propuesto son similares a los obtenidos en el esquema convencional del filtro Kalman, como se observa en la Tabla 1.Sin embargo, de los resultados obtenidos en la Tabla 2 (tiempo promedio de estimación) y de la Fig 6. (la reconstrucción obtenida de las señales reales de EEG), se puede apreciar que se obtiene una reducción en el tiempo de estimación mediante la estimación de parámetros con la adaptación de covarianza de factor de olvido en el Filtro de Kalman.

Así, se logró reducir el costo computacional en la estimación de parámetros del modelo, pasando de estimar 1200 parámetros por cada iteración a estimar tan solo 576 del modelo MVAR (3). Adicionalmente se consideraron dos estrategias para mejorar la velocidad de convergencia del filtro de Kalman, utilizando factor de olvido y recocido simulado para la adaptación de las matrices de covarianza. Las pruebas hechas mostraron que el estimador tiene menor error cuadrático medio cuando no hay adaptación de las covarianzas, pero al adicionar un factor de olvido, el estimador converge de forma más rápida, manteniendo el error de reconstrucción. Se encontró que el Filtro de Kalman tiene un mejor desempeño en la estimación de parámetros del modelo MVAR tanto en las señales sintéticas como en las señales reales, sin embargo, tiene más costo computacional que las estrategias de adaptación de covarianza, con las que se logran resultados similares. Como trabajo futuro se pretende adecuar el modelo para reducir aún más el costo computacional, sin que se vea afectada la calidad de la estimación de los parámetros, con el fin de aplicar estos modelos para la ubicación de focos epilépticos.

AGRADECIMIENTOS

Los autores agradecen al Instituto de Epilepsia y Parkinson del Eje Cafetero, Pereira por su ayuda en el acceso a la base de datos. Este trabajo se realiza bajo el apoyo de los proyec-tos: Desarrollo de un sistema automático de mapeo cerebral y monitoreo intraoperatorio cortical y profundo: Aplicación a la neurocirugí a - 111045426008 COLCIENCIAS y Análisis dinámico de relevancia en bioseñales, 20201007075, Universidad Nacional de Colombia.

REFERENCIAS

Chen, C., & Davis, R. (1996). Order Determination for Multivariate Autoregressive Processes using Resampling Methods. Journal of multivariate analysis, 57(2), 175-190.         [ Links ]

Gómez-Herrero G., Atienza, M. Egiazarian, K., & Cantero, J. L. (Nov, 2008) Measuring directional coupling between EEG sources. Neuroimage, 3(43), 497-508.         [ Links ]

Grewal, M. & Andrews, A. (2008) . Kalman Filtering: theory and practice using MATLAB. (3ra ed). New York, Wiley.         [ Links ]

Harrison, L., Penny, W. & Friston K. (Aug, 2003). Multivariate autoregressive modeling of fMRI time series. Neuroimage, 19 (4), 1477-1491.         [ Links ]

Haykin, S. Kalman filtering and neural networks. (2001). New York: Wiley Interscience.         [ Links ]

Hytti, H., Takalo, R., & Ihalainen, H. (Apr, 2006). Tutorial on multivariate autoregressive modeling, Journal of Clinical Monitoring and Computing, 2 (20), 101-108.         [ Links ]

Kitagawa, G., & Gersch., W. (1996). Smoothness priors analysis of time series. Lecture notes in statistics. Berlin, Springer.         [ Links ]

Milde, T., Leistritz, L., Astolfi, L., Miltner, W., Weiss, T., Babiloni, F., & Witte, H., (2010). A new Kalman filter approach for the estimation of high-dimensional time-variant multivariate AR models and its application in analysis of laserevoked brain potentials, Neuroimage, 50(3), 960-969.         [ Links ]

Neumaier, A. &Schneider, T. (2001). Estimation of parameters and eigenmodes ofmultivariate autoregressive models, ACM Transactions on Mathematical software, 27(1), 25-57.         [ Links ]

Pappas, S., Leros, K. & Katsikas, K.. (2006) Joint order and parameter estimation of multivariate autoregressive models using multi-model partitioning theory, Digital Signal processing, 16 (6),782-795.         [ Links ]

Poulimenos, A. G,, & Fassois, S. D. (May, 2006) Parametric Time-domain Methods for Non-stationary Random Vibration Modelling and Analysis: A Critical Survey and Comparison. Mech. Syst. Signal Pr.20 (4), 763-816.         [ Links ]