CC BY
19
ОПТИМИЗАЦИЯ НАБОРА БАЗИСНЫХ ФУНКЦИЙ ПРИ РЕКОНСТРУКЦИИ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПО КОРОТКИМ ВРЕМЕННЫМ РЯДАМ ЭЭГ ВО ВРЕМЯ ЭПИЛЕПТИЧЕСКОГО ПРИПАДКА
А.Н. Северюхина
Решается задача построения компактной математической модели по коротким записям электроэнцефалограмм больных абсанс эпилепсией во время припадка. Полученная эмпирическая модель может быть полезна при решении целого ряда задач, включая сегментацию временных рядов, кластеризацию полученных фрагментов и пр. Для решения поставленной задачи предлагается использовать методы оптимизации реконструируемых многокомпонентных моделей. Показано, что использование методов оптимизации позволяет получить адекватную математическую модель, значительно сократив при этом количество базисных функций, входящих в состав ее уравнения.
Ключевые слова: Модельные уравнения, экспериментальные данные, временные ряды, оптимизация набора базисных функций, электроэнцефалограмма, абсанс эпилепсия.
Введение
Одновременно с усовершенствованием вычислительной техники активно развивались подходы к созданию эмпирических математических моделей на основе экспериментальных данных. Такие модели нашли свое применение во многих отраслях науки и прошли свой путь от простейших функций вида х = / (Ь), аппроксимирующих множества экспериментальных точек на плоскости (х,Ь) [1], к сложным многомерным системам нелинейных разностных и дифференциальных уравнений.
При реконструкции модельных уравнений типичной ситуацией является отсутствие априорной информации о структуре функций, входящих в уравнение. Поэтому для их аппроксимации используются различные универсальные формы. Например, в соответствии с теоремой Вейерштрасса любую непрерывную функцию можно аппроксимировать полиномом [2]. Однако на практике подобные модели зачастую оказываются неработоспособными - в подобном громоздком базисе почти всегда присутствуют лишние элементы. Строго говоря, лишней следует считать такую базисную функцию, коэффициент перед которой отличен от нуля только из-за присутствия шумов. Присутствие таких базисных функций немного уменьшает
ошибку аппроксимации моделью тренировочного ряда, однако приводит к большим ошибкам во всех остальных областях восстановленного фазового пространства. Одним из возможных подходов в решении данной проблемы - исключение из базиса лишних функций. Такой подход может существенно расширить область фазового пространства, где построенная модель будет хорошо описывать поведение объекта.
Объектом исследования являются электроэнцефалограммы детей больных абсанс эпилепсией. Экспериментальные данные предоставлены нам коллегами из Института высшей нервной деятельности и нейрофизиологии РАН (Москва). Данные представляют собой набор временных реализаций девятнадцатиканальной поверхностной ЭЭГ. Временные ряды записаны с частотой 256 Гц; каждая запись содержит по одному эпилептическому припадку, один из которых представлен на рис. 1. В работе использовались временные реализации 7 пациентов, всего 26 припадков.
Следует отметить, что временные ряды ЭЭГ характеризуются высокой степенью нестационарности [3,4]. Изменения вида ЭЭГ происходят почти постоянно и могут быть вызваны переходами между различными физиологическими и эмоциональными состояниями [4-6]. Особенно заметны изменения во время эпилептического припадка [7]. Период времени, в течение которого статистические (в частности спектральные) свойства ЭЭГ можно считать неизменными, по разным оценкам составляет от 0.2 секунды до одной минуты [8,9]. По этой причине для анализа таких данных очень важно уметь строить компактные математические модели.
Целью работы является построение компактной эмпирической модели сигнала электроэнцефалограммы во время эпилептического припадка при абсанс эпилепсии.
Рис. 1. а - временная реализация ЭЭГ-сигнала; б - увеличенный фрагмент, содержащий эпилептический припадок
Построенную модель планируется использовать для сегментации ЭЭГ-сигнала на квазистационарные фрагменты с их последующей классификацией.
За последние несколько лет было опубликовано большое число работ, посвященных исследованию эпилептической активности человека и животных с помощью различных походов нелинейной динамики, включая методы авторегрессионного анализа [10-13], а также применение непрерывного вейвлет-преобразования [14-17] и пр.
В данной работе по экспериментальным данным реконструировалось модельное отображение, аппроксимирующая функция которого имела вид степенного полинома. В первую очередь, решалась задача подбора оптимального набора параметров модельного уравнения. Однако использование критерия насыщения ошибки аппроксимации, а также критерия минимума целевой функции оказалось неэффективным. По этой причине было решено использовать другой подход для нахождения оптимального вида модельного отображения - оптимизацию набора базисных функций.
1. Методика исследования
В данной работе модельное уравнение строилось в виде разностного уравнения
Xi+DP = G(xi,xi-Iag,...,xi-(D-i)iag). (1)
Здесь (xi, xi-iag,...,xi-(D-i)iag) - вектор состояния, lag - расстояние между точками вектора состояния (лаг), DP - дальность прогноза, D - размерность отображения. G - аппроксимирующая функция, представляющая собой линейную комбинацию базисных функций fk
G = ак fk(xi,Xi-iag,...,xi-(D-i)iag), (2)
k=l
где функции fk (Xi, Xi-lag,...,Xi—(D—l)lag) имеют вид
DD
fk = П xPik, где £ Pik < S. (3)
Здесь в - порядок полинома, или степень нелинейности. Коэффициенты подбирали методом наименьших квадратов.
Как уже отмечалось выше, одной из основных проблем при реконструкции уравнений является чрезмерная громоздкость аппроксимирующих функций [18]. Так, для модели вида (1) характерно быстрое увеличение числа слагаемых, (соответственно, коэффициентов и базисных функций) согласно формуле
« (в + О)! СОе// = ^Ш! • (4)
К примеру, в модельное уравнение с параметрами О = 2 и в = 2 будут входить 6 слагаемых
С = а1 + а2х\ + а3х1 + а4х\х2 + а5х2 + абх^
Базисные функции, входящие в состав данного модельного уравнения,
/l = 1, /2 = XI, /3 = x\, /4 = Х1Х2, /5 = X2, /e = x\.
При D = 5 и s = 5 количество слагаемых, входящих в состав модельного отображения, будет уже 252.
Таким образом, подбор оптимальных параметров (размерности D, порядка s, расстояние между точками lag и др.) является важным этапом при реконструкции модельного уравнения.
Существует ряд методов автоматического выбора параметров многочлена, которые опирается на идею минимизации целевой функции вида [18]
Ф(сое// ) = gi(e2)+ g2(coe//), (5)
где g1, g2 - возрастающие функции своих аргументов. Среднеквадратичная ошибка
аппроксимации определяется по формуле
1 2
= iXi+DP ~ G(Xi,Xi-lag, ..., Xi_{D_ 1)lag)) . (6)
i=1
Первое слагаемое определяет вклад ошибки аппроксимации (6), а второе - размера модели (4). Наличие минимума целевой функции (5) следует ожидать зачастую при промежуточных значениях размера модели, так как при малом размере многочлена слишком велика эмпирическая ошибка, а при большом - слишком велико становится второе слагаемое. Таким образом, целевая функция обеспечивает некоторый компромисс между ошибкой аппроксимации и размером модели (количеством коэффициентов сое//).
В качестве такой целевой функции в данной работе использовался критерий Шварца [19]
N'^ _2 I ln(N')_
Sch(coeff) = J lneJ + -1—Lcoeff, (7)
где эффективная длина ряда
N' = N - DP - (D - 1)lag. (8)
Для поиска оптимальных параметров модельного уравнения (1) был осуществлен простой перебор порядков s € {1, 5} и размерностей D € {1, 5}. Однако функция Шварца (7) либо вообще не имела минимума, либо минимум достигался при высоких размерностях D модельного отображения, когда количество коэффициентов начинало приближаться к количеству точек в анализируемом временном ряде. Однако в таком случае неизбежно возникает проблема переобучения модели [20].
Возможный подход к решению данной проблемы - процедура оптимизации набора базисных функций.
Существует два распространенных подхода к оптимизации [21]: на основе удаления базисных функций из полного полиномиального базиса (Backward Elimination) и на основе добавления к изначально пустому базису (Forward Selection). В данной работе были протестированы оба метода, при этом оба подхода давали одинаковый результат. В дальнейшем для решения задачи оптимизации использовали метод Backward Elimination. В качестве критерия для нахождения оптимального базиса использовали минимум функции Шварца (7).
Для оценки адекватности полученной после оптимизации модели проводился анализ регрессионных остатков на коррелированность [18]
— Хг+ДР — О(Хг,Хг-1ад,...,Хг-(Д-1)гад).
(9)
На основе данного критерия было установлено, что структуру сигнала во время эпилептического припадка наилучшим образом описывают многомерные модели (когда в вектор состояния входит порядка 100 элементов).
К таким многомерным моделям и применяли методы оптимизации набора базисных функций. Рис. 2 иллюстрирует процесс оптимизации для линейной модели с размерностью О — 100, то есть в структуру модели входит 100 базисных функций
1,Х1,Х2,Хз, ...,Х99,Хюо,
а само модельное уравнение выглядит следующим образом:
О — ао + а1Х1 + а2Х2 + азХз + ... + а99Х99 + аюоХюо-
Напомним, что Хк, к € {1,100}, - точки, взятые из анализируемого временного ряда и формирующие вектор состояния. Таким образом, для линейной модели вектор состояния и набор базисных функций фактически полностью совпадают, за исключением /1 — 1.
На рис. 2, а приведена диаграмма, которая описывает процесс оптимизации: по оси ординат отмечены шаги оптимизации (от 1 до 100), по оси абсцисс отмечены номера базисных функций (также от 1 до 100). На каждом шаге оптимизации из набора удаляется одна базисная функция. Таким образом, каждый столбец показывает, на каком шаге оптимизации из набора удалена конкретная базисная функция (например, базисная функция Х35 была удалена на 80-м шаге процесса оптимизации).
30 40 50 60 Сдвиг по времени
Рис. 2. а - диаграмма, описывающая процесс оптимизации; б - зависимость величины функции Шварца от шага оптимизации; в - автокорреляционная функция, соответствующая участку временного ряда ЭЭГ, на котором осуществляется процесс оптимизации
Сплошная линия на рисунке соответствует номеру шага, на котором был достигнут минимум функции Шварца, изображенной на рис. 2, б. Столбцы, оказавшиеся выше этой горизонтальной линии, соответствуют номерам базисных функций, которые остались в результате процесса оптимизации и составляют оптимальный набор.
Анализ набора оставшихся по итогам оптимизации базисных функций и графика автокорреляционной функции (рис. 2, в) показал, что лучшая аппроксимация соответствует модели, вектор состояния которой составляют точки временного ряда, отстоящие друг от друга на неравные интервалы (так называемые неравномерные вложения). При этом большая часть точек, по которым делаем прогноз, лежит рядом с предсказываемой точкой и еще несколько точек лежат через характерный период Т. Так, например, для процесса оптимизации, представленного на рис. 2, оптимизированный набор будет выглядеть следующим образом:
1, Х2 ,Хз ,Х4, Х5, Х7, Х10, Х13, Х15, Х18, Х21 ,%23, Х27, Х47, Х77, Х95.
Рис. 3. Остатки (а) и автокорреляционная функция остатков (б) для оптимизированной модели
На рис. 3, б приведен график автокорреляционной функции остатков для оптимизированной модели. Значения автокорреляционной функции практически не выходят за границы 95-% квантиля для случая независимых нормально распределенных случайных величин [22]. Такая модель хорошо описывает корреляции между точками. Важно, что данный результат был достигнут для всех использованных реализаций.
Отметим, что похожий результат был получен для крысиных ЭЭГ в работе [10], где вектор состояния формировался из координат, взятых с задержкой 1ад1, и одной дополнительной точки, взятой через характерный период (1ад2 = Т).
Заключение
Известно, что временные ряды ЭЭГ характеризуются высокой степенью нестационарности [3], поэтому для таких данных очень важно иметь возможность построения компактных математических моделей. Такие модели могут быть применены не только для анализа участков ЭЭГ во время эпилептического приступа, но и, в принципе, для любой ЭЭГ, например, для разделения ЭЭГ на квазистационарные сегменты и дальнейшей их кластеризации с целью нахождения различных режимов.
В рамках данной работы показано, что применение методов оптимизации позволяет значительно уменьшить размер модельного уравнения. В итоге получается компактная модель (примерно 10-15 слагаемых) с неравномерными вложениями, включая точки, взятые в районе характерного периода Т. Полученная модель хорошо описывает динамику электроэнцефалограммы во время эпилептического припадка.
Приношу благодарность профессору, доктору физико-математических наук Без-ручко Борису Петровичу за постоянное внимание и руководство работой, а также кандидату физико-математических наук Диканеву Тарасу Викторовичу за помощь, которая была оказана при проведении исследования, результаты которого изложены в статье.
Библиографический список
1. Демиденко Е.З. Линейная и нелинейная регрессия. М.: Финансы и статистика, 1981. 302 с.
2. Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. М.: Физматлит, 2001. Т. 3. С. 582.
3. Deistler M., Prohaska O., Reschenhofer E., Vollmer R. Procedure for different stages of EEG background activity and its application to the detection of drug effects // Electroenceph. clin. Neurophysiol. 1986. Vol. 64, № 4. P. 294.
4. Palus M. Nonlinearity in normal human EEG: Cycles, temporal asymmetry, non-stationarity and randomness, not chaos // Biol. Cybern. 1996. Vol. 75. P. 389.
5. Каплан А.Я. Нестационарность ЭЭГ: Методологический и экспериментальный анализ // Успехи физиологических наук. 1998. Т. 29, № 3. C. 35.
6. Gribkov D., Gribkova V. Learning dynamics from nonstationary time series: Analysis of electroencephalograms // Phys. Rev. E. 2000. Vol. 63, № 6. P. 6538.
7. Jefferys J.G.R. Basic mechanisms of focal epilepsies // Exp. Neurol. 1990. Vol. 75. P. 127.
8. Каплан А. Я. Проблемы сегментного описания ЭЭГ человека // Физиология человека. 1999.Vol. 25, № 1. P. 125.
9. Barlow J.S. Methods of analysis of nonstationary EEGs, with emphasis on segmentation techniques: A comparative review//J. Clin. Neurophysiol. 1985. Vol. 2. P. 267.
10. Сысоева М.В., Сысоев И.В. Математическое моделирование динамики электроэнцефалограммы во время эпилептического припадка // Письма в ЖТФ. 2012. Т. 23, вып. 3.
11. Филина Е.В. Динамика локальных потенциалов мозга при абсанс-эпилепсии: Эмпирическое моделирование // Изв. вузов. Прикладная нелинейная динамика. 2011. Т. 19, №4. С. 109.
12. Сысоева М.В., Диканев Т.В., Сысоев И.В., Безручко Б.П. Анализ связей между отведениями электроэнцефалограмм крыс до и во время эпилептического припадка с помощью предсказательных моделей // Вестник Нижегородского университета им. Н.И. Лобачевского. 2013. № 1-1. С. 73.
13. Чепурнов С.А., Аристов А.В., Бердиев Р.К., Аббасова К.Р., Чепурнова Н.Е. Перспектива современных подходов к анализу ЭЭГ при эпилепсии с целью поиска предикторов приступов // Успехи физиологических наук. 2010. Т. 41, № 4. С. 27.
14. Ситникова Е.Ю., Короновский А.А., Храмов А.Е. Анализ электрической активности головного мозга при абсанс-эпилепсии: Прикладные аспекты нелинейной динамики // Изв. вузов. Прикладная нелинейная динамика. 2011. Т. 19, № 6. С. 173.
15. Ситникова Е.Ю., Грубов В.В., Храмов А.Е., Короновский АЛ. Возрастные изменения частотно-временгой структуры сонных веретен на ЭЭГ у крыс с генети-
ческой предрасположенностью к абсанс-эпилепсии (линия WAG/Rij) // Журнал высшей нервной деятельности им. И.П. Павлова. 2012. Т. 62, № 6. С. 733.
16. Грубов В.В., Ситникова Е.Ю., Назимов А.И., Руннова А.Е., Храмова М.В., Храмов А.Е. Возрастная динамика частотно-временных особенностей сонных веретен на ЭЭГ крыс с генетической предрасположенностью к эпилепсии // Вестник Тамбовского университета. Серия: Естественные и технические науки. 2013. Т. 18. № 4-1. С. 1288.
17. Грубов В.В., Овчинников А.А., Ситникова Е.Ю., Короновский А.А., Храмов А.Е. Вейвлетный анализ сонных веретен на ЭЭГ и разработка метода их автоматической диагностики // Изв. вузов. Прикладная нелинейная динамика. 2011. Т. 19, №4. С. 91.
18. Безручко Б.П., Смирнов Д.А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ «Колледж», 2005. 320 с.
19. Gideon Schwarz. Estimating the dimension of a model // The Annals of Statistics. 1978. Vol. 6, № 2. P. 461.
20. Кобзарь А.И. Прикладная математическая статистика. М.: Физматлит, 2006. С. 658.
21. Larose D. T. Data mining methods and models. Wiley: IEEE Press, 2006.
22. Гмурман В.Е. Теория вероятностей и математическая статистика: Уч. пособ. для вузов. Изд. 4-е доп. М.: Высшая школа, 1972. С. 327.
Саратовский госуниверситет Поступила в редакцию 8.04.2013
им. Н.Г. Чернышевского После доработки 15.07.2013
OPTIMIZATION OF BASIS FUNCTION SET FOR MODEL MAP RECONSTRUCTION OF SHORT ELECTROENCEPHALOGRAM TRACINGS DURING EPILEPTIC SEIZURE
A. N. Severyuhina
The problem of compact mathematical model reconstruction of short electroencephalogram tracings during epileptic seizure is solved. This kind of model map can be usefUl in many applications, for example, in time series segmentation with following clustering of obtained fragments. Optimization methods are proposed as a solution. It is shown that application of optimization methods allows to obtain adequate model at that time decreasing number of modeling map basis functions.
Keywords: Model map, experimental data, time series, optimization of basis function set, electroencephalogram, absence seizure.
Северюхина Александра Николаевна - родилась в Саратове (1988). Окончила Саратовский государственный университет (2011). Аспирант кафедры динамического моделирования и биомедицинской инженерии СГУ. Область научных интересов - численные методы, моделирование по временным рядам, нелинейная динамика.
410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: severyuhina_alexandra@mail.ru
