Хабрахабр

Простой классификатор P300 на открытых данных

Сегодня в продолжение темы нейроинтерфейсов мы используем один из открытых датасетов, записанных на игре, использующей механику P300, чтобы визуализировать сигнал ЭЭГ, посмотреть структуру вызванных потеницалов, построить основные классификаторы, оценить качество, с которым мы можем предсказать наличие такого вызыванного потенциала. Мой коллега Рафаэль Григорян eegdude недавно написал статью о том, зачем человечеству потребовалась ЭЭГ и какие значимые явления могут быть зарегистрированы в ней.

Обычно он используется для построения современных BCI. Напомню, что P300 — это вызванный потенциал (ВП), специфический отклик мозга связанный с принятием решений и и различением стимулов (что он из себя представляет мы увидим ниже).

Для того, чтобы заняться классификацией ЭЭГ, можно позвать друзей, написать игру про Енотов и Демонов в VR, записать собственные реакции и написать научную статью (об этом я расскажу как-нибудь в другой раз), но по счастью, учёные со всего мира уже провели некоторые эксперименты за нас и осталось только скачать данные.

Разбор способа построения нейроинтерфейса на P300 с пошаговым кодом и визуализациями, а также ссылку на репозиторий можно найти под катом.

В статье приведены только основные моменты из кода, полную воспроизводимую версию в jupyter notebook искать тут

Способов вызвать его известно множество, например, если концентрироваться на одном предмете, а он в случайный момент активируется (изменит форму, цвет, яркость или отпрыгнет куда-то). С точки зрения ЭЭГ P300 это всего лишь всплеск в определённое время в определённых каналах. Вот как это было реализовано в стародавние времена

Человек выбирает один из них и сосредотачивается на нём (хороший способ — считать количество активаций), далее каждый объект мигает в случайном порядке. В общем виде схема следующая: в поле зрения человека есть несколько (обычно от 3 до 7) стимулов. Поскольку человек концентрировался только на одном стимуле, то и пик должен быть один. Зная время активации каждого стимула мы теперь можем посмотреть на следующую за этим ЭЭГ и установить, был ли в ней характерый пик (мы увидим его на визуализациях ниже). Если вариантов больше семи, можно положить их на сетку и свести задачу к выбору строки+столбца. Таким образом в этих ваших нейроинтерфейсах происходит выбор одного из нескольких вариантов (букв для написания, действий в игре и Бог знает чего ещё). Так получился классический матричный P300 speller, показанный выше.

Выглядело это примерно так В случае датасета, рассматриваемого сегодня, визуальная часть (как и название) была позаимствованы у известной игры space inviders.

Brain Inviders gameplay

По-сути, это тот же спеллер, только буквы заменены на игровых инопланетян.
Также сохранилось видео процесса игры и технические отчёты

Данные состоят из 16 каналов ЭЭГ и одного канала событий, показывающего в какие моменты были активированы целевой (загаданный игроком), и нецелевые стимулы, с ними мы и будем работать. Так или иначе, в интернете появились данные, собранные с помощью этой игры и мы можем получить к ним доступ.

Большинство датасетов для BCI было записано нейрофизиологами, а это ребята, которые не очень заботятся о совместимости, поэтому форматы данных отличаются большим разнообразием: от разных версий .mat файлов до "стандартных" форматов .edf и .gdf.
Самое главное, что нужно знать про эти форматы — вы не хотите их парсить или работать с ними напрямую.
По счастью группа энтузиастов из NeuroTechX написала загрузчики для некоторых датасетов напрямую в numpy.
Эти загрузчики являются частью проекта moabb претендующего на универсализацию решений для BCI.

Загрузка сырого датасета

import moabb.datasets
sampling_rate = 512 m_dataset = moabb.datasets.bi2013a( NonAdaptive=True, Adaptive=True, Training=True, Online=True,
)
m_dataset.download()
m_data = m_dataset.get_data()

Это структура из пакета mne, им обычно пользуются биологи для взаимодействия с сигналами: в этой структуре есть встроенные методы для фильтрации, визуализации, хранения меток и мало ли чего ещё. На данном этапе мы получили структуру RawEDF, содержащую записи ЭЭГ. интерфейс пакета имеет свойство быть нестабильным (текущая версия 0. Но мы не будем идти этим путём т.к. 17 т.к. 19, но мы будем использовать 0. датасет уже не читается новой версией) и слабо документированным, через это наши результаты могут стать невоспроизводимыми.

