Авторегрессионный анализ в электроэнцефалографии Текст научной статьи по специальности «Математика»

АВТОРЕГРЕССИОННЫЙ АНАЛИЗ В ЭЛЕКТРОЭНЦЕФАЛОГРАФИИ И.О. Жаринов

Рассматривается применение для анализа и математической аппроксимации ЭЭГ-сигналов параметрических моделей авторегрессии-скользящего среднего (АРСС) - вычислительные методы, статистические аспекты, границы применимости.

Введение. Постановка задачи

Использование для анализа и математической аппроксимации ЭЭГ-сигналов параметрических моделей авторегрессии-скользящего среднего (АРСС) находится сейчас в стадии интенсивного развития: совершенствуются вычислительные методы, уточняются статистические аспекты, выясняются границы применимости.

АРСС-анализ базируется на предположении, что текущие значения сигнала имеют статистическую связь с его "предысторией". АРСС-модель ЭЭГ представляет значения дискретных отсчетов у(' )[и] /-го квазистационарного участка линейным

соотношением вида

„(')

„(')

У

(')

[п]+ 2

а

(') ,,(')

к

■ У

)[п - к ] = ■ х[п] + 2 ЪкР ■ х[п - к ] , '=1,2,.. .М,

к=1

к=1

где {х[п]} — порождающий процесс — последовательность независимых, нормально распределенных случайных величин с нулевым математическим ожиданием

и

единичной дисперсией. Коэффициенты {а['),а2'),...,а('()/) } и {>о'),Ь(),. .,Ь(()) }, а также

величины и д(') являются параметрами модели. Параметры р(-/ и определяют порядок модели авторегрессии (АР) и скользящего среднего (СС) соответственно, а

,(')

,(')

величина (р(') + д(')) определяет порядок АРСС-модели. Индекс (') означает описание фрагмента сигнала моделью '-го класса из М возможных взаимоальтернативных классов.

Адекватную параметрическую модель ЭЭГ удается получить, полагая

равными нулю, при этом уравнение определяет

коэффициенты

] иарашыи

К'".-.. ь% I

авторегрессионную модель

р(')

У

(')г

[п] + 2 ак° ■ У(0[п - к] = Ь^ ■ х[п] , 7=1,2,...М .

к=1

В последующем через параметры АРСС-модели могут быть выражены модель сигнала (см. рис.1,а), оценки частотных и корреляционных функций, определены отрезки квазистационарности сигнала, решены задачи сегментации и классификации.

В частности, при известных параметрах АРСС-модели оценка спектральной плотности (см. рис.2, б) реализации '-го класса имеет вид

£ (0( / ) = а 2 М-

„(О

Ь0° +2 Ьк° ■ ехР(- ]2фЫ)

к=1

р(')

1 + 2 а^ ■ ехр(- ]2п/кА0

к=1

, ¡=1,2,...М

где Аt - период дискретизации сигнала; а2 - дисперсия порождающего БШ {х[п]}; у -мнимая единица.

Уравнению авторегрессии отвечает линейный цифровой формирующий фильтр (ФФ), который входной последовательности {х[п]}0^ 1 ставит в соответствие выходную

2

2

последовательность {у(г)[и]}^ . Передаточная функция (ПФ) этого фильтра имеет следующий вид:

Ж (0( г) =

( (■) ф

IЬ?

.-к

Л /Г рО

к=0

1 + 1 а)У ■ г

к=1

(0 ■.-к

к

, 7=1,2,...М•

; 1 1

'/ / : 1 \ \ 1 1

А1 / ^ Д -+-V - г / \ \ 1 ~ V " "V ~ 1 | - - + — + - - 41 - -1----

уУ 1 /' V у 1 1 \ х\ . - .1 \ Л^-Ол - л - -+ + V 1 -----

0.2 --1--

(б)

Рис.1. Параметрические модели временных реализаций (а) и спектральных плотностей (б) стандартных ЭЭГ-сигналов: 1. Низкочастотная активность на ЭЭГ; 2. Замедленный вариант а-ритма; 3. Нормальный а-ритм;

4. Убыстренный вариант а-ритма; 5. Пароксизмальная активность на ЭЭГ.

Тогда, если гауссовский стационарный некоррелированный порождающий процесс {[и]} пропустить через ФФ, то на его выходе будет реализация сигнала 7-го класса, обладающая достаточным сходством с исходной физиологической кривой. Причем специалисты-электрофизиологи отмечают, что исходные и синтезированные с помощью ФФ реализации визуально практически не отличаются друг от друга.

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

Условие стационарности определяется условиями устойчивости ФФ. Фильтр устойчив, если полюса его ПФ находятся внутри круга единичного радиуса на

Z-плоскости. Иными словами, характеристического уравнения

корни

,0)

к=1,2,..., р

7 = 1,2,...М

-р(° + £4) • гр(1)-к = 0 , '=1,2,...,М

должны удовлетворять неравенству

(к)

/о < г1

к=1

< 1, к=1,2,..., р(1) 1=1,2,.. ,,М

Нестационарный характер модели ЭЭГ может учитываться зависимостью коэффициентов {4°}, к=1,2,..., р(1 и {>к')}, к=0,1,..., д(1у), от времени: {г^И} и

Ък' )[п]}.

Условие обратимости формулируется следующим образом: дискретный сигнал обратим только тогда, когда корни г^, к=1,2,..., 1=1,2,...М характеристического уравнения

1 + + ^^-2 + - + ^= 0 ■ —

удовлетворяют неравенству

(1)

< 1, к=1,2,...,д(1), 1=1,2,.. .М

Диапазоны изменения значений весовых коэффициентов АРСС-моделей малых порядков, соответствующие условию устойчивости и обратимости, приведены в табл. 1.

Вид модели Диапазоны изменения коэффициентов

АР-модель (р(1) = 1) -1< а[1) <1

АР-модель (р(1) = 2) 1+а^ + а2 )>0; 1- а?) + а21)>0; -1< а2°<1

АРСС-модель (р(1) = 2), (д(1) = 1) 1+а}') + а2°>0; 1-а^ + а£ )>0; -1<а?)<1; -1<Ь((1)/)<1

Таблица 1. Диапазоны изменения значений коэффициентов моделей ЭЭГ Задача формирования параметрического описания ЭЭГ

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

Рис. 2. Диаграмма параметрического синтеза моделей ЭЭГ-сигналов

Концептуальный выбор класса моделей. Выбор класса параметрических моделей

требует некоторых предварительных сведений о возможной форме функции спектральной плотности (см. рис.1) для наблюдаемых реализаций аппроксимируемого сигнала. Если необходимо сформировать модели сигналов с острыми спектральными пиками, но без глубоких впадин (нулей), то наиболее подходящей является АР-модель. СС-модель удовлетворительно аппроксимирует сигнал, если необходимы спектры с глубокими нулями, но без острых пиков. Что же касается АРСС-модели, то она может применяться в обоих этих предельных случаях. В тех случаях, когда одинаково пригодна любая из трех моделей, следует использовать ту из них, которая имеет наименьшее число параметров. Этот принцип экономии связан с тем, что получить оценки с хорошими статистическими свойствами удается тогда, когда число оцениваемых параметров минимально.

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

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

Оценивание параметров модели. Традиционно в задачах оценивания параметров рассматривается ряд свойств оценок, важнейшими из которых являются: состоятельность, несмещенность, достаточность и эффективность.

Подтверждение качества модели. Существует, по крайней мере, два различных подхода к решению задачи проверки адекватности модели. В первом подходе проверяется адекватность допущений, положенных в основу модели. Для классов параметрических АР- и АРСС-моделей во многих случаях единственное существенное допущение состоит в том, что порождающий процесс |х[и]} — это последовательность независимых случайных величин с нулевым средним. Используя выбранную модель и доступные наблюдения, можно получить оценки так называемых остатков или ошибок аппроксимации, и проверить предположение об их независимости с помощью теории проверки статистических гипотез.

Сущность второго подхода состоит в следующем. Чтобы удостовериться в адекватности модели, аппроксимирующей экспериментальные данные, необходимо непосредственно сравнить характеристики ее выходного сигнала, такие как коррелограммы, спектральные плотности, моментные характеристики и т. п. с соответствующими характеристиками исходных эмпирических данных. Осуществить такое сравнение можно аналитически или посредством математического моделирования. Модель принимается, если расхождение между двумя множествами характеристик не превышает одного-двух стандартных отклонений.

Оценивание параметров моделей

Оценивание параметров {акг)} и {Ьк(г)} линейного уравнения регрессии по имеющимся эмпирически полученным данным представляет собой одну из основных

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

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

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

Наибольшее распространение для анализа и статистической аппроксимации классов сигналов электроэнцефалограмм получили параметрические не АРСС-, а АР-модели. Связано это не только с характерной колоколообразной формой функции СПМ отдельных классов ЭЭГ, но и с тем, что сама процедура оценивания АРСС-параметров по экспериментальных данным достаточно трудоемка и требует решения систем нелинейных уравнений высоких порядков с одновременным (или последовательным) оцениванием АР- и СС-параметров АРСС-модели. Для решения таких уравнений применяются сложные итеративные численные методы оптимизации, основанные на использовании оценок максимального правдоподобия, или субоптимальные методы, которые, как правило, базируются на критерии наименьших квадратов.

Авторегрессионные модели ЭЭГ формируются в основном методом решения систем линейных уравнений и позволяют в подавляющем большинстве технических приложений проводить исследования (сегментация ЭЭГ, классификация ЭЭГ и т.д.) с достаточной для практики точностью. Существует достаточно большое число различных методов, позволяющих оценивать параметры АР-моделей по экспериментальным данным. Широкое распространение получил известный метод Юла-Уолкера, основанный на решении одноименной системы уравнений я (1)[0] я (1)[-1] ••• я (1)[-р(1)] я (°[1] я (1)[0] ••• я (1)[-р(1) +1]

_я(1)[р(1)] я(1)[р(1) -1] ••• я(1)[0]

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

1 (1) № )2

X Я] = 0

(') а р(1) 0

, '=1,2,.,.,М ,

стационарности и обратимости сигналов. Оценки параметров модели получают рекуррентным образом по нижеследующим выражениям:

а

(') ,(7)

= а (■) - а (■) "„('■ )к р(7)+1 р(7) +

а

(' )

, к=1,2,..., р(') , = 1,2,...М

р(7)+1,к р(' ),к р(' )+1,р(7 )+1 р(7),р(7)+1-к

,С)

а

(')

р (') +1,р (') +1

Г 1 р(7)

р(') +1]-1 а

к=1

(') ■ Г(7) р(7),к Г

[ (') + 1 - к

р (')

1 -I а«и ■ г(7) [к]

к=1 р ,к

р(7),к

г (7)[ к ]=

Я (')[0]

где Я(' )[к ] - смещенная оценка порядка к для автокорреляционной функции выборки. При использовании несмещенных оценок АКФ автокорреляционная матрица Я может оказаться неположительно-определенной, и, следовательно, эквивалентный ФФ будет неустойчивым.

Нулевой компонент скользящего среднего рассчитывается с использованием уже известных оценок коэффициентов корреляции и параметров модели

Ь (') = и0

1_^^^

1 - а}') ■г (')[]-а2') ■ г(') [2]-... - а('1) ■ г(7) [р(7) ]

, = 1,2,.„М

Сформированная таким образом модель '-ого класса ЭЭГ представляет собой косвенное математическое описание реализации сигнала, заданное набором (вектором) вещественных параметров.

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

модели, так и их порядками р

(')