CC BY
61
электронное научно-техническое издание
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС 77 - 30569. Государственная регистрация №0420900025. ISSN 1994-040S
Адаптивный классификатор многомерных нестационарных сигналов на основе анализа динамических паттернов
# 08, август 2010
авторы: Трофимов А. Г., Скругин В. И.
УДК 519.25+004.9
atrofimov@¡list. т
Национальный исследовательский ядерный университет "МИФИ" Введение
Задачи анализа и классификации многомерных динамических данных получили широкое распространение в различных областях науки и техники, среди которых медицинская диагностика, распознавание речи, идентификация динамических систем, эконометрика и другие. В ряде случаев исследователю, работающему в этих областях, приходится иметь дело с принципиально нестационарными динамическими данными. Существующие методы цифровой обработки сигналов в основном ориентированы на анализ и обработку одномерных временных рядов, в то время как многомерным динамическим данным уделяется существенно меньшее внимание. Методы анализа и классификации, ориентированные на одномерные сигналы, зачастую оказываются неэффективными для многомерных данных или вообще не могут быть обобщены на многомерный случай. С другой стороны, нестационарность временных рядов накладывает до-
полнительные ограничения на возможность применения известных методов цифровой обработки сигналов. Как правило, методы обработки одномерных нестационарных сигналов сводятся к использованию скользящих оконных преобразований (в частности, вейвлет-преобразования) или к сегментации с последующим анализом выделенных квазистационарных участков. [1-3]
Построение классификатора временных рядов предполагает выполнение следующих этапов:
1) Предобработка сигналов. На этом шаге осуществляется фильтрация сигналов от шумов, устранение артефактов или другие операции, позволяющие на ранних этапах отсеять неинформативные данные, используя априорно известную информацию о природе сигнала.
2) Выделение характерных признаков сигнала. Целью этого этапа является переход от исходного пространства динамических данных к пространству многомерных статических данных. Как только такой переход сделан, задача классификации сигналов сводится к известной задаче классификации многомерных данных.
3) Классификация в пространстве характерных признаков сигналов. На этом этапе могут использоваться классические алгоритмы классификации многомерных статических данных, например метод дискриминантных функций, нейросетевой аппарат и другие [4,5].
В настоящей работе вопросы предобработки сигналов рассматриваться не будут. При работе с конкретными временными рядами исследователь, обычно, знает, какие алгоритмы предобработки нужно применить к этим сигналам, чтобы увеличить их информативность. Как правило, это какие-либо стандартные методы обработки сигналов или известные специальные методы для рядов рассматриваемой природы. Наибольший интерес в данной работе представляет проблема выделения характерных признаков сигнала для классификации.
Динамическая система, продуцирующая многомерный временной ряд, может находиться в различных функциональных состояниях или работать в различных режимах. Работа системы в том или ином функциональном состоянии находит отражение в наблюдаемых сигналах в виде паттернов. Как правило, моментам смены функционального состояния системы соответствуют
моменты резкого изменения динамических характеристик наблюдаемого временного ряда, т.е. моменты смены паттернов.
В данной работе предлагается методика расчёта характерных признаков многомерных нестационарных сигналов, основанная на выделении функциональных состояний динамической системы, продуцирующей сигнал. Анализ функциональных состояний и переходов между ними, а также соответствующих им паттернов, может послужить основой при классификации сигналов, а также позволить выявить скрытые закономерности, присутствующие в сигнале, и использовать их для лучшего понимания природы наблюдаемой динамической системы.
В работе рассматривается практическое применение предлагаемой методики для анализа и классификации многомерных нестационарных сигналов электроэнцефалограмм (ЭЭГ). Классификатор сигналов ЭЭГ является основой телекоммуникационных систем нового поколения -интерфейсов "мозг-компьютер" (ИМК). ИМК - системы, позволяющие осуществлять управление внешними техническими устройствами непосредственно с помощью электрической активности нейронов мозга, минуя мышцы [6]. Такие системы могут использоваться в медицинских целях - для организации канала коммутации с парализованными пациентами, а также в быту и промышленности - для предоставления нового канала взаимодействия с внешним миром. Современные разновидности этих устройств, основанные на анализе сигналов ЭЭГ, не поддерживают необходимых для промышленного внедрения скоростей передачи информации от мозга к электронному устройству, что в высокой степени обусловлено недостаточной эффективностью существующих алгоритмов анализа сигналов ЭЭГ, отражающих электрическую активность мозга.
Статья состоит из 4 разделов. В первом разделе приводится математическая постановка задачи классификации многомерных нестационарных сигналов. В разделе 2 вводится понятие динамического паттерна и излагаются теоретические основы предлагаемого алгоритма выделения характерных признаков сигнала для классификации на основе анализа динамических паттернов.
В разделе 3 описывается математическая модель нейросетевого классификатора в пространстве выделенных характерных признаков. Раздел 4 посвящён результатам применения предложенной методики для классификации сигналов ЭЭГ, соответствующих различным типам мыслительной активности человека, и обсуждению полученных результатов.
1. Постановка задачи
В связи с тем, что обработка данных осуществляется преимущественно в рамках классических ЭВМ, будем считать, что рассматриваемые сигналы либо имеют дискретную природу, либо дискретизованы. В этом случае значение каждого многомерного сигнала в момент времени t является вектором вида:
х^) = (х ,х2 ,...,х, (ф X, t = й, (1.1)
где t = 1, Т - номер такта дискретного времени, Т - число дискретных отсчетов, Ь - размерность сигнала, х1 (t) - значение 1-ого компонента сигнала |х (t), t = 1, Т|, в момент времени t, I = 1, Ь ,
X е RL - множество значений временного ряда.
Пусть в результате Р наблюдений получено множество многомерных временных рядов в общем случае различной длины, каждый из которых отнесён к одному из К классов. Таким образом, каждый класс представлен конечным числом сигналов. Обозначим Р]с - число наблюдений, относящихся к ^му классу, k = 1, К . Таким образом, общее число рассматриваемых вре-
К
менных рядов Р = ^ Р1с .
k=1
Значение временного ряда, относящегося к ^му классу, полученное в результате р-го наблюдения в момент времени t, будем обозначать
/,р (0 = (х^,р (0,х^,р (0,...,хЬ,р (/))е X, t = 1,т; , р = 1,Рк , k = 1,К (1.2)
где Тр - число дискретных отсчетов р-го сигнала к-го класса, хк,р (г) - значение 1-ого компонента сигнала |хк,р (г), г = 1, Трк |, в момент времени г, I = 1, L .
Классификатор сигналов запишем как оператор, сопоставляющий сигналу {х (г), г = 1, Т | номер класса к, к которому он относится:
С^{х(г),г = 1Т}] = к, к е{1,2,...,К}. (1.3)
Задача построения классификатора временных рядов, таким образом, сводится к задаче нахождения оператора С на основе наблюдений вида (1.2), для которых результат действия этого оператора, т.е. номер класса, известен. В случае классификации статических данных построение оператора С, как правило, сводится к аппроксимации дискриминантных функций классов или границ между классами либо в пространстве самих данных, либо в некотором новом пространстве, рассчитанном в результате предобработки исходных данных [4]. Использование такого подхода к временных рядам, т.е. рассмотрение временной последовательности как вектора признаков, не представляется возможным в силу ряда причин, среди которых высокая размерность этого пространства (длина сигналов может достигать нескольких тысяч), а также различная длина наблюдаемых сигналов. Практически все методы классификации временных рядов предполагают сведение задачи классификации динамических данных к задаче классификации статических данных. Для такого сведения необходимо осуществить переход от временного ряда
{х (г), г = 1, Т} к некоторому вектору признаков V размерности М, характеризующему этот временной ряд. Обозначая оператор этого перехода через запишем
^[{х(г),г = 1Т}] = г, г е ЯМ . (1.4)
Основная проблема состоит в том, что не существует единого алгоритма такого перехода и в каждом конкретном случае исследователь сталкивается с необходимостью поиска оператора ¥
и пространства признаков для эффективной классификации. Особенно остро эта проблема встаёт для случая нестационарных временных рядов.
Учитывая (1.3) и (1.4), задача классификации временных рядов, таким образом, состоит в поиске операторов ¥ и с:
Отметим, что задача нахождения этих операторов на основе результатов наблюдений является неоднозначной. Неточность в оценке оператора ¥ может быть компенсирована за счёт увеличения сложности оператора с, что, однако, может привести к ухудшению обобщающих способностей классификатора.
Качество построенного классификатора напрямую связано с его обобщающими способностями, которые могут быть оценены по результатам классификации данных тестовой выборки, не использовавшихся при нахождении операторов ¥ и с. Также при построении адаптивного классификатора нередко применяется валидационная выборка, назначение которой состоит в контроле процесса настройки классификатора по обучающим примерам, в частности, валида-ционная выборка может участвовать в критерии останова итеративной процедуры обучения. В связи с вышесказанным, перед построением классификатора разобьём множество наблюдаемых временных рядов на обучающую, тестовую и валидационную выборки в заданной пропорции и перенумеруем сигналы внутри этих выборок так, что каждый из них может быть представлен в форме (1.2).
2. Алгоритм выделения характерных признаков многомерного сигнала на основе анализа
(15)
с[г] = ^ k е{1,2,...,К}.
динамических паттернов
2.1. Выделение моментных характерных признаков сигнала
Оператор ¥ в модели (1.5) ставит в соответствие временному ряду {х (г), г = 1, T} некоторый
вектор характерных признаков г, не зависящий от времени г. Это означает, что все интересующие особенности динамики значений временного ряда на рассматриваемом промежутке времени должны быть интегрированы в конечное множество статических показателей. Эти показатели могут быть рассчитаны как по исходным значениям временного ряда, так и по некоторым производным показателям, возможно, лучше отражающим эти особенности.
Введём оператор ф, определённый на множестве рассматриваемых временных рядов вида
(1.1) и возвращающий в каждый момент времени г, г = 1, Т, вектор, характеризующий рассматриваемый временной ряд в этот момент:
ф[{х(т),т = 1Т},г = v(г), v(t)е V, г = 1Т. (1.6)
Вектор v(г), характеризующий временной ряд {х (г), г = 1, Т} в момент времени г, г = 1, Т, будем
называть вектором моментных характерных признаков (МХП). Важно подчеркнуть, что значение вектора v(г) в каждый момент времени г формируется в общем случае на основе значений
сигнала {х (г), г = 1, Т} на всём рассматриваемом интервале времени.
Последовательность векторов ^(г), г = 1, Т} сама является многомерным временным рядом и
определяет фазовую траекторию сигнала {х(г), г = 1,Т} в пространстве V моментных характерных признаков.
Цель перехода от пространства X значений исходного временного ряда к пространству МХП V состоит во введении в рассмотрение дополнительных параметров исходного временного ряда, характеризующих его динамические свойства в каждый момент времени. Выбор конкретного типа оператора ф осуществляется с учетом специфики анализируемых данных. Так, для сигналов гармонической природы в качестве элементов вектора МХП v(г) могут быть выбраны зна-
чения спектральной плотности мощности, рассчитанные по результатам оконного преобразования Фурье с центром окна в момент времени в задачах выделения специфических графоэле-ментов элементами вектора v(t) могут являться коэффициенты непрерывного вейвлет-
преобразования в заданном частотном диапазоне. В простейшем случае v(t) = х^), t = 1, Т.
Отметим, что размерность вектора v(t) может быть выше, чем размерность значений исходного временного ряда х(^.
2.2. Понятие статического и динамического паттерна временного ряда
В связи с тем, что элементы вектора МХП v(t) связаны с динамическими особенностями сигнала х(^ в момент времени t, можно предположить, что близость векторов МХП v(t1) и v(t2) в пространстве V означает близость динамических свойств сигнала в моменты времени Ь и Так, на квазистационарных участках динамические свойства сигнала с течением времени меняются слабо, что приводит к группированию соответствующих значений МХП в некоторой локальной области или появлению регулярности переходов между несколькими локальными областями. В то же время изменение динамических свойств сигнала приведёт к появлению новых локальных областей в пространстве МХП и/или изменению регулярности в переходах между существующими областями.
Отметим, что обратное утверждение, вообще говоря, неверно. Близость динамических свойств сигнала в моменты времени t1 и ^ может не означать близость векторов v(t1) и v(t2) в пространстве МХП. Вектор МХП стационарного сигнала с неизменными динамическими свойствами может иметь в разные моменты времени существенно различные значения МХП.
Пусть пространство Vнекоторым образом разбито на Q областей S1,...,SQ таких, что
Q
П^2 =0, я = 1,Q, 42 = 1,Q, Я *42;
Эти области в пространстве МХП V далее будем называть статическими паттернами (СП). В каждый момент времени I, I = 1, Т, вектор МХП v(t) будет принадлежать одной из областей 51,...,50. Будем говорить, что статический паттерн Sq проявляется в сигнале х(^ в момент времени t, если соответствующий ему вектор МХП v(t) е . Такой паттерн будем называть активным в момент времени t. Таким образом, каждый сигнал ), t = 1,Т| вида (1.1) может быть представлен как последовательность номеров активных статических паттернов
), t = 1Т} .
С течением времени активные статические паттерны для сигнала ), t = 1, Т} некоторым
образом сменяют друг друга. Долю времени, в течение которого СП был активным, будем называть индексом активности р статического паттерна
Т ($а) —
Р(5,) = , , = 1,0, (1.7)
где Т(5,) - число тактов дискретного времени, в течение которого СП был активным.
Таким образом, для временного ряда ), t = 1, Т} может быть составлен вектор
р = (р(51),..., р(5в )) индексов активностей всех рассматриваемых СП.
Классификация временного ряда может быть проведена по результатам анализа индексов активностей СП, присутствовавших в сигнале. Так, среди множества СП могут быть найдены паттерны, специфичные для временных рядов одного класса и неспецифичные для временных рядов всех остальных классов. Эта специфичность паттерна может выражаться в существенном превышении индекса его активности для сигналов одного из классов над другими.
Однако, как показывают эксперименты и как показано в работах [7,8], в ряде случаев классификация сигналов лишь на основе индексов активности СП не оказывается эффективной. Связано это с тем, что сигналы различных классов могут иметь одни и те же СП с примерно
одинаковыми индексами активности. Различие же между классами сигналов проявляется в последовательности смены этих статических паттернов.
Обобщим понятие индекса активности (1.7) на цепочку статических паттернов £Я Я = ,...,Sq ) длины г, которую далее будем называть динамическим паттерном (ДП) г-го порядка:
Т($в,...,Sq ) —
Р(ЕЯ1,...,Яг ) = P(Sq1,..., ^ ) = Т- г + 1Я , Я' G{1,..., Q}, ' = 1, Г , (18)
где Т^,...,Sq ) - количество вхождений цепочки (я1,...,Яг) в последовательность номеров СП ^Я^), t = 1, Т|, соответствующую рассматриваемому сигналу ), t = 1, Т|. Индекс активности
динамического паттерна фактически представляет собой относительную частоту встречаемости соответствующей последовательности переходов между статическими паттернами на рассматриваемом интервале времени.
С ростом порядка г число возможных динамических паттернов комбинаторно растёт, в то время как индексы активностей динамических паттернов, наблюдающихся в сигнале, стремятся к нулю - с увеличением длины цепочки вероятность её возникновения уменьшается. Для ДП порядка г = 2 индексы активности (1.8) образуют матрицу Р размерности Q*Q. Визуально ДП 2-го порядка могут быть представлены в виде ориентированного графа с взвешенными ребрами, каждая вершина которого соответствует статическому паттерну, а веса рёбер равны индексам соответствующих ДП 2-го порядка.
2.3.Алгоритм формирования вектора характерных признаков сигнала
Индексы активности динамических паттернов являются интегральными характеристиками временного ряда на рассматриваемом интервале времени. Однако использование индексов активности всех ДП (при г = 2 - всех элементов матрицы Р) в качестве характерных признаков сигнала для последующей классификации не представляется рациональным в силу нескольких
причин. Во-первых, некоторые ДП являются общими для сигналов всех классов, т.е. отражают некоторые общие свойства динамической системы, продуцирующей сигналы, и поэтому не представляют интереса с точки зрения классификации. Присутствие других ДП (с невысокими индексами активности) может являться индивидуальной особенностью конкретного сигнала или быть обусловлено действием различного рода случайных факторов и не может являться характеристикой класса сигналов. В связи с этим такие ДП также целесообразно исключить из дальнейшего рассмотрения. Во-вторых, общее число ДП г-го порядка может быть довольно большим, особенно при высоких порядках г, что усложняет процесс их классификации.
Введем понятие специфичности динамического паттерна £ по отношению к паре классов (k2), ^ е {1,...,К}, k2 е {1,...,К}, ^ * k2. Под специфичностью ДП £ будем понимать степень
успешности классификации значений индексов активности этого ДП на обучающих сигналах, относящихся к классам ^ и с помощью некоторого классификатора. В качестве показателя специфичности динамического паттерна £ по отношению к паре классов (k2 ), таким образом, выберем точность работы классификатора:
(£) = 2
N N
k1 е{1,...,К}, к2 е{1,...,К}, кх * k2, (1.9)
где и " - число правильно классифицированных примеров классов ^ и k2 соответственно по значению индекса активности ДП £, а и Nk - общее число примеров обучающей выборки соответствующих классов. Возможные значения показателя ¡^ к (£) лежат в отрезке [0; 1]. Значение ¡^^ (£) = 1 означает разделимость классов ^ и k2 с помощью выбранного классификатора.
Наиболее простым классификатором индексов активности р динамического паттерна £ на два класса ^ и ^ является линейный классификатор:
С[р(Е)Як" есл""(Е)£Р" или с[р(Е)] = {*" еслиР(Е)£(1.10)
[к2, в ином случае, [к^ в ином случае,
где значение разделяющего порога р0 выбирается таким образом, чтобы обеспечить максимальное значение критерия (1.9).
Отсортируем динамические паттерны, встречающиеся в обучающих сигналах классов к1 и к2, по значению показателя специфичности (1.9) в порядке убывания. Обозначим ££\ - ДП,
имеющий '-ый номер в таком отсортированном ряду, а ^ = Е^ - наиболее специфичный
динамический паттерн. Предполагая, что индекс ДП с максимальной специфичностью в наибольшей степени характеризует тип классифицируемого сигнала, будем использовать его в качестве характерного признака сигнала ), t = 1, Т| для его классификации на два класса: к1 и к2:
2к1кг = р (Ек1,к2). (1л1) Отметим, что для большей устойчивости классификации размерность вектора гкк может быть увеличена:
*М2 =(р(Ек1)к2),...,Р(Ек",)к2)) , (1.12)
где X - число используемых для классификации наиболее специфичных ДП. Далее будет рассматриваться вектор характерных признаков гкк вида (1.11).
Задача линейной классификации сигналов на К > 2 классов может быть сведена к построению ряда линейных классификаторов на два класса вида (1.10), каждый из которых осуществляет классификацию в пространстве соответствующих признаков гкк , к1 = 1, К , к2 = 1, К, к1 ^ к2. Таким образом, вектор признаков г для классификации сигнала на К классов будет со-
ставлен из всех векторов гкк , к1 = 1, К -1, к2 = к1 +1, К :
г = (^¡2,..., 2к-1,к ) . (1.13)
Размерность вектора z не превышает К*(К-1) / 2. Именно этот вектор предполагается использовать в качестве вектора характерных признаков сигнала ), t = 1, Т} для его классификации.
2.4. Алгоритм оптимизации специфичности динамических паттернов
Для некоторых классов к1 и к2, кх = 1, К , к2 = 1, К, кх * k2, выделенные наиболее специфичные динамические паттерны г-го порядка, образующие вектор гкк , могут иметь значения показателя специфичности (1.9), не позволяющие проводить классификацию с помощью выбранного классификатора с требуемой точностью. Причиной этому может быть неудачное разбиение пространства МХП V на статические паттерны. Индексы активности, а, следовательно, и специфичность ДП г-го порядка £, , = (5,,...,) зависят от формирующих его статических паттернов ,...,, т.е. областей в пространстве МХП V. Изменение вида этих областей будет отражаться на значении показателя специфичности I^к (£, , ) , к1 = 1, К , к2 = 1, К, к1 * к2.
Поставим задачу разбиения пространства Vна статические паттерны 51,...,таким образом,
чтобы максимизировать точность работы выбранного классификатора в пространстве признаков г, определяемом выражением (1.13). В связи с тем, что точность классификации на К классов связана с точностью попарной классификации каждой пары классов, введём следующий
показатель точности:
I = Шл (£и ). (1.14)
2
к =1
к2 =1,
Значение показателя I, равное 1, будет означать разделимость всех пар классов в пространстве характерных признаков г с помощью выбранного классификатора.
Задача поиска оптимальных статических паттернов состоит в разбиении пространства признаков МПХ Vна Q областей S1,...,Sq таким образом, чтобы максимизировать значение точности классификации (1.14): I ^ max ,
(^ -,SQ )
и = V, (1.15)
5,1 П 5,2 =0, ,1 = Ш ,2 = Ш ,1 * ,2;
Задачу (1.15) можно свести к задаче разбиения множества векторов МХП сигналов обучающей выборки, каждый из которых представляет собой точку в пространстве V, на Q групп. В связи с тем, что число точек в пространстве V может достигать десятков или даже сотен тысяч, перед решением оптимизационной задачи (1.15) сократим её размерность. Для этого проведём предварительную кластеризацию данных в пространстве V на N кластеров (Ы >> Q), используя стандартные методы кластеризации, например, метод к-средних или самообучающиеся карты Ко-хонена [5,9].
Обозначим 51,...,sN - кластеры, полученные в результате кластеризации данных в пространстве V, которые далее будем называть микросостояниями. Таким образом, задача (1.15) сводится к задаче разбиения N кластеров на Q групп (задаче формирования Q статических паттернов из N микросостояний) таким образом, чтобы максимизировать критерий (1.14). Для решения этой комбинаторной задачи могут быть применены различные эвристические методы. В данной работе для её решения предложено использовать генетические алгоритмы (ГА) [10].
На рис. 1 представлена последовательность операций, выполняемых при формировании вектора характерных признаков г для временного ряда ), t = 1, Т}. Эта последовательность фактически определяет алгоритм построения оператора ¥ в модели (1.5).
{хк,р (t)
хк,р (t), t = 1,Тк, р = 1,Рк, к = 1,К\
> г
Расчет векторов МХП
> р,р ^), t = 1, Тр;, р =
Выделение микросостояний
> к — % \ г
Формирование статических паттернов
> г {^^Л}
Формирование динамических паттернов
> к=1, к , к2=1, г
Выделение наиболее специфичных динамических паттернов
> к = 1, К, кг = 1, 1
Формирование вектора характерных признаков
1. спектральный анализ,
2. вейвлет-анализ
1. метод ^средних,
2. самообучающиеся карты Кохонена
1. генетический алгоритм,
2. эвристические методы комбинаторной оптимизации
о}
, р = к = 1К\
Рис. 1. Схема формирования вектора характерных признаков гк'р для многомерного временного ряда |хк,р ^), t = 1, Тр) \ из обучающей выборки, р = 1, Рк , к = 1, К . Справа указаны методы,
которые могут быть использованы при выполнении соответствующей операции.
\
3. Классификация в пространстве характерных признаков
После выбора пространства характерных признаков г, следующим шагом алгоритма классификации является построение собственно классификатора данных с [ г ] в этом пространстве. В
связи с тем, что при формировании вектора 2 уже был выбран и использован классификатор для разделения пар классов по значениям индексов ДП, на данном этапе для разделения К классов тип классификатора целесообразно сохранить. Так, при использовании линейного классификатора (1.10) в расчётах показателя (1.9) классификатор с [ г ] также целесообразно выбирать ли-
нейным.
В настоящей работе в качестве классификатора с [ г ] предложено использовать многослойный персептрон [5]. Число входов персептрона определяется размерностью вектора г, рассчитываемого по формуле (1.13), число выходов определяется числом классов. Обучающая выборка, таким образом, представляет собой совокупность векторов
{(гк,р ;ок,р ), р = к = 1К},
где гк'р - вектор характерных признаков для сигнала {хк,р^), t = 1,Тр} }, ок,р =(ок'р,..,<зкКР) -
биполярный вектор желаемых выходов персептрона, определяющий принадлежность входного вектора к-му классу:
окр =
1, если г = к, -1, иначе.
Математическая модель многослойного персептрона имеет вид [5]:
* = **, sг = /г (И),
И = Ж т-1,
5г = (1; sг), s0 = г, г = ,
где R - число слоёв многослойного персептрона, 2 - вектор входов, у - вектор выходов, ^ - вектор выходов '-го слоя нейронов, И - вектор потенциалов '-го слоя нейронов, РР1 - расширенная матрица синаптических коэффициентов т-го слоя нейронов, £ - активационная характеристика 1-го слоя нейронов.
В результате работы классификатора подаваемый на его вход пример будет отнесён к одному из К классов либо будет принято решение о невозможности классификации. Построение классификатора проводится на основе данных обучающей выборки, оценка качества его работы - на основе данных тестовой выборки. Для оценки качества классификации будем использовать показатель
1 к п
' -кЫ> (116>
где пк - число верно распознанных примеров к-ого класса при предъявлении Ык сигналов этого класса. Значение показателя (1.16) равное 1 означает правильную классификацию всех предъявленных на вход классификатора примеров.
4. Результаты экспериментальных исследований
4.1. Описание исходных данных
Экспериментальные исследования предложенных алгоритмов классификации многомерных нестационарных временных рядов проводились на сигналах электроэнцефалограмм (ЭЭГ), полученных в лаборатории высшей нервной деятельности человека Института высшей нервной деятельности и нейрофизиологии Российской Академии Наук. Каждая ЭЭГ соответствует процессу решения испытуемым задачи одного из К = 4 типов, связанных с различными типами мышления: планиметрическое мышление (тип 'G'), пространственное мышление (тип 'С'), составление предложений (тип и исключение лишнего слова из перечня (тип ' РР").
Испытуемый решил примерно по 70 задач каждого типа, потратив на каждую от четырех секунд до нескольких десятков секунд. Все сигналы ЭЭГ были оцифрованы с частотой дискретизации шд = 250 Гц и очищены от артефактов. Запись ЭЭГ велась по 22 каналам, из которых информативными являются L = 19. Все сигналы были разбиты на обучающую и тестовую выборки, объёмы которых соответственно равны Робуч = 48+55+42+53 = 198, Ртест = 21+24+19+23 = 87 (Р = Р1+ Р2+ Р3+ Р4).
4.2. Расчет векторов моментных характерных признаков
В связи с тем, что электрическая активность мозга имеет гармоническую природу [11], для формирования моментных характерных признаков многомерного сигнала ЭЭГ предложено использовать непрерывное вейвлет-преобразование (НВП) с материнской вейвлет-функцией Морле, имеющей аналитический вид [12]:
у^) = /2 ег 2л*.
Расчёт вейвлет-коэффициентов сигнала {х^), t = 1, Т}, производился по формуле
1 ^ _('ТТ *(т-ь
Ж (а, Ь) = У х(0 Г у* I — I dт
ыа г=0 ¿Т V а )
где АТ - интервал дискретности, а, Ь - параметры масштаба и временного сдвига соответственно, - материнская вейвлет-функция, оператор * означает комплексное сопряжение.
Для представления результатов НВП сигнала 1-го канала ЭЭГ в памяти компьютера используется матрица Ж1 = } вейвлет-коэффициентов, рассчитанных в узлах дискретной сетки на
плоскости (а, Ь): wl.. = Ж1 (аг, Ь.), I = 1, L .
г' ]-
На рис. 2 показан пример визуального представления матрицы квадратов вейвлет-коэффициентов (вейвлетограммы), рассчитанной для одного из каналов ЭЭГ. Каждый столбец матрицы отражает распределение энергии сигнала по частотам в фиксированный момент вре-
мени. Каждая строка отражает распределение энергии сигнала на фиксированной частоте по времени.
Рис. 2. Вейвлетограмма сигнала ЭЭГ в диапазоне частот от 5 Гц до 30 Гц; по оси абсцисс отложены значения временных сдвигов Ь, по оси ординат - масштабных сдвигов а, соответствующих частотам.
Поскольку сигналы ЭЭГ являются многомерными, вектор МХП ЭЭГ составлен из значений энергии сигнала каждого из каналов ЭЭГ в заранее заданных частотных диапазонах О^... Од в момент времени ^
V(t) = (V (t),...,V, (t)), V,(0 = (у„ (0,...,^ (t)), I = й,
Чт (t) = 2 1 = 1, Ь т = ^
Размерность М вектора МХП v(t) равна д*Ь, где д - число выбранных частотных полос, Ь -число каналов (Ь = 19). В работе исследовался диапазон частот от 5 до 30 Гц, что соответствует диапазонам альфа- и бета-ритма. В качестве элементов вектора МХП выбраны моментные значения энергии в частотных полосах шириной 0.5 Гц в рассматриваемом диапазоне. Таким образом, размерность вектора МХП v(t) равна М = 50*19 = 475.
4.3. Выделение микросостояний
Для каждого обучающего сигнала в каждый момент времени были рассчитаны векторы МХП jvk'p (t), t = 1, Tpk I, p = 1, Pk , к = 1, K . Общее число векторов в пространстве МХП оказалось
примерно равным 182000. Для сокращения их числа была проведена их кластеризация на N = 200 кластеров с использованием самообучающихся карт Кохонена с механизмом утомления CWTA (Consciousness Winner Takes All) [5]. В процессе самообучения уровень соседства экспоненциально менялся от Xmax = 100 до Xmin = 0,001, коэффициент обучения линейно менялся от Цтах = 1 до Цтп = 0,001. Самообучение продолжалось в течение 100 эпох.
В результате процесса кластеризации векторов МХП были сформированы N = 200 микросостояний, из которых далее формировались статические паттерны. На рис. 3 показано распределение числа точек в кластере по кластерам после завершения процедуры кластеризации.
Рис. 3. Распределение числа точек в кластере по N = 200 кластерам после завершения процедуры кластеризации.
4.4. Формирование статических паттернов
Формирование статических паттернов S1,...,SQ из микросостояний s1,...,sN предполагает решение оптимизационной задачи (1.15) разбиения N микросостояний на Q групп, при котором обеспечивается оптимальное значение критерия (1.14). Для решения этой оптимизационной задачи был использован генетический алгоритм.
В связи с тем, что каждое микросостояние должно входить ровно в один статический паттерн, сформируем хромосому следующим образом:
g = (&,...,gN), ^ е{1,...,Q}, 1 = , где gi - 1-ый ген хромосомы g, представляющий собой целое число в диапазоне от 1 до Q. Значение gi = q означает, что микросостояние si принадлежит статическому паттерну Sq. Каждая хромосома, таким образом, однозначно определяет разбиение множества микросостояний
{л1,..., SN
} на статические паттерны S1,...,SQ . Генетические операции над такими хромосомами
определяются аналогично генетическим операциям для бинарных генетических алгоритмов.
Отметим, что никакой априорной информации о числе статических паттернов в пространстве МХП нет, в связи с этим определение числа статических паттернов, при которых наблюдается приемлемая точность классификации, определялось экспериментально. На рис. 4 показаны графики зависимости значений приспособленности лучшей хромосомы в популяции и средней приспособленности в популяции от номера популяции (такта работы генетического алгоритма) для заданного числа статических паттернов Q, рассчитанные для сигналов обучающей и тестовой выборок. Отметим, что лучшие хромосомы выбирались из всех хромосом популяции, при этом лучшие хромосомы на обучающей и тестовой выборках в общем случае являются разными хромосомами.
Исследование поведения на тестовой выборке хромосомы, обладающей лучшей приспособленностью на обучающей выборке, показало, что практически во всех случаях её приспособленность на тестовой выборке не является лучшей, иными словами в популяции практически
всегда существует хромосома с лучшей приспособленностью на тестовой выборке, при этом не обладающей лучшей приспособленностью на обучающей выборке. На рис. 5 представлены графики зависимостей значений приспособленности на обучающей и тестовой выборках для двух хромосом - показывающей лучший результат на обучающей выборке и лучший результат на тестовой выборке при Q = 15. Подобная картина наблюдается и при других значениях Q.
В работе использованы следующие параметры генетического алгоритма. Длина хромосомы -N = 200, численность популяции - 40, алгоритм формирования родительского пула по методу рулетки, кроссинговер равномерный с вероятностью обмена генами 0.5, элитная стратегия селекции, вероятность мутации - 0.05.
0 05_I_I_I_I_I_I_I о_I_I_I_I_I_I_I
'о 50 100 150 200 250 300 150 0 50 100 150 200 250 300 350
рори1айоп рори1айоп
<* 0.4
<л
в) I
г- tr^in
k 1, test
Ам J TlfTlr №
г..........
0.5
train
test
i
population population
Рис. 4. Графики зависимости значений приспособленности лучшей хромосомы в популяции (слева) и средней приспособленности в популяции (справа) от номера популяции (такта работы генетического алгоритма) для заданного числа статических паттернов Q = 5 (а), 10 (б), 15 (в). Значения рассчитаны для сигналов обучающей (train) и тестовой (test) выборок.
Рис. 5. Графики зависимости значений приспособленности на обучающей (train) и тестовой (test) выборках для двух хромосом - хромосомы с лучшей приспособленностью на обучающей выборке (графики train 1 и test 1) и хромосомы с лучшей приспособленностью на тестовой выборке (графики train 2 и test 2). Число статических паттернов Q = 15.
Из графиков видно, что с ростом числа статических паттернов значение функции приспособленности лучшей хромосомы, полученной в результате работы генетического алгоритма, растёт. Отметим, что в процессе работы генетического алгоритма, несмотря на рост приспособленности лучшей хромосомы и средней приспособленности хромосом в популяции на обучающей выборке, на тестовой выборке эти значения практически перестают изменяться уже после 100200 итераций генетического алгоритма. Кроме того, значения функции приспособленности, рассчитанные по данным тестовой выборки, для хромосомы, имеющей лучшую приспособленность на обучающей выборке, также практически перестают изменяться, что говорит о неспособности эволюционного процесса к дальнейшему улучшению обобщающих свойств найденного решения. Значения критерия (1.14) после окончания работы генетического алгоритма на обучающей выборке при Q = 5, 10, 15 равны соответственно I = 0.37, 0.50, 0.62, рассчитанные для этих же хромосом значения приспособленности на тестовой выборке - I = 0.11, 0.11, 0.14 соответственно.
Для дальнейшего анализа была отобрана лучшая хромосома (Гбуч = 0.62, 1тест = 0.14). Эта хромосома определяет разбиение N = 200 микросостояний на Q = 15 статических паттернов. На рис. 6 приведено распределение числа микросостояний по найденным статическим паттернам, рассчитанное для найденного разбиения.
Рис. 6. Распределение числа микросостояний по статическим паттернам для лучшей хромосомы. По оси абсцисс отложен номер статического паттерна, по оси ординат - отсортированное относительное число микросостояний. Общее число микросостояний N = 200.
4.5. Выделение наиболее специфичных динамических паттернов
Как только сформированы статические паттерны S1,...,SQ, далее из них формируются динамические паттерны 2-го порядка Е = , S ), q1 е {1,..., Q}, q2 е {1,..., Q}. Для каждого из полученных динамических паттернов на основе сигналов обучающей выборки рассчитывается показатель специфичности (1.9) по отношению к классам ^ и k2, ^ е{1,...,K}, k2 е{1,...,K}, k1 Ф k2. В качестве классификатора двух классов использован линейный классификатор вида
(1.10). На рис. 7 приведены графики отсортированных значений показателя специфичности для некоторых пар классов в зависимости от порядкового номера динамического паттерна 2-го порядка. Общее число ДП 2-го порядка Q*Q = 15*15 = 225.
dynamical pattern dynamical pattern dynamical pattern
а) б) в)
Рис. Графики отсортированных значений показателя специфичности динамических паттернов 2-го порядка для классов 'C', 'S' (а), 'C', G (б), 'G', W (в), рассчитанные для сигналов обучающей выборки. По оси абсцисс отложен порядковый номер динамического паттерна. Число соответствующих статических паттернов Q = 15.
В табл. 1 приведены наиболее специфичные ДП Е^ к для каждой пары классов (к1, к2), рассчитанные по данным обучающей и тестовой выборок, и указаны значения показателя специфичности для них.
Таблица 1
Наиболее специфичные ДП 2-го порядка и значения показателей специфичности
С5 CG СЖ SG GW
ДП (9,9) (9,3) (9,9) (12,12) (2,11) (15,12)
обуч. 0.9792 0.8782 0.9527 0.9589 0.8735 0.8913
тест. 0.8728 0.7473 0.8928 0.6751 0.5958 0.7179
Из табл. 1 видно, что наиболее хорошо разделимы классы 'С', '£', наименее хорошо - классы 5', Ж. Этот результат согласуется с нейрофизиологическими представлениями - когнитивные задачи на составление предложений (тип '£') и исключение лишнего слова из перечня (тип оба являются вербальными и задействуют схожую мыслительную активность, в то время как при решении задач на пространственное мышление (тип 'С') и на составление предложений (тип '£') активизируются принципиально различные участки коры головного мозга.
4.6. Формирование вектора характерных признаков и классификация
Вектор характерных признаков г в соответствии с (1.13) формируется из значений индексов активности наиболее специфичных ДП, приведённых в табл. 1. В связи с тем, что для классов 'С', '5' и 'С', 'Ж оказался один и тот же наиболее специфичный ДП, размерность вектора 2 равна 5. По построению каждый элемент вектора 2 является оптимальным с точки зрения попарного линейного разделения сигналов на соответствующую пару классов. На рис. 8 приведено рас-
пределение проекций вектора 2, рассчитанных для сигналов обучающей выборки, на прямую, соответствующую выбранной паре классов. На рис. 9 приведено аналогичное распределение проекций вектора 2 на плоскость.
Рис. 8. Распределение проекций вектора г, рассчитанных для сигналов обучающей выборки, на прямую, соответствующую паре классов 'С', '5'. Фактически по оси абсцисс отложены значения индексов ДП (9,9) для сигналов класса 'С' (помечены 'о') и сигналов класса '5' (помечены 'х').
Рис. 9. Распределение проекций вектора г, рассчитанных для сигналов обучающей выборки, на плоскость индексов ДП (9,9) - ось абсцисс, и (12,12) - ось ординат. ДП (9,9) является наибо-http://technomag.edu.ru/doc/151934.html Страница 27
лее специфичным для классов 'С', ДП (12,12) - для классов 'О'. Символами 'о', 'х', '•' помечены индексы классов 'О', 'С' соответственно.
Из рис. 9 видно, что классы 'С', \У хорошо разделимы вертикальной прямой (по значениям индексов ДП (9,9)), а классы 'О', \У хорошо разделимы горизонтальной прямой (по значениям индексов ДП (12,12)). Однако классы 'С', 'О' в указанной плоскости разделимы плохо - наиболее специфичный ДП для этих классов в соответствии с табл. 1 - (9,3).
Векторы характерных признаков гк'р, рассчитанные для всех сигналов обучающей выборки, используются для обучения многослойного персептрона, предназначенного для их классификации. Объёмы обучающей и тестовой выборок соответственно равны р°буч = 198, ртест = 87. Обучение продолжалось по методу Левенберга-Маркардта в течение 2000 эпох.
Число слоёв и нейронов в слоях многослойного персептрона подбиралось экспериментально. Наилучшие результаты по точности работы персептрона на тестовой выборке наблюдались при числе слоёв д = 3 и распределении нейронов по слоям: N1 = 8, N2 = 5, N3 = 4. Активационные характеристики скрытых слоёв - сигмоиды, выходного слоя - линейная. Достигнутые значения показателя точности (1.16) на обучающей и тестовой выборках соответственно Рбуч = 0.81, 1тест = 0.56. Доля числа правильно распознанных примеров каждого класса (значения п^к в формуле (1.16)) обученным нейросетевым классификатором на обучающей и тестовой выборках приведена в табл. 2.
Таблица 2
Доля числа правильно распознанных примеров каждого класса обученным нейросетевым классификатором на обучающей и тестовой выборках.
Класс 'С О Я 'Г
обуч. 0.9792 0.6182 0.8810 0.7547
тест. 0.7619 0.2917 0.5263 0.6522
На рис. 10 приведено распределение значений выходов обученного нейросетевого классификатора при подаче на вход сигналов обучающей и тестовой выборок различных классов.
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 О
I
I
Ж
с в э и
а)
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
№
б)
Рис. 10. Распределение значений выходов обученного нейросетевого классификатора при подаче на вход сигналов обучающей (а) и тестовой (б) выборок различных классов. Каждый рисунок содержит 4 гистограммы (по числу классов), каждая из которых содержит 4 столбика, соответствующие классам 'С', 'G', 5', Ж. Высота 1-го столбца к-го класса представляет собой относительное число сигналов к-го класса, отнесённых классификатором к 1-му классу (значения пМк).
Из табл. 2 видно, что наиболее точно распознаются сигналы класса 'С', что говорит о высокой специфичности сигналов ЭЭГ этого класса. Наименьшая точность распознавания наблюдается для сигналов класса 'G', при этом, как видно из рис. 10, ошибки распознавания сигналов класса 'G' на тестовой выборке практически равномерно распределились "в пользу" всех остальных классов, что говорит о невысокой специфичности сигналов ЭЭГ данного класса по отношению ко всем остальным классам. Ошибки распознавания сигналов класса '5' в основном приходятся
"в пользу" класса }Ж, что соответствует нейрофизиологическим представлениям о деятельности мозга при решении когнитивных задач этих типов.
Заключение
Результаты применения предложенного алгоритма для выделения характерных признаков многомерных нестационарных сигналов ЭЭГ и их классификации показали его работоспособность на имеющихся экспериментальных данных. Точность нейросетевого классификатора в построенном пространстве признаков составила Гбуч = 0.81 на обучающей выборке и 1тест = 0.56 на тестовой. Установлено, что наиболее специфичными являются сигналы ЭЭГ класса 'С' (соответствующие ментальной активности при решении задач на пространственное мышление), точность распознавания сигналов этого класса достигает 98% на обучающей выборке и 76% на тестовой. Наименее специфичными являются сигналы класса 'О' (соответствующие задачам на планиметрическое мышление), которые вносят наибольший вклад в результирующую ошибку работы классификатора.
Несмотря на достигнутый результат, предложенный алгоритм требует проведения ряда исследований, которые можно объединить в следующие группы:
- исследование различных пространств моментных характерных признаков;
- исследование процедуры кластеризации векторов МХП и формирования множества микросостояний;
- исследование процедуры формирования статических паттернов из микросостояний;
- исследование возможности и результатов применения динамических паттернов высоких порядков;
- исследование различных вариантов формирования вектора характерных признаков на основе индексов активности динамических паттернов;
- исследование процедуры обучения многослойного персептрона для классификации выделенных векторов характерных признаков.
Результатом этих исследований должны стать рекомендации по выбору числа микросостояний и числа статических паттернов, определение порядка вводимых в рассмотрение динамических паттернов, достаточного для обеспечения требуемой точности классификации, а также ряд других выводов по применимости и эффективности предложенного алгоритма формирования характерных признаков для классификации многомерных нестационарных временных рядов.
Работа выполнена при финансовой поддержке Совета по Грантам Президента РФ для поддержки молодых российских учёных (Грант МК-4245.2009.8), а также в рамках гранта РФФИ 10-08-00816 и ФЦП "Научные и научно-педагогические кадры инновационной России" на 20092013 годы.
Литература
1. M.B. Priestley. Nonlinear and Nonstationary Time Series Analysis. Academic Press, 1988.
2. Nonlinear and Nonstationary Signal Processing. Edited by W. J. Fitzgerald. University of Cambridge. 2001.
3. R. Dahlhaus. Fitting time series models to nonstationary processes // The Annals of Statistics. 1997. Vol. 25. № 1. Pp. 1-37.
4. G. McLachlan. Discriminant Analysis and Statistical Pattern Recognition. New York: John Wi-ley&Sons, 1992.
5. Осовский С. Нейронные сети для обработки информации. - М.: Финансы и статистика, 2004.
6. Wolpaw J.R., Birbaumer N., McFarland D.J., Pfurtscheller G., Vaughan T.M. Brain-computer interfaces for communication and control // Clinical Neurophysiology, 2002, V.113, P.767-791.
7. В.И. Скругин, А.О. Роик, Р.А. Наумов, Трофимов А.Г. Алгоритм классификации сигналов ЭЭГ на основе анализа в частотно-временной области // Сборник научных трудов XII Всероссийской научно-технической конференции "Нейроинформатика-2010". - М.: НИЯУ МИФИ, 2010. Т.1. С.266-276.
8. В.И. Скругин, А.О. Роик, Р.А. Наумов, Трофимов А.Г. Алгоритм выделения пространственно-частотных паттернов в структуре сигнала ЭЭГ // Материалы избранных научных трудов по теме: «Актуальные вопросы нейробиологии, нейроинформатики и когнитивных исследований». - М.: НИЯУ МИФИ, 2010. С.213-223.
9. С.А. Айвазян, Мхитарян В.С. Прикладная статистика и основы эконометрики: Учебник для вузов. - М.: ЮНИТИ, 1998.
10. Емельянов В. В., Курейчик В. В., Курейчик В. М. Теория и практика эволюционного моделирования. - М: Физматлит, 2003.
11. Зенков Л.Р. Клиническая электроэнцефалография (с элементами эпилептологии). Руководство для врачей. - М.: МЕДпресс-информ, 2004.
12. Витязев В.В. Вейвлет-анализ временных рядов: Учеб. пособие. - СПб.: Изд-во С.-Переб. ун-та, 2001.