Это международная схема расположения электродов на голове человека, созданная для того, чтобы учёные могли соотносить зоны мозга и расположения каналов ЭЭГ. Что мы возьмём из полученной структуры — это метки каналов в системе 10-20. Ниже приведено расположение электродов в системе 10-10 (отличается от 10-20 вдвое большей плотностью разметки) и красным отмечены те каналы, которые были записаны в этом датасете.
система 10-10

print(m_data[1]['session_1']['run_1'])
# <RawEDF | 1.gdf, n_channels x n_times : 17 x 159232 (311.0 sec), ~20.7 MB, data loaded> channels = m_data[1]['session_1']['run_1'].ch_names[:-1]
channels # ['FP1', 'FP2', 'F5', 'AFz', 'F6', 'T7', 'Cz', 'T8', 'P7', 'P3', 'Pz', 'P4', 'P8', 'O1', 'Oz', 'O2']

Сначала мы из загруженных данных для каждого испытуемого выделяем массивы непрерывного ЭЭГ по 16 секунд и все метки для этого промежутка (в данных это просто ещё один канал, в котором отмечены начала интересующих нас событий).

На данном этапе мы сохраняем максимальную длину непрерывного ЭЭГ для того, чтобы при дальнейшей фильтрации не встречаться с краевыми эффектами.

raw_dataset = [] for _, sessions in sorted(m_data.items()): eegs, markers = [], [] for item, run in sorted(sessions['session_1'].items()): data = run.get_data() eegs.append(data[:-1]) markers.append(data[-1]) raw_dataset.append((eegs, markers))

Фильтрация и разделение на эпохи

Также не очень давно вышел более свежий обзор нейросетевых метедов. В целом для обзора методов предобработки и классификации ЭЭГ могу порекомендовать отличный обзор от мэтров нейрокомпьютерных интерфейсов.

Минимальная предобработка сигнала ЭЭГ для классификации включает 3 шага:

  • децимация
  • фильтрация
  • масштабирование

Для реализации этих шагов мы воспользуемся старым добрым sklearn-ом и его парадигмой трансформеров и пайплайнов для того, чтобы наш препроцессинг мог быть легко расширяемым.
Код трансформеров вынесен в отдельный файл, ниже опишем некоторые детали.

Децимация

Это совершенно некорректно (почему — см. Почему-то в некоторых статьях и примерах обаработки я встречал понижение частоты сигнала простым выкидыванием отсчётов в стиле eeg = eeg[:, ::10]. Мы используем стандартную реализацию из scipy. любую книжку по обработке сигналов).

Фильтрация

Частоты отсечки — от 0. Здесь мы также опираемся на фильтры из scipy, выбрав bandpass фильтр Баттерворта 4 порядка и применив его в прямом и обратном направлении (filtfilt) для сохранения фазы. 5 до 20 Гц, это стандартный диапазон для нашей задачи.

Масштабирование

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

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

