В.И. СКРУГИН, А.Г. ТРОФИМОВ, А.О. РОИК¹,

Р.А. НАУМОВ*

Национальный исследовательский ядерный университет «МИФИ», Москва

¹Институт высшей нервной деятельности и нейрофизиологии РАН, Москва

atrofimov@list.ru
АЛГОРИТМ КЛАССИФИКАЦИИ СИГНАЛОВ ЭЭГ

НА ОСНОВЕ АНАЛИЗА В ЧАСТОТНО-ВРЕМЕННОЙ

ОБЛАСТИ
Рассматривается задача классификации электроэнцефалограмм, полученных при решении испытуемым задач на разные типы мышления. Предложен алгоритм расчёта вектора характерных признаков электроэнцефалограммы на основе результатов непрерывного вейвлет-преобразования. Для классификации этих векторов предложено использовать многослойный персептрон. Приводятся результаты экспериментальных исследований построенного классификатора и оценка его точности.
Введение
Электроэнцефалография является одним из основных неинвазивных методов регистрации электрической активности нейронов мозга. В последнее десятилетие наблюдается рост интереса к анализу нейрофизиологических данных, прежде всего, к электроэнцефалограмме (ЭЭГ) и томограмме, что связано с попыткой использования этой информации об активности нервных клеток мозга для управления техническими устройствами. Такого рода устройства, называемые нейрокомпьютерными интерфейсами (Brain-Computer Interface, BCI), могут применяться в медицине, военной сфере, промышленности, а также являться основой коммуникационных систем нового поколения.

В основе систем нейрокомпьютерного интерфейса лежит математический анализ нейрофизиологических данных и алгоритмы их классификации. Количественным анализом сигналов ЭЭГ с применением математических методов и вычислительной техники занимается отдельное направление электроэнцефалографии – компьютерная электроэнцефалография [1]. Обычно в компьютерной электроэнцефалографии задача распознавания сигнала ЭЭГ включает два этапа:

1) выделение характерных признаков сигнала ЭЭГ;

2) классификация вектора выделенных характерных признаков.

При решении первой задачи наибольшее распространение получили методы спектрального анализа, в частности, преобразование Фурье, с помощью которого осуществляется переход от временного представления сигнала к частотному.

Недостатком спектрального анализа является потеря информации о временной локализации частот, что особенно важно при анализе нестационарных сигналов. Частично этот недостаток может быть устранён при использовании оконного преобразования Фурье. Развитием оконных методов спектрального анализа являются вейвлет-преобразования [2].

В настоящей работе предлагается подход к расчёту вектора характерных признаков ЭЭГ на основе вейвлет-анализа сигналов. В связи с тем, что характер ритмической активности ЭЭГ специфичен для каждого индивида [3], классификатор векторов характерных признаков ЭЭГ должен адаптивно настраиваться на индивидуальные особенности конкретного испытуемого. Учитывая это, для дальнейшей классификации характерных признаков ЭЭГ предложено использовать нейросетевой подход [4].
Исходные данные и постановка задачи
В качестве исходных данных для классификации были использованы электроэнцефалограммы, полученные в лаборатории высшей нервной деятельности человека Института высшей нервной деятельности и нейрофизиологии Российской академии наук. Каждая ЭЭГ соответствует процессу решения испытуемым задачи одного из четырёх типов, связанных с различными типами мышления: планиметрическое мышление (тип ‘G’), пространственное мышление (тип ‘C’), составление предложений (тип ‘S’) и исключение лишнего слова из перечня (тип ‘W’).

Испытуемый решил примерно по 70 задач каждого типа, потратив на каждую от четырех секунд до нескольких десятков секунд. Все сигналы ЭЭГ были оцифрованы с частотой дискретизации _д = 250 Гц и очищены от артефактов. Запись ЭЭГ велась по 22 каналам, из которых информативными являются L = 19.

Таким образом, результат i-го эксперимента k-го класса представляет собой многомерный временной ряд , , где t – номер дискретного отсчёта, продолжительностью от одной тысячи до нескольких тысяч тактов.

В работе рассматривается задача построения классификатора многомерных временных данных ЭЭГ , , , где N_k – число сигналов ЭЭГ k-го класса, K – число классов. Первым шагом при решении этой задачи является формирование вектора признаков , характеризующего i-й сигнал ЭЭГ k-го класса, , .

