CC BY
25
БИОФИЗИКА И МЕДИЦИНСКАЯ ФИЗИКА
DOI: 10.5862/JPM.253.1 0 УДК 61:577.3; 004.415.2.031.43
A.A. Иванков, А.В. Кутузов, Л.А. Криворучко
Санкт-Петербургский политехнический университет Петра Великого
размерности аттракторов детерминированных слагаемых
аддитивной модели сигналов, регистрируемых в транскраниальной ультразвуковой доплерографии
Разработаны алгоритм и математическое обеспечение для исследования аттракторов гармоник, формирующих пульсовую волну в сигналах, которые регистрируются в ходе доплерографического обследования внутримозгового кровообращения в организме человека. Приведены оценки корреляционной размерности этих аттракторов. Оценки получены в результате анализа реальных клинических данных для двух различных состояний обследуемых: полного покоя и в процессе тестирования на гипокапнию.
МОДЕЛЬ ДИНАМИЧЕСКОГО ХАОСА, КОРРЕЛЯЦИОННАЯ РАЗМЕРНОСТЬ АТТРАКТОРА, ЦЕРЕБРАЛЬНАЯ ГЕМОДИНАМИКА, АЛГОРИТМ.
Введение
Транскраниальная ультразвуковая до-плерография — это неинвазивная методика, с помощью которой изучают состояние систем регуляции внутримозгового кровообращения в организме человека. Суть ее состоит в том, что с помощью импульсных ультразвуковых сигналов исследователь имеет возможность регистрировать линейную скорость кровотока (ЛСК) в магистральных артериях мозга. Как правило, ЛСК регистрируют одновременно с динамикой системного артериального давления (АД). Затем эти данные используют для того, чтобы оценить состояние системы регуляции внутримозгового кровообращения.
Первым шагом при регистрации какого-либо физиологического показателя является попытка найти набор информативных признаков — вектор, элементы которого определяют положение системы в ее фазовом пространстве. такое описание физиологического состояния организма позволяет впоследствии сделать следующий шаг
на пути от категориальных оценок к количественным — построить математические модели, более сложные по своей структуре. Однако они более содержательны и легко интерпретируемы специалистами предметной области. Например, физиологи и кардиологи, изучая механизм иннервации сердечной мышцы, перешли, наконец, от методов последовательной бинарной классификации электрокардиограммы (ЭКГ) на основе ряда характерных паттернов (например, зубцов Р, О) или их сочетаний (РО5) к анализу ЭКГ методами теории динамических систем и стохастическими методами теории случайных функций [1].
Подобным образом дело обстоит и с изучением системы кровообращения методами доплерографии. В тех случая, когда обследуются пациенты с серьезными патологиями, опытные специалисты способны уже при регистрации сигналов поставить по их виду предварительный диагноз. Более тщательное обследование требует и более тонких подходов, а значит, и более сложных математических моделей исследу-
емых систем. В современной кардиологии размерности аттракторов, оценки которых получают в ходе анализа ЭКГ, уже давно пытаются интерпретировать как скалярную характеристику состояния системы иннервации миокарда. Ряд функциональных и органических патологий мозга удалось диагностировать только на основе оценок корреляционных интегралов и корреляционной размерности аттракторов, построенных в ходе анализа электроэнцефалограмм.
В связи с изложенным выше, нам хотелось бы подчеркнуть, что в настоящей работе мы не предлагаем конкретной методики, которую можно было бы непосредственно предложить практикующим клиницистам. В данном случае мы касаемся скорее вопросов методологии.
Дело в том, что в ходе исследований квазипериодических сигналов вопросы выбора адекватного метода их анализа тесно связаны с выбором априорной модели для самого сигнала (либо процесса). Если ограничиться наиболее популярными и (как полагают применяющие их исследователи) универсальными решениями, то, прежде всего, это методы спектрального анализа. Оценки спектра получают, анализируя как стационарные, так и нестационарные случайные процессы; к ним прибегают и в ходе изучения динамических систем. Однако такие оценки не позволяют найти ответ на один из самых существенных вопросов: имеют ли эти квазипериодические сигналы случайный характер или здесь наблюдается хаотическая динамика системы со сравнительно небольшим числом степеней свободы.
В своем исследовании мы попытались отыскать ответы именно на такие вопросы и руководствовались соображениями, что именно сердечные сокращения определяют, в конечном счете, основные изменения квазипериодических составляющих ЛСК в сосудах, а значит и в артериях, которые снабжают мозг кровью. следовательно, изучение тех периодических составляющих, гармоник в спектре регистрируемой ЛСК, частоты которых vd соизмеримы с частотой сердечных сокращений или кратны ей, представляет интерес для широкого
круга специалистов. Еще раз заметим, что эти гармоники вносят наибольший вклад в мощность регистрируемых сигналов и наиболее существенные изменения в регистрируемую ЛСК.
следующее предположение, на котором построен наш алгоритм анализа ЛСК, состоит в том, что вышеупомянутые периодические составляющие с частотами vd — это проекции траекторий динамической системы и указанные траектории демонстрируют свойства детерминированного хаоса. Условия проведения доплерографического обследования таковы, что обследуемый организм находится в состоянии практически полного покоя. Частота его сердечных сокращений изменяется незначительно, по крайней мере, в течение нескольких десятков секунд. В то же время иерархически организованная система регуляции кровообращения не может не проявлять себя на частотах vd. Эту гипотезу мы также попытались проверить в ходе нашего исследования.
В настоящее время одной из востребованных методик анализа сигналов, регистрируемых в ходе исследований медико-биологических систем, является построение оценки корреляционной размерности Х>2 аттракторов этих систем в их фазовых пространствах [2]. Указанная скалярная характеристика уже сама по себе может служить критерием в решающих правилах, построенных для идентификации различных физиологических состояний исследуемой системы. По этой причине в настоящей работе мы сделали акцент на алгоритме и математическом обеспечении для оценки корреляционной размерности аттракторов гармоник (их спектральной мощности), которые вносят наибольший вклад в мощность регистрируемых сигналов.
Основное внимание уделено анализу оценок плотности спектральной мощности сигналов, а не самих траекторий сигналов по двум причинам.
Во-первых, аналоговый сигнал, поступающий на вход аналого-цифрового преобразователя (АЦП), не обеспечивает достаточно широкого динамического диапазона ординат на выходе АЦП. Это свойство не
самого аналогового сигнала, а скорее процедуры калибровки аппаратного обеспечения непосредственно перед началом сеанса. Аппаратура калибруется таким образом, чтобы предвосхитить нестационарные эффекты, обнаруживающие себя в ходе регистрации сигнала. Другими словами, перед началом сеанса прибор калибруется так, чтобы будущие (в ходе сеанса) существен-
ные изменения среднего уровня сигнала не приводили бы к выходу сигнала за пределы заданного диапазона напряжений. Как следствие, в ходе обработки исходного сигнала мы получали многочисленные артефакты; для большинства значений параметров алгоритма не удавалось получить оценки корреляционной размерности.
Во-вторых, обработка сигнала в частот-
а)
350
ззо
310
290
* 270
Е 250
230
210
190
ч . /'ч Ч. ......*■•■-<■
>ч/ 5
----- / ~~ Ч ч N
N 2 \\ / ~ \ /--Ч 4
1 " \ \ \ ч \ ч ч N
\ \ \ "" \ \ ч ^
61
81
101
121 141 Ите, б
161
181
201
Рис. 1. Оценки первой (а) и второй (Ь) моментных функций сигналов ЛСК, зарегистрированных левым (ВРУ/, нечетные номера кривых) и правым (ВБУг, четные номера) доплеровскими датчиками в ходе обследований здорового добровольца; даны результаты первого сеанса (1 — 4) и второго (5 — 8), проведенного через полгода
ной области существенно проще. Алгоритмы фильтрации и устранения высокочастотного шума в частотной области хоть в какой-то мере обеспечивают робастность финальных оценок.
Постановка задачи
Прежде всего, следует отметить, что все регистрируемые сигналы представляют собой траектории нестационарных процессов. Эти процессы не стационарны в отношении как первой моментной функции (математического ожидания), так и меры рассеяния ординат. На рис. 1, а и 2, а приведены оценки первой моментной функции для типичных сигналов, зарегистрированных в ходе обследования здорового добровольца и пациента с патологиями внутримозгового кровообращения. На рис. 1, Ь и 2, Ь при-
ведены оценки второй моментной функции для этих же сигналов.
Отсылаем читателя к нашей предыдущей работе [3], где эти свойства оценок ЛСК и АД обсуждаются более подробно.
Мы рассматриваем исследуемые процессы как композицию случайных и детерминированных составляющих. Разделение слагаемых, если рассуждать об этой композиции как об аддитивной модели сигнала, — это самостоятельная задача. Изучаемый процесс имеет квазипериодический характер, причем периодическая составляющая, по всей видимости, имеет детерминированный характер. Случайные составляющие мы может отнести к классу нестационарных и/или квазистационарных случайных процессов. С точки зрения программно-аппаратной реализации алгоритма регистрации, сигнал —
а)
290
230
41
П-1-1-1-1-1-1-1-1-г
81 121 161 201 "Пте, в
241
1800
41 81 121 161 201
Пте, в
241
Рис. 2. Такие же оценки функций сигналов ЛСК (см. рис. 1), полученные в ходе обследований пациента с патологиями внутримозгового кровообращения. Представлены результаты одного сеанса
это процесс Xind(t) с дискретным аргументом t е Z+ (Z+ = 0, 1, 2, ...), или, другими словами, это временной ряд. Верхний индекс ind е {BFV/,BFVr,BP}. Семантика элементов этого множества, как и в предыдущих наших работах, следующая [3]:
BFV/, BFVr — линейные скорости кровотока, зарегистрированные левым (left) и правым (right) доплеровскими датчиками; BP — артериальное давление, зарегистрированное с помощью датчиков CNAP (Continuous Noninvasive Arterial Pressure), датчики фиксируются на левой руке (A. Digitalis).
В последующем изложении аргумент t фигурирует как индекс ординаты процесса. В тех случаях когда необходимо воспользоваться аргументом, имеющим размерность, т. е. определить в секундах моменты времени, которым соответствуют конкретные ординаты процесса Xind(t), достаточно разделить индексы этих ординат на частоту дискретизации АЦП.
В исследованиях, результаты которых приводятся в настоящей работе, восьми-канальный АЦП обеспечивал 125 ординат в секунду по каждому каналу (напомним, что в ходе обследования регистрировались как минимум две оценки ЛСК и динамика АД).
В качестве исходной примем следующую аддитивную модель плотности спектральной мощности (или периодограммы) регистрируемого сигнала Xind(t) на его участке
„ind
длины T
ip =IP (0)+1 ip (v)+ £ip (), (i)
где
I ( V ) =
i
N +1
£ X(t)e"
ipd (v) = T
£ Xind (t )eivt
tе[((d ),((d ))-1]
— оценка плотности спектральной мощности для гармоники частоты V; Т — длина участка, на котором построена оценка; vd, V я — частоты детерминированных и стохастических гармоник соответственно;
Jind '
,(0) -
вой моментной функции сигнала на этом рП-м участке длины Т.
Длина Т выбиралась таким образом, чтобы относительные изменения второй моментной функции оценки плотности спектральной мощности гармоник с частотами Vd на таких участках не превышали нескольких процентов (т. е. чтобы первые две значащие цифры оценки дисперсии были инвариантны).
Напомним (см. [3]), что множество {рП}, построенное нами таким образом, удовлетворяет следующему условию:
р1 : иР -Х (?), Р -
1
-{Хп(/):/ е[(?)Ар'1) + Т]},
где t(p'Jnd) — смещение участка по отношению к началу траектории Хш(().
В данном исследовании мы начали с построения этих участков как пересекающихся подмножеств Хш(() последовательным сдвигом на 125 точек (1 с) вдоль траектории Х'"\/) окна размером в Т = 1000 точек (8 с). Затем мы повторили оценивание, выполнив построение пересекающихся подмножеств Х'"\/) последовательным сдвигом на 125 точек (1 с) вдоль траектории Хш(() окна размером в Т = 4000 точек (32 с).
Для построения состоятельных и асимптотически несмещенных оценок (1) в нашем программном обеспечении реализован метод непараметрического усреднения (метод Уэлча).
Стохастические гармоники мы рассматриваем как физически осуществимый белый шум, который безусловно несет какую-то информацию об исследуемой системе, но ее интерпретация в настоящее время не представляется возможной. Поэтому для выделения из периодограммы детерминированных слагаемых, соответствующих гармоникам с частотами vd, потребовалось скорректировать зарегистрированный сигнал, а затем построить оценку периодограммы и отфильтровать стохастические слагаемые, соответствующие частотам vs.
Далее мы будем опускать и верхний
ш \
оценка плотности спектральной (ind е {BFV/,BFVr,BP}), и нижний ((d)
мощности, соответствующая оценке пер-
индексы, если в нашем изложении модель
2
t =0
2
распространяется на каждый из временных рядов и на каждый их участок.
корректировка исходного сигнала и фильтрация стохастических слагаемых
Алгоритм устранения из сигнала детерминированной составляющей, обеспечивающий в результате исключение слагаемого 1(0) Т = 1000 из модели (1), был изложен нами в статье [3].
Здесь мы более подробно остановимся лишь на алгоритме фильтрации стохастических слагаемых. Они составляют большую часть спектра скорректированных сигналов. Фильтрация этих гармоник осуществляется нашим программным обеспечением с помощью решающего правила, в основу которого положена проверка статистической гипотезы о согласии с равномерным распределением фиксированного подмножества оценки плотности спектральной мощности сигнала (траектории Хш(()). Мы упорядочиваем множество
{Iу)}, V е {1,[N /2]}
(квадратные скобки выражают операцию получения целой части своего аргумента) и получаем выборку
{I(п)}, п е {1,^ /2]}, I(п) > I(и+1),
сохраняя биекцию между индексами п и частотами V.
Затем из выборки {I(п) последовательно, начиная с п = 1, исключаем ее элементы и проверяем сложную статистическую гипотезу (Н0), о согласии оставшейся части выборки {I } с равномерным распределением. В основе решающего правила для проверки Н0 лежит вычисление статистики Морана [4].
Гипотеза проверяется при уровне значимости 0,05. Табл. 1 иллюстрирует процесс фильтрации стохастических слагаемых согласно этому алгоритму.
В каждой строке таблицы представлены первые 16 оценок плотности спектральной мощности, полученные нами для более чем
ВБУ/
двух десятков смежных участков р} траектории сигнала, зарегистрированного с помощью левого доплеровского датчика в ходе обследования здорового добровольца.
В ячейках таблицы указаны частоты гармоник в герцах. Отношение порядка элементов в пределах каждой строки определяет вклад мощности этих гармоник в мощность сигнала в порядке убывания, т. е. гармоники расположены в каждой строке в порядке убывания оценок их мощности. Серым фоном выделены частоты }, соответствующие детерминированным слагаемым; на белом фоне — часть стохастических слагаемых (подмножество {V^), которые попали в подмножество оценок, представленных в табл. 1.
Далее мы рассматриваем в качестве детерминированных слагаемых аддитивной модели (1) те элементы подмножества М с {I}, которые были последовательно удалены из выборки {Iчтобы нулевая гипотеза о согласии с равномерным распределением не отвергалась.
Воспользуемся биекцией между индексами п и следующими частотами vd :
{vd} = {0,750, 0,875, 1,0, 1,125} и
^{1,875, 2,0} и {2,750, 2,875, 3,0};
получаем соответствующие им оценки плотности спектральной мощности для каждого анализируемого участка р°Пс1 зарегистрированного сигнала.
Следует отметить, что для детерминированных слагаемых мы, как правило, получали в частотной области три непересекающихся подмножества:
{0,750, 0,875, 1,0, 1,125},
{1,875, 2,0, 2,125, 2,250}
{2,750, 2,875, 3,0, 3,125}.
Они соответствуют частоте сердечных сокращений (элементы первого подмножества) или почти кратны ей. Для каждого из этих подмножеств мы проводили анализ размерности аттракторов как отдельных гармоник, вносящих наибольший вклад в спектральную мощность сигнала на этом участке, так и различных сочетаний гармоник, принадлежащих этому подмножеству. кроме того, нами были построены оценки и для различных сочетаний гармоник (попарно) из различных подмножеств. При
таблица 1
Пример результатов фильтрации стохастических слагаемых
1,000 0,875 2,000 1,875 3,000 2,875 2,750 1,125 0,750 1,750 0,125 3,125
1,000 0,875 1,875 2,000 3,000 2,875 2,750 1,125 0,750 1,750 0,125 3,125
1,000 0,875 1,875 2,000 2,875 3,000 1,125 2,750 0,750 1,750 0,125 0,250
1,000 0,875 1,875 2,000 2,875 3,000 0,750 1,125 2,750 1,750 0,125 1,250
1,000 0,875 1,875 2,000 2,875 3,000 1,125 0,750 2,750 1,750 0,125 2,625
1,000 0,875 1,875 2,000 2,875 0,750 1,125 2,750 3,000 1,750 0,125 2,625
1,000 0,875 1,875 2,875 2,000 2,750 1,125 0,750 3,000 1,750 0,125 2,625
1,000 0,875 1,875 2,875 2,000 2,750 0,750 1,125 3,000 1,750 0,125 3,625
1,000 0,875 1,875 2,875 2,750 2,000 1,125 0,750 3,000 1,750 0,125 3,625
1,000 0,875 1,875 2,875 2,750 2,000 0,750 1,125 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,750 2,875 1,125 2,000 0,750 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,750 2,875 0,750 1,125 2,000 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,750 2,875 1,125 0,750 2,000 1,750 3,000 0,125 3,625
0,875 1,000 1,875 2,750 2,875 0,750 1,125 2,000 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,750 2,875 1,125 0,750 2,000 3,000 1,750 0,125 3,625
1,000 0,875 1,875 2,750 2,875 0,750 1,125 2,000 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,875 2,750 1,125 0,750 2,000 3,000 1,750 0,125 3,625
1,000 0,875 1,875 2,875 2,750 0,750 1,125 2,000 3,000 1,750 3,625 0,125
1,000 0,875 1,875 2,875 2,750 0,750 1,125 2000 3,000 1,750 3,625 0,125
1,000 0,875 1,875 2,875 2,750 0,750 2,000 1,125 3,000 1,750 3,625 0,125
1,000 0,875 1,875 2,875 2,750 2,000 1,125 0,750 3,000 1,750 0,125 3,625
1,000 0,875 1,875 2,875 2,750 1,125 0,750 2,000 1,750 3,000 0,125 3,625
1,000 0,875 1,875 2,875 2,750 1,125 2,000 0,750 3,000 1,750 3,625 0,125
Примечания. 1. В ячейках таблицы указаны частоты гармоник в герцах.
2. Сигналы зарегистрированы с помощью левого доплеровского датчика (ВРУ/), в ходе обследования здорового волонтера.
этом мы преследовали единственную цель — обнаружить значимые различия оценок размерности аттракторов детерминированных слагаемых (1) для двух различных физиологических состояний.
Алгоритм оценивания корреляционной размерности детерминированных слагаемых
Обобщенная размерность аттрактора системы в ее фазовом пространстве определяется по Реньи [5] как
1
• Иш-
1п £ р! (8)
Б = ____
! ! -1 8 >0 1п 8
где р1 (8) — вероятность попадания ордина-
ты траектории системы в г-ю ячейку (кубик с длиной ребра, равной 8) ее фазового пространства; N — число таких ячеек как функция 8; ! — параметр, определяющий, какая конкретно оценка размерности нам требуется.
Очевидная связь такого определения со случайной мерой, заданной на множестве ординат процесса в его фазовом пространстве, позволяет достаточно просто пояснить смысл такой оценки при ! = 2. Это корреляционная размерность, которая характеризует вероятностную связь между двумя ординатами траектории процесса.
На практике для построения оценки корреляционной размерности используют
следующее определение:
ln C (s) ln s
D2 = lim:
2 s^0
где C( s) — корреляционный интеграл вида
Z H(s — ||x,. - xj\|)
C (s) = j-2-.
Ns2
(2)
причем {х.}, V/ е {1,2, ...,N} — реконструированный или восстановленный аттрактор динамической системы; Ns — мощность множества {х.} (параметр алгоритма оценивания размерности аттрактора).
Согласно теореме Такенса [6], ставшей уже классической, аттрактор динамической системы, одна из проекций траектории которой регистрируется в эксперименте, имеет те же метрические свойства, что и восстановленный аттрактор, построенный на основе этой проекции.
Конструктивно восстановленный аттрактор определяется как
х,. = (/„ (v), / „ +x(v),..., / „ +Т(Я— Vi е [p?d},
(v )),
где т — размерность лагового пространства, или пространства вложения (параметр алгоритма оценивания размерности аттрактора).
На основе оценок спектральной плотности гармоники из множества мы строим восстановленный аттрактор как множество т-мерных векторов {х.}. Элементы каждого вектора х. получены как подмножество из т равноотстоящих оценок спектральной плотности выбранной гармоники. Расстояние т между ними (время задержки) — это еще один параметр алгоритма оценивания размерности аттрактора.
Единого критерия выбора параметров е, N , т, т не существует. Авторы работ, в которых получены конкретные оценки размерностей аттракторов различных характеристик медико-биологических систем, практически единодушны в том, что априори можно указать лишь домен К с Я4 возможных значений этих параметров. Оценки Л2 можно затем получить в качестве решения задачи о нахождении локального
оптимума на этом компакте K из четырехмерного евклидова пространства.
Задача определения домена допустимых значений величины задержки т решалась нами в следующей последовательности.
1. Сначала мы оценили первые корни автоковариационных функций анализируемых временных рядов {I nd (v)}, Vv е {vd}. Эти корни можно использовать в качестве супремума домена допустимых значений т. Инфимум мы выбрали равным единице, руководствуясь теми же соображениями, что и авторы статьи [7]. В ходе анализа реальных данных установлено, что оптимальные оценки т локализованы в окрестности инфимума.
2. При выборе домена допустимых значений размерности вложения (лагового пространства) m, мы исходили из тех соображений, что аттракторы исследуемых нами динамических систем почти наверное являются мультифракталами. следовательно, важно выбрать прежде всего супремум для этого домена значений. Наш первоначальный алгоритм построения оценки этого супремума был реализацией хорошо известного алгоритма k ближайших соседей (k-nearest neighbors algorithm, k-NN) [8]. Доставляемые им оценки m ~ 20 в ходе последующего анализа реальных данных были снижены вдвое. В нашем алгоритме m е [1,2, ...,10}. При обсуждении результатов мы покажем, что такой компакт вполне адекватен нашим реальным данным, так как оценки размерности аттракторов не превышали его медиану.
3. При построении домена допустимых значений s критичным является его ин-фимум. именно он определяет дисперсию оценки D2. Если провести рассуждение на языке геометрии, то при слишком малых значениях инфимума s мы получаем большое количество равных нулю ординат оценки C(s). Причина состоит в том, что при малом количестве ординат (малое количество участков p"nd, на которых выполнялось оценивание) числитель оценки из интеграла (2), а именно
Z H (s — I Iх,-— xj II)'
обращается в нуль для большого множества значений 8. Как следствие, для таких значений 8 мы не можем получить оценку С (8), т. е. вынуждены иметь дело с данными, которые в теории оценивания называют данными с пропусками.
Поскольку для получения оценки Б2 требуется идентифицировать ту область значений 8, где достигается плато зависимости 1п С (8) от 1п 8, получаемые оценки С (8) (см. выражение (2)) с пропусками (мы трактуем как пропуски те точки, в которых оценки С(8) не были получены) увеличивают дисперсию оценки Б2. В ходе исследования мы остановились на следующем компакте: 8е[0,1 , 1], шаг был выбран порядка 0,02. При меньших значениях инфи-мума и/или шага возникает вышеуказанная проблема оценивания Б2 по данным с пропусками.
После того, как мы определили домен допустимых значений 8, выбрать домен для N8 не составляло труда, поскольку в подавляющем большинстве случаев нам пришлось анализировать практически весь ряд оценок, которые мы получали по всей траектории процесса. Здесь мы столкнулись с фундаментальной проблемой, когда при огромных объемах входных данных (количество зарегистрированных ординат сигнала ~ 105 — 106) мы, тем не менее, не можем располагать достаточным количеством информации для изучения размерности аттракторов.
4. Заключительный этап алгоритма построения оценки Б2 — это осреднение оценок Б2(8). Вычислив Б2 для каждого т е {2,3, ...,10} и 8е[0,1 , 1], получим зависимости Б2 от е, т, геометрическое представление которых — так называемый график Раппа (рис. 3). Он позволяет получить финальную оценку Б посредством осреднения Б2( 8) при каждом фиксированном т на участке 8Ь < 8 < 8е. Участок, по которому выполняется осреднение, должен удовлетворять следующему условию: оценка коэффициента наклона для линейной регрессионной зависимости Б2( 8) на этом участке статистически незначима. Следуя тем же рассуждениям, что и авторы [9], запишем
выражения для оценки среднего значения
Б 2(т, 8„, 8е):
8е
£ Б( т8*)
Б2(т, 8ь , 8е ) = ^-, (3)
£ 1п
8к=8ь к-1
где
1п
С (т, 8^)
Б(т,8^) = . (4)
Сравнивая графики а и Ь на рис. 3, нетрудно заметить, что вариабельность оценок весьма значительная. Прежде всего, она определяется психоэмоциональным состоянием обследуемого волонтера. То обстоятельство, что он находился в положении лежа на спине, в состоянии полного покоя, вовсе не означает, что это необходимое и достаточное условие, при котором системы регуляции кровообращения обеспечивают стационарное состояние всей системы.
Результаты и их обсуждение
Разработанный нами алгоритм тестировался на модельных данных: аттракторы Икеды, Хенона и Лоренца. Объем модельных траекторий выбирали соизмеримым с объемами оценок плотности спектральной мощности, которые мы анализировали впоследствии (~ 102). Уклонения оценок Б2, полученных по нашему алгоритму, отличались от теоретических не более чем на 5 %.
нами получены оценки плотности спектральной мощности регистрируемых сигналов для двух вариантов выходных данных алгоритма дискретного преобразования Фурье (ДПФ):
с прямоугольным окном и последующим осреднением непараметрическим методом Уэлча;
с окном Хэмминга.
Результаты, полученные для двух вариантов таких данных, сопоставимы. Оценки для здоровых добровольцев приведены в табл. 2.
Для пациентов с патологиями внутри-мозгового кровообращения наблюдался
Б
к -1
а)
2.00 1.80 1.60 1.40 1.20 1.00 0.80 0.60 0.40 0.20
Л
•А' *
м
I 'V
v
V
М
0.00
0.89 0.76 0.65 0.56 0.48 0.41 0.35 0.30 0.26 0.22 £
Ь)
1.60
Рис. 3. Графики Раппа, полученные по итогам анализа сигналов ЛСК, зарегистрированных
в ходе первого (а) и второго (Ь) обследований одного из здоровых волонтеров. Номера кривых соответствуют размерностям вложения. Временной интервал между двумя сеансами
составлял полгода
более существенный сдвиг оценок корреляционной размерности. Например, на частоте 1,250 Гц значения указанных оценок изменялись от 1,82 и 2,15 в состоянии покоя до 1,62 и 1,32 в ходе теста на гипокап-нию. Этот вывод нельзя обобщать на всех
пациентов. Для одних не было выявлено значимых изменений в оценках размерности, в то время как для других пациентов, с патологией, оценки размерности могли увеличиться почти в два раза. Поскольку эти результаты не давали ответа на наш
Примечания. 1. В ячейках таблицы указаны две оценки для сигналов ЛСК: слева — зарегистрированных левым (ВРУ/), справа — правым (ВРУг) доплеровскими датчиками. 2. Прочерки соответствуют случаям, когда имели место данные с пропусками.
таблица 2
Результаты оценки размерности аттракторов (оценки построены для двух гармоник)
Состояние здорового волонтера Окно Оценки размерности (ВРУ// ВРУг)
1,000 Гц 1,125 Гц
Покой Прямоугольное 1,50/2,10 1,40/2,00
Хэминга 2,20/1,90 - / -
Во время теста на гипокапнию Прямоугольное 1,40/1,50 1,50/1,50
Хэминга 1,55/1,88 - / -
главный вопрос, а именно — насколько ро-бастны оценки корреляционной размерности, доставляемые нашим алгоритмом, мы не ограничили наше исследование анализом оценок ) 2, построенных на основе оценки плотности спектральной мощности только одной гармоники из множества {^}. Далее изучались оценки ) 2, построенные для пар гармоник из этого множества (табл. 3).
как следует из наших результатов, ро-бастность оценок намного ниже, чем ожидалось. Наш вывод: расхождение оценок в 10 — 20 % нельзя рассматривать как значимое. В среднем по всему множеству, в ходе теста на гипокапнию у здоровых добровольцев наблюдается уменьшение оценки корреляционной размерности. При этом следует помнить, что наиболее существенные изменения в спектре сигналов в ходе тестов на гипокапнию наблюдаются в ин-франизкой области, на частотах порядка 0,125 Гц (изменяется частота дыхательных экскурсий). Другими словами, несмотря на удаленность гармоник {^}, подвергавшихся анализу, и гармоник, несущих основную информацию о конкретном физиологическом состоянии обследуемого, изменения состояния в ряде случаев можно зарегистрировать и на частотах }. Это наблюдение никак нельзя объяснить эффектом растекания мощности, поскольку мощность гармоник на частотах порядка 0,125 Гц по крайней мере на два порядка меньше мощности гармоник из множества }.
Заключение
Построенный нами алгоритм и вычисленные с его помощью оценки корреляционной размерности детерминированных гармоник позволяют сделать следующие выводы.
Несмотря на кажущиеся огромными входные потоки данных (в нашем случае — оценки ЛСК и АД), их оказывается недостаточно для построения робастных оценок. Причина еще и в том, что в ходе сеансов количественные характеристики физиологического состояния обследуемых могут изменяться очень значительно (см. рис. 1), даже когда их оставляют в полном покое. Поведение первой и второй моментных функций может послужить критерием для исключения из рассмотрения участков траекторий сигналов, на которых не выполняется гипотеза о неизменности наблюдаемой динамической системы.
В последующих исследованиях мы планируем построить входной поток оценок плотности спектральной мощности таким образом, чтобы свести к минимуму эффект растекания мощности, который является неизбежным следствием применения алгоритмов ДПф (см. табл. 1).
Что же касается алгоритма построения компакта допустимых значений параметров 8, , т, т, то, на наш взгляд, он обеспечивает решение поставленной задачи, и в ближайшее время вряд ли будет нами пересмотрен.
Согласно полученным в настоящей работе результатам, в качестве априорной
Таблица 3
Результаты оценки размерности аттракторов (оценки построены для пар гармоник)
Пара гармоник, Гц Оценки размерности (ВРУ// ВРУ г)
Состояние покоя Тест на гипокапнию
Прямоугольное окно Окно Хэминга Прямоугольное окно Окно Хэминга
Данные для здорового волонтера
1,000 1,125 1,52/1,99 2,42/2,19 1,65/1,66 1,67/1,60
1,000 2,000 1,90/2,17 2,32/2,90 2,00/2,00 2,22/2,90
1,875 2,000 1,61/1,98 2,00/2,78 1,76/1,73 2,10/1,78
2,875 3,000 2,17/2,67 2,16/2,00 1,86/1,74 1,72/1,58
1,000 3,000 2,00/2,94 2,50/2,50 2,15/2,80 2,18/2,13
Данные для пациента с патологией
1,000 1,125 1,62/2,00 - / - 1,98/2,17 - / -
1,125 2,125 1,92/2,18 - / - 1,14/ - - / -
2,125 2,250 2,54/1,50 - / - - /2,35 - / -
3,125 3,250 2,30/2,30 - / - 2,79/2,42 - / -
1,125 3,125 2,55/2,00 - / - - / - - / -
Примечания. 1. В ячейках таблицы указаны две оценки для сигналов ЛСК: слева — зарегистрированных левым (ВРУ/), справа правым (ВРУ) доплеровскими датчиками.
2. Оценки получены для результатов ДПФ в двух вариантах: с окном Хэмминга, с прямоугольным окном и последующим осреднением непараметрическим методом Уэлча.
3. Обследования проводились как в состоянии покоя человека, так и во время теста на гипокапнию.
модели для основных детерминированных гармоник модели (1) (эти гармоники вносят основной вклад в мощность сигнала и формируют пульсовую волну ЛСК) можно выбрать модель динамической системы со сравнительно небольшим числом степеней свободы. Но свойства оценок, получаемых в ходе обследования пациента, находятся в очень сложной зависимости от многочисленных факторов, в том числе от психоэмоционального состояния обследуемого. та-
кие оценки следует рассматривать скорее в качестве исходных данных при построении референсных интервалов для конкретной группы обследуемых.
Для реализации изложенных в данной статье алгоритмов нами использовалось свободно распространяемое программное обеспечение с открытым исходным кодом (среда разработки СодеВ1оскз, компиляторы §сс).
Работа выполнялась как инициативная без какой-либо финансовой поддержки.
СПИСОК ЛИТЕРАТУРЫ
[1] Advanced Methods and Tools for ECG Data Analysis. Eds. G.D. Clifford, F. Azuaje, P.E. McSharry. London, Boston: Artech House, 2006. 384 p.
[2] Advanced Biosignal Processing. Advances in Neuroelectric and Neuromagnetic Methods. Ed. by Nait-Ali A. Berlin, Heidelberg: Springer-Verlag, 2009. 378 p.
[3] Ivankov A.A., Samarin A.D. On stationarity of transcranial blood flow velocity registered with Doppler sonography methods // Университетский научный журнал. Физико-математические, технические и биологические науки. 2015. № 13. C. 31-41.
[4] Moran P.A.P. The random division of an interval. Part II // J.R. Statist. Soc. Ser. B. 1951. Vol. 13. No. 2. Pp. 147-150.
[5] Renyi A. on measures of information and entropy // Proceedings of the 4th Berkeley Symposium on Mathematics, Statistics and
Probability. 1960, 1961. Pp. 547-561.
[6] Takens F. On the numerical determination of the dimension of an attractor // Dynamical systems and bifurcations (Eds. B.L.J. Braaksma, H.W. Broer and F. Takens). Lect. Notes in Math. 1125. Heidelberg: Springer, 1985. Pp. 99-106.
[7] Widman G., Schreiber T., Rehberg B., et al. Quantification of depth of anesthesia by nonlinear time series analysis of brain electrical activity // Phys. Rev. E. 2000. Vol. 62. No. 4. Pp. 4898-4903.
[8] Kennel M.B., Brown R., Abarbanel H.D.I. Determining embedding dimension for phase-space reconstruction using a geometrical construction // Phys. Rev. A. 1992. Vol. 45. No. 6. Pp. 3403-3411.
[9] майоров о.Ю., Фенченко в.Н. Вычисление корреляционной размерности и энтропии ЭЭГ сигналов на кластерных вычислительных системах // Клиническая информатика и телемедицина. 2014. Т. 10. Вып.11. С. 10-20.
СВЕДЕНИЯ ОБ АВТОРАХ
ивАНков Алексей Александрович - кандидат физико-математических наук, доцент кафедры прикладной математики Санкт-Петербургского политехнического университета Петра Великого. 195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 [email protected]
куТуЗов Антон валерьевич - студент Института прикладной математики и механики Санкт-Петербургского политехнического университета Петра Великого.
195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 [email protected]
кРивоРУЧко Людмила анатольевна - лаборант кафедры прикладной математики Санкт-Петербургского политехнического университета Петра Великого.
195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 [email protected]
Ivankov A.A., Kutuzov A.V., Krivoruchko L.A. THE ADDITIVE MODEL OF USDG RECORDS: ATTRACTOR DIMENSIONS OF DETERMINATED ADDENDS.
The paper presents the algorithm and the mathematical support to estimate the correlation dimensions of attractors manifested themselves in the most powerful harmonics of spectra retrieved during Doppler sonography assessment of human cerebral circulation. The objective of the present study was to estimate the correlation dimensions for power spectral densities of the most powerful harmonics. The results of real clinical data analysis are given. The estimates of the correlation dimensions were formed for the different examined body's states: in the rest and under test for hypocapnia. The significant relation was found between those estimates of the correlation dimensions and psychoemotional state of volunteer.
CORRELATION DIMENSION, CHAOTIC ATTRACTOR, FOURIER SPECTRUM, POWER SPECTRAL DENSITY.
REFERENCES
[1] Advanced Methods and Tools for ECG Data Patrick E. McSharry, Artech House, London, Analysis, Eds. Gari D. Clifford, Francisco Azuaje, Boston, 2006.
[2] Advanced Biosignal Processing, Advances in Neuroelectric and Neuromagnetic Methods, Pd. by Nait-Ali A. Springer-Verlag, Berlin Heidelberg, 2009.
[3] A.A. Ivankov, A.D. Samarin, On stationarity of transcranial blood flow velocity registered with Doppler sonography methods, Humanities and Science University Journal. Physical and Mathematical, Biological and Technical Science. No. 13 (2015) 31-41.
[4] P.A.P. Moran, The random division of an interval, Part II , J. R. Statist. Soc. B. 13 (2) (1951) 147-150.
[5] A. Renyi, on measures of information and entropy, Proceedings of the 4th Berkeley Symposium on Mathematics, Statistics and Probability, 1960 (1961) 547-561.
[6] F. Takens, on the numerical determination of the dimension of an attractor, Dynamical systems
and bifurcations (Pds. B.L.J. Braaksma, H.W. Broer and F. Takens). Lect. Notes in Math. 1125, Springer, Heidelberg, 1985. Pp. 99-106.
[7] G. Widman, T. Schreiber, B. Rehberg, et al., Quantification of depth of anesthesia by nonlinear time series analysis of brain electrical activity, Phys. Rev. P. 62 (4) (2000) 4898-4903.
[8] M.B. Kennel, R. Brown, H.D.I. Abarbanel, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Phys. Rev. A. 45 (6) (1992) 3403-3411.
[9] O.Yu. Mayorov, V.N. Fenchenko, Vychisleniye korrelyatsionnoy razmernosti i entropii PPG signalov na klasternykh vychislitelnykh sistemakh [Calculation of correlation dimension and entropy of PPG signals on the cluster computer system], Klinicheskaya informatika i Telemeditsina. 10 (11) (2014) 10-20.
THE AuTHORS
IVANKOV Alexei A.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
KUTUZOV Anton V.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
KRIVORUCHKO Lyudmila A.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
© Санкт-Петербургский политехнический университет Петра Великого, 2016