reload(transformers) decimation_factor = 10
final_rate = sampling_rate // decimation_factor
epoch_duration = 0.9 # seconds
epoch_count = int(epoch_duration * final_rate) eeg_pipe = make_pipeline( transformers.Decimator(decimation_factor), transformers.ButterFilter(sampling_rate // decimation_factor, 4, 0.5, 20), transformers.ChannellwiseScaler(StandardScaler()),
)
markers_pipe = transformers.MarkersTransformer(labels_mapping, decimation_factor)

Будем называть эпохой промежуток времени сразу после активации стимула с характерной продолжительностью 0. Далее мы применим пайплайн препроцессинга к нашим данным и порежем непрерывный сигнал ЭЭГ на эпохи. 5-1 секунды, в нашем случае продолжительность составляет 900 мс, хотя её можно и сократить.

В нашем датасете 16 каналов ЭЭГ, после применения децимации частота понизится до 50 Гц, таким образом одна эпоха будет описываться матрицей (16, 45) — 900 мс на 50 Гц это 45 временных отсчётов.

Метки в данном датасете только бинарыне — они маркируют целевой (загаданный игроком, активный, 1) и не целевой (пустой, 0) сигнал.

for eegs, _ in raw_dataset: eeg_pipe.fit(eegs) dataset = [] for eegs, markers in raw_dataset: epochs = [] labels = [] filtered = eeg_pipe.transform(eegs) markups = markers_pipe.transform(markers) for signal, markup in zip(filtered, markups): epochs.extend([signal[:, start:(start + epoch_count)] for start in markup[:, 0]]) labels.extend(markup[:, 1]) dataset.append((np.array(epochs), np.array(labels)))

dataset[0][0].shape, dataset[0][1].shape
# ((1308, 16, 45), (1308,))

С такой структурой мы можем как провести кросс-валидацию внутри данных одного человека, так и протестировать переносимость классификатора между разными людьми (т.н. Таким образом у нас получился датасет в стиле Pytorch в котором первый индекс отсчитывает различных людей. Данные одного человека состоят из массива эпох и меток классов. transfer learning, calibration-less prediction). Количество эпох у каждого человека немного разнится в силу особенностей записи.

Исследование и визуализация данных

Для начала взглянем на один из непрерывных сигналов до нарезки на эпохи.

Не смотря на то, что он уже отфильтрован, на нём не видно глазом каких-либо активаций и больше походит на какой-то шум.

Filtered EEG signal of one run of the game

Это и есть наш искомый вызванный потенциал P300. Если же мы рассмотрим только одну целевую эпоху из нашего датасета, то увидим характерный подъём на промежутке 400-600 мс.

p300 epoch

У каждого человека примерно от 1300 до 1750 (это связано с тем, что кто-то сбивал пришельцев быстрее, а кто-то медленнее). Всего в нашем датасете около 35 тысяч эпох, то есть активаций стимулов.

у нас в матрице 6 строчек и столбцов и только по одному из них — целевые. Также в классах есть ощутимый дисбаланс: 1 к 5 в пользу пустых стимулов т.к. Позже мы вернёмся к этому при обсуждении полученных метрик.

class balance

Теперь настало время взглянуть на отличие целевого сигнала от нецелевого

dataset mean response

На левом графике можно увидеть, что средние сигналы разнятся очень сильно, причём у обоих есть неспецифический отклик в районе 180мс, но у целевого он гораздо амплитуднее, также у целевого есть характерный горб от 250 до 500 мс — это и есть пресловутый P300.

И это не смотря на то, что соотношение сигнал-шум для P300 считается одним из самых высоких в нейрофизиологии. С таким отличием в сигнале наша задача может показаться плёвой, но если мы добавим на график стандартное отклонение в каждой точке, то увидим, что картина не так уж радужна — сигнал довольно сильно зашумлён.

пустой сигнал усредняется по впятеро большему числу разных семплов, таким образом случайные отклонения давятся больше, но как мы видим по дисперсии одного порядка, это не слишком помогает) (На самом деле эти графики построены не вполне честно т.к.

Person mean response

Также полезно взглянуть на средние сигналы одного человека.

Также пик P300 у одного человека более высокий из-за меньшего усреднения. Тут находит подтвержедние предыдущая ремарка про "нечестное" усреднение — пустой сигнал заметно амплитуднее, чем при усреднении по всем.

Межличностная вариативность нейрофизиологических реакций довольно высока, мы ещё увидим влияние этого фактора в работе классификаторов. Важно отметить ещё одну особенность сигнала одного человека — он имеет несколько другую форму, чем обобщённый. Впрочем, внутриличностные различия (один человек в разном настроении, уровне стресса, усталости) также достаточно велика.

png

Точка зрения тут совпадает с картинкой выше, на которой изображены положения электродов — нос сверху и т.д. Далее мы видим поканальную развёртку сигналов.

На Fp1,2 ярко выражены два отрицательных пика, предшествующих положительному пику. Отклик каждой части головы разный. Также в некоторых каналах положительных пика два, а в некоторых — один или что-то переходное между.

Примененнием этих методов мы займёмся в другой раз. Разные каналы имеют разную важность для определения наличия P300, её можно оценить разными методами — вычислением взаимной информации (mutual information) или методом добавления-удаления (aka stepwise regression).

Понятно, что при наличии 16 электродов точность такой карты оставляет желать лучшего, но какое-то понимание сформировать. Стоит помнить, что электродами мы измеряем разность потенциалов между электродами, а значит мы можем по измененрям напряжения в отдельных точках построить карту напряжений для всей головы в определённые моменты времени. (mne по умолчанию ожидает увидеть микровольты, но мы уже применили масштабирование, так что абсолютные величины указаны не верно)

png

png

Классификация

Наконец пришло время применить методы машинного обучения к нашей выборке.

регрессия, метод опорных векторов (SVM) и несколько методов, использующих корреляционный анализ из пакета pyriemann (детали работы каждого метода можно найти в документации), стоит отметить, что эти методы специально разрабатывались для применения к ЭЭГ и с их помощью было выиграно несколько соревнований на kaggle. В качестве классификаторов были выбраны несколько базовых — лог.

clfs = , ), 'LDA': ( make_pipeline(Vectorizer(), LDA(shrinkage='auto', solver='eigen')), {}, ), 'SVM': ( make_pipeline(Vectorizer(), SVC()), {'svc__C': np.exp(np.linspace(-4, 4, 9))}, ), 'CSP LDA': ( make_pipeline(CSP(), LDA(shrinkage='auto', solver='eigen')), {'csp__n_components': (6, 9, 13), 'csp__cov_est': ('concat', 'epoch')}, ), 'Xdawn LDA': ( make_pipeline(Xdawn(2, classes=[1]), Vectorizer(), LDA(shrinkage='auto', solver='eigen')), {}, ), 'ERPCov TS LR': ( make_pipeline(ERPCovariances(estimator='oas'), TangentSpace(), LogisticRegression()), {'erpcovariances__estimator': ('lwf', 'oas')}, ), 'ERPCov MDM': ( make_pipeline(ERPCovariances(), MDM()), {'erpcovariances__estimator': ('lwf', 'oas')}, ),
}