Формирование моментных характерных признаков ЭЭГ
Вектором моментных характерных признаков (МХП) v(t) будем называть вектор, отражающий состояние сигнала ЭЭГ x(t), а следовательно, мозга испытуемого, в некоторый фиксированный момент времени t. Учитывая возможности применения вейвлетов для анализа в частотно-временной области, предложено формировать вектор МХП на основе результатов вейвлет-преобразования сигналов ЭЭГ. Согласно [2], для анализа временных рядов наиболее эффективным оказывается непрерывное вейвлет-преобразование (НВП), позволяющее выделить локализованные во времени (моментные) гармонические составляющие сигнала.

На практике расчеты вейвлет-коэффициентов сигнала x(t), заданного дискретным набором значений, t = 0,…,T–1, производят по формуле [2]:

, (1)

где T – интервал дискретности, a, b – параметры масштаба и временного сдвига соответственно, ψ(t) – материнская вейвлет-функция, оператор * означает комплексное сопряжение. В качестве материнского вейвлета ψ(t) выбран комплексный вейвлет Морле, имеющий следующий аналитический вид:

.

В [2] показано, что НВП синусоиды вейвлетом Морле представляет собой гряду локальных максимумов, ориентированную вдоль оси времени. Причем ордината, общая для всех максимумов, равна периоду детектируемой гармоники. Сумма квадратов абсолютных значений вейвлет-коэффициентов служит оценкой энергии сигнала в соответствующем частотном и временном диапазоне.

Для представления результатов НВП сигнала l-го канала ЭЭГ в памяти компьютера используется матрица вейвлет-коэффициентов, рассчитанных в узлах дискретной сетки на плоскости (a, b): , . Величина шага вдоль оси b выбрана равномерной. Учитывая обратно пропорциональную связь между значением масштаба a и соответствующей ей частотой , изменение шага сетки вдоль оси a выбрано гиперболическим. Значения исследуемых масштабов определяются по формуле:

,

где [ω_min, ω_max] – исследуемый частотный диапазон, N_a – разрешение матрицы W по масштабу. Таким образом, рассчитанной сетке на оси масштабов будет соответствовать равномерная сетка на оси частот.

Рис. 1. Вейвлетограмма сигнала ЭЭГ в диапазоне частот от 5 Гц до 30 Гц;

по оси абсцисс отложены значения временных сдвигов b, по оси ординат –

масштабных сдвигов a
На рис. 1 показан пример визуального представления матрицы квадратов вейвлет-коэффициентов (вейвлетограммы), рассчитанной для одного из каналов ЭЭГ. Каждый столбец матрицы отражает распределение энергии сигнала по частотам в фиксированный момент времени (моментную спектральную плотность мощности). Каждая строка отражает распределение энергии сигнала на фиксированной частоте по времени.

Поскольку сигналы ЭЭГ являются многомерными, вектор МХП ЭЭГ v(t) составлен из значений энергии сигнала каждого из каналов ЭЭГ в заранее заданных частотных диапазонах ₁,… _M в момент времени t:

, ,

. (2)

Размерность вектора МХП v(t) равна M*L, где M – число выбранных частотных полос, L – число каналов.

Формирование вектора характерных признаков ЭЭГ
В каждый момент времени t i-й сигнал ЭЭГ k-го класса характеризуется вектором МХП , . Таким образом, этот сигнал может быть описан множеством характерных векторов . Обозначим пространство этих векторов через V.

В работе выдвигается предположение, что распределения векторов МХП ЭЭГ в пространстве V в пределах одного класса схожи, в то время как для ЭЭГ из разных классов имеют значимые различия. Подтверждением этому могут служить результаты визуального анализа распределений проекций векторов МХП на различные оси и плоскости. На рис. 2 представлены примеры некоторых из этих распределений.

Для выделения областей в пространстве V, характерных для сигналов ЭЭГ одного класса и нехарактерных для сигналов ЭЭГ всех остальных классов, а также областей, общих для всех классов, используется кластерный анализ. Выделение кластеров должно проводиться среди векторов МХП всех имеющихся сигналов ЭЭГ всех классов. Именно в этом пространстве могут быть обнаружены наиболее характерные области.

Пусть в результате кластеризации множество векторов разбито на кластеры , где P – заранее заданное число. Далее возникает задача определения того, какие из этих кластеров являются наиболее характерными для сигналов ЭЭГ каждого из классов.