crossvalidate record

сначала нужно, чтобы человек какое-то время концентрировался на заранее указанных стимулах и только после этого мы предсказываем его выбор. Наиболее часто встречающаяся схема нейроинтерфейсов это "калибровка+работа" т.е. Этот подход обладает очевидным минусом занудого начального этапа.

Для оценки работоспособности наших методов в этом режиме проведём кроссвалидацию внутри эпох одного человека.
Метрика accuracy в данном случае не релевантна в силу разбалансированности датасета (бейзлайн тут 5/6 ~ 83%), так что я предпочитаю смотреть на тройку precision-recall-f1.

В целом, перфоманс лучших моделей достаточно высок по сравнению с тем, что мы в Neiry имеем в "полевых" условиях парка развлечений (напомню, что этот датасет записывался в лаборатории). Чтобы обозреть весь датасет, усредним результаты такой кроссвалидации по всем людям.

В целом же нам нужно решить мультиклассовую задачу выбора одного из стимулов (она, кстати, сбалансирована т.к. В этом датасете присутствуют только бинарные метки для данных. Для её решения обычно фиксируют число активаций каждого стимула (например, 6 стимулов по 5 активаций) и все стимулы в случайном порядке активируются (30 раз), получается 30 эпох и для каждого стимула складываются вероятности его активаций быть целевыми, после этого стимул, набравший максимальную сумму признаётся целевым. каждый стимул активируется одиаковое количество раз). Релизацию этого подхода мы продемонстрируем в одном из будущих постов на подходящем датасете.

crossvalidate_dataset

Дело в том, что когда мы делаем калибровку, то мы фактически переобучаемся на форму пика одного человека, поэтому можем неплохо предсказывать её же в последующих тестах. Вторая схема проведения называется transfer learning — то есть перенос классификатора между людьми. В случае отсутствия калибровки, предобученный классификатор должен уметь выделять концепцию P300 не зная наперёд форму сигнала конкретного человека.

Мы проведём два эксперимента — обучим классификатор на одном человеке, и предскажем пятерых, а потом увеличим тренировочную выборку до 10 человек и сравним результаты чтобы убедиться, что модели смогли повысить свою обобщающую способность

Тренировка на 1 человеке
transfer_validate_1

Тренировка на 10 людях
transfer_valiidate_10

23 до 0. Итак f1 поднялся с 0. 4 для лучшего классификатора (в обоих случаях это логрегрессия с одинаковой регуляризацией).

Исходя из нашего опыта при таких метриках бинарной задачи достаточно 5 активаций каждого стимула, для достижения accuracy мультиклассовой задачи около 75%. Это означает, что предсказательная способность повысилась от "никакой" до "приемлемой".

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

Заключение

Рекомендую заинтересовавшимся самостоятельно открыть ноутбук (находится в репозитории) и поэкспериментировать с вариантами визуализаций и классификаторов. Сегодня мы ближе познакомились с вызванным потенциалом P300 и построили простой пайплайн для нейроинтерфейса.

Продолжение следует... Имея базовое представление о методах работы с сигналом ЭЭГ мы в дальнейшем сможем более глубоко рассмотреть эту тему — применить продвинутые методы предобработки, а также нейросети для решения задач построения нейроинтерфейсов.

Показать больше

Похожие публикации

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»