(а) (б)

Рис. 2. Распределение проекций векторов МХП двух ЭЭГ (обозначены ‘x’ и ‘о’)

из одного (а) и разных классов (б) на плоскость

Каждый кластер представляет собой некоторую область в пространстве V. Будем говорить, что кластер _p активен в момент времени t, если вектор МХП ЭЭГ . Таким образом, каждый сигнал ЭЭГ может быть описан последовательностью номеров активных кластеров .

Для определения кластеров, специфичных для сигналов ЭЭГ каждого класса, могут быть предложены различные критерии, связанные с анализом активности кластеров. Введём следующие временные характеристики p-го кластера, :

– время (число тактов дискретного времени), в течение которого
p-й кластер был активен для i-го сигнала ЭЭГ k-го класса;

– относительное время активности p-го кластера для i-го сигнала ЭЭГ k-го класса – индекс активности p-го кластера.

Введём критерий, позволяющий выявить кластеры, для которых наблюдается наибольшее различие между сигналами ЭЭГ двух классов:

, (3)

где n₁, n₂ – число экспериментов 1-го класса, отнесённых ко 2-му классу, и наоборот соответственно, при использовании оптимального линейного классификатора. Критерий оптимальной линейной классификации сигнала ЭЭГ x, имеющего индекс активности в p-м кластере _p, имеет вид:

где значение определяется, исходя из минимально возможного значения критерия I(p) при фиксированном p.

Фактически значение I(p) для p-го кластера характеризует вероятность ошибочной линейной классификации сигнала ЭЭГ на основе значения индекса активности данного кластера. Индекс активности кластера, лучшего с точки зрения этого критерия, выбран в качестве характерного признака сигнала ЭЭГ для случая двух классов. Несмотря на то, что критерий (3) может быть легко обобщён на случай K классов, более перспективным представляется формирование вектора характерных признаков ЭЭГ на основе результатов попарного разделения классов.

Пусть – номер кластера, оптимального с точки зрения критерия (3) при разделении классов k₁ и k₂, , . Тогда вектор характерных признаков ЭЭГ составим из индексов кластеров, оптимальных при разделении всех пар классов:

. (4)

Если при разделении различных пар классов оказывается оптимальным один и тот же кластер, то повторное значение из вектора z удаляется. Таким образом, максимальная размерность вектора z равна K*(K 1)/2.

Для построения классификатора векторов характерных признаков z предложено использовать нейросетевой подход. В качестве архитектуры нейросетевого классификатора выбран многослойный персептрон. Обучающая выборка для него имеет вид , где  – вектор характерных признаков вида (4) сигнала ЭЭГ , подаваемый на вход персептрона,  – индикаторный вектор выходов:

Общее число примеров, используемых для обучения и тестирования классификатора, равно числу имеющихся сигналов ЭЭГ всех классов.

Результаты экспериментальных исследований
Исследовались два диапазона частот: от 7 до 12 Гц, что соответствует диапазону альфа-ритма, и от 5 до 30 Гц. В качестве элементов вектора МХП выбраны моментные значения энергии в частотных полосах шириной 1 Гц в рассматриваемых интервалах, рассчитанные по формуле (2). В первом случае размерность вектора МХП v(t) равна 19*5=95, во втором – 19*25=475.

Для кластеризации векторов МХП ЭЭГ применён алгоритм нейронного газа [4]. В процессе самообучения уровень соседства экспоненциально менялся от λ_max = 100 до λ_min = 0,001, коэффициент обучения линейно менялся от μ_max = 1 до μ_min  = 0,001. Графики распределения значений критерия I(p) по кластерам для различных пар классов ЭЭГ показаны на рис. 3.

(а) (б)

Рис. 3. Распределение значения критерия I(p) по кластерам: а) диапазон исследуемых частот [ω_min, ω_max] = [7;12] Гц, P = 50, б) [ω_min, ω_max] = [5;30] Гц, P = 200. Сплошная линия соответствует классам ‘C’, ’G’, пунктирная – ‘C’, ’S’, точечная – ‘G’, ’S’. Номера кластеров упорядочены по значению критерия
Из рисунка видно, что наилучшее разделение наблюдается для классов ‘C’ и ‘S’ при исследовании обоих частотных диапазонов, причём при использовании единственного признака (значение активности 1 го кластера) погрешность оптимального линейного классификатора составляет около 90% на обучающей выборке. Отметим, что для некоторых кластеров индекс активности всех сигналов ЭЭГ рассматриваемой пары классов равен нулю («пустые кластеры»). Для таких кластеров значение критерия I(p) не определено.

На рис. 4 показано распределение значений индекса активности  оптимального с точки зрения критерия (3) кластера для ЭЭГ классов ‘C’ и ‘S’. Из графика видно, что для ЭЭГ класса ‘S’ индекс активности рассматриваемого кластера в большинстве экспериментов отличается от нуля, в то время как практически для всех ЭЭГ класса ‘C’ он равен нулю, что и позволяет использовать данный показатель при классификации.

(а) (б)

Рис. 4. Распределение значений (полигон частот) индекса активности оптимального кластера, найденного для пары классов 'C' (сплошная линия) и 'S' (пунктирная линия): а) диапазон исследуемых частот [ω_min, ω_max] = [7;12] Гц, б) [ω_min, ω_max] = = [5;30] Гц. По оси абсцисс отложены значения индекса активности, по оси ординат – доля экспериментов, для которых индекс активности принимает соответствующие значения
В табл. 1 приведены характеристики оптимальных кластеров при P = 200, исследуемый диапазон частот [ω_min, ω_max] = [5; 30] Гц. Индексы активности этих кластеров далее выбраны в качестве элементов вектора характерных признаков z.

Для классификации вектора z использован многослойный персептрон с числом входов – 4, числом выходов – 4. В результате исследований установлено, что персептрон с 2 мя скрытыми слоями, число нейронов в скрытых слоях – 15 и 12, показывает лучшие результаты по точности. Активационные характеристики нейронов скрытых слоёв – сигмоиды, нейронов выходного слоя – линейные. Из всех обучающих векторов для обучения случайно отобраны 70%, остальные примеры использованы для тестирования. Показатели точности обученной сети на обучающей и тестовой выборках приведены в табл. 2.

Таблица 1
Характеристики оптимальных кластеров

№ п/п	Пара классов	Опт. кластер, p	I(p)	0p
1	'C', 'G'	143	0.155	0.001148
2	'C', 'S'	143	0.089	0.001297
3	'C', 'W'	103	0.158	0.015991
4	'G', 'S'	64	0.206	0.018765
5	'G', 'W'	64	0.324	0.037332
6	'W', 'S'	89	0.253	0.005097

Таблица 2

Показатели точности обученной сети

Выборка	‘C’	‘G’	‘S’	‘W’
обуч.	95%	91%	98%	92%
тест.	80%	78%	75%	74%

Проведены исследования зависимости точности обучения классификатора от числа кластеров P. Для каждого из значений P, приведённых в табл. 3, был сформирован вектор характерных признаков z, построен свой нейросетевой классификатор и определены показатели его точности на тестовой выборке.

Таблица 3

P	‘C’	‘G’	‘S’	‘W’
10	77%	68%	78%	74%
50	77%	60%	69%	69%
100	75%	61%	72%	61%
200	80%	78%	75%	74%
500	68%	60%	64%	66%

Заключение
В результате проведённых исследований установлено, что предложенный алгоритм распознавания сигналов ЭЭГ, основанный на анализе результатов непрерывного вейвлет-преобразования, показывает на имеющихся данных точность около 75%. Предложенный подход к формированию вектора характерных признаков ЭЭГ и их классификации являеюся основой для дальнейших теоретических исследований.

Работа выполнена при финансовой поддержке Совета по Грантам Президента РФ для поддержки молодых российских учёных (Грант МК 4245.2009.8).

Список литературы

1. 
   Зенков Л.Р. Клиническая электроэнцефалография (с элементами эпилептологии). Руководство для врачей. – М.: МЕДпресс-информ, 2004.
   
2. 
   Витязев В.В. Вейвлет-анализ временных рядов: Учеб. пособие. – СПб.: Изд-во
   С.-Пб. ун-та, 2001.
   
3. 
   Иваницкий Г.А. Распознавание типа решаемой в уме задачи по нескольким секундам ЭЭГ с помощью обучаемого классификатора // Сборник научных трудов VIII Всероссийской научно-технической конференции "Нейроинформатика-2006". – М.: МИФИ, 2006. Т.3. С.217-224.
   
4. 
   Осовский С. Нейронные сети для обработки информации. – М.: Финансы и статистика, 2002.