Анализ фазовых портретов биомедицинских сигналов

Алгоритмы нелинейной динамики находят применение как при моделировании работы сердечнососудистой системы, головного мозга, так и для диагностики состояний различных физиологических систем¹².

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

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

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

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

Введем некоторые определения из теории нелинейной динамики^{34 5}.

Основные понятия и определения

Если система описывается конечным набором n параметров, то ее целесообразно рассматривать в абстрактном пространстве, оси которого представлены переменными x₁, . ..,x_n. Такое n-мерное пространство получило название фазовое пространство. Каждому состоянию динамической системы соответствует точка в этом пространстве, и каждой точке из этого пространства соответствует единственное состояние системы. Изменение состояния системы можно интерпретировать как движение точки в фазовом пространстве. Траектория такой точки, то есть последовательное изменение ее положения в фазовом пространстве, называется фазовой траекторией. Множество траекторий в фазовом пространстве дают общую картину поведения исследуемой системы, т. е. ее фазовый портрет. В общем случае фазовые траектории с течением времени стягиваются к некоторой области фазового пространства, называемой аттрактором. Аттракторы бывают статические (фиксированная точка), периодические (предельный цикл), квазипериодические (тор) и хаотические. Последний тип аттрактора называют странным аттрактором — он характеризуется дробной фрактальной размерностью.

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

В хаотических системах фазовые траектории с течением времени стягиваются к странному аттрактору. Но на самом аттракторе движение является неустойчивым: любые две траектории системы расходятся экспоненциально, оставаясь на странном аттракторе. При этом поведение их очень похоже на случайное, хотя динамика детерминирована и воспроизводима при условии точного повторения начальных условий. Классическими примерами таких систем являются отображение Хенона и аттрактор Лоренца⁶.

Динамические системы, имеющие несколько аттракторов, при незначительных внешних воздействиях могут «перескакивать» с одного аттрактора на другой. Такой скачок называется бифуркацией, а момент времени, в который он произошел, — точкой бифуркации. При бифуркациях происходит значительная качественная перестройка динамики системы. Анализу бифуркаций в нелинейной динамике

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

Наиболее широко распространенная величина, используемая при описании хаотической динамики — размерность аттрактора. В общем случае она равна минимальной размерности пространства, содержащего все фазовые траектории системы. Другими словами — это число степеней свободы динамики системы. Существует большое количество численных методов нахождения размерности аттрактора. Самый используемый алгоритм предложен Grassberger и Procaccia ⁷⁸⁹. Он позволяет вычислить так называемую корреляционную размерность D_с аттрактора, характеризующую степень сложности динамики исследуемой системы. Важность этой характеристики состоит еще в том, что она может быть определена непосредственно по экспериментальным данным. При этом путем обработки только одного временного ряда можно оценить минимальное число переменных, которое должна включать математическая модель исследуемой динамической системы. Кроме того, в результате специальной обработки регистрируемого сигнала становится возможной оценка размерности фазового пространства и построение фазового портрета, по которому проводится анализ состояния физиологической системы.

Реконструкция фазового пространства

Рассмотрим метол реконструкции фазового пространства для системы, состояние которой оценивается только по одной реализации. Аттрактор для такой временной последовательности может быть построен с помощью метода задержек, основанного на теореме Такенса¹⁰. Согласно этой теореме, можно восстановить некоторые свойства аттрактора в фазовом пространстве по временной последовательности только одной составляющей, т. е. наблюдая за изменением лишь одного параметра этой системы.

Пусть состояние некоторой динамической системы оценивается только по одной из фазовых координат, которая представлена временным рядом:

x(1),х(2)…..x(N),                    (2.11)

где N — количество измерений.

Восстановление аттрактора системы сводится к построению псeвдоаттрактора, в котором компонентами вектора, соответствующего точке x(t), являются элементы анализируемой последовательности, но полученные с некоторой временной задержкой:

x(t), x(t + τ)…..x(t) + (m – 1)τ),

где т – величина задержки (дискретная величина); m — размерность пространства вложения; t = ((m — 1)τ + 1, …, N — дискретное время.

В векторном пространстве этих точек (x(t), x(t + τ), …, x(t + (m – 1)τ)) можно построить аттрактор таким образом, что он будет сохранять важнейшие свойства и динамику исходного аттрактора. При этом размерность т определяется из условия

m≥2[d| + l,                     (2.12)

где d – фрактальная размерность аттрактора, [d] — обозначает целую часть фрактальной размерности¹¹.

Основной задачей анализа последовательности (2.10) для последующей пссвдофазовой реконструкции является определение величин временной τ задержки и размерности m.

В работах ¹²¹³ подробно рассмотрено несколько способов выбора временной задержки: метод автокорреляционной функции, метод взаимной информации, метод среднего отклонения. Поскольку компоненты вектора [x(f), x(t + τ) … x(t+(m-1)τ], характеризующего систему, должны быть независимы, величина τ может быть получена как точка первого пересечения с нулевой линией автокорреляционной функции временного ряда. Данный метод является достаточно простым и часто применяется в практических задачах.

Для оценки величины d используют алгоритм Grassberger-Procaccia, позволяющий вычислить корреляционную размерность D_c. Процедура расчета этого показателя включает вычисление корреляционного интеграла С(r) по формуле (2.8), которую в общем случае применяют для векторов Х(i) и X(j), описывающих поведение системы в фазовом пространстве. Величина С (r) равна нормированному числу пар точек на аттракторе, расстояние между которыми не превышает величины r. На специально построенном графике, показывающем логарифмическую зависимость In С(r) = f In(r), выделяют линейный участок (область скейлинга) и методом наименьших квадратов строят линейную зависимость. Корреляционную размерность, определяемую соотношением D_c = lim In С(r) / In(r), вычисляют как тангенс угла наклона этой прямой. Полученная величина и будет оценкой размерности d.

В случае анализа временного ряда (2.10) данная процедура применяется для множества последовательностей

X(t) = (x(i), x(i — τ)…..x(i- (m – 1)τ)), i = ((m – l)τ+ 1)…..N

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

Рис. 1 иллюстрирует построение линейных зависимостей корреляционного интеграла для расчета D_c(m) на примере отображения Хенона. Вычисленные значения корреляционной размерности для разных m сведены в таблицу 1.

Рис. 3. Представление отображения Хенона в виде: а — последовательности отсчетов; б — фазовой траектории на плоскостиТаблица 1

Для определения размерности вложения необходимо последовательно вычислять корреляционную размерность для каждого m = 2,3, … и затем следует построить графическую зависимость D_c = f(m). В начальной области графика при добавлении новых компонент псевдоаттрактора и, соответственно, увеличении m, корреляционная размерность растет. Это означает, что    m компонент недостаточно для сохранения основных свойств аттрактора. Начиная с некоторой размерности m, кривая перестает изменяться, т.е. наблюдается насыщение величины D_c(m) на некотором уровне d. Эта величина d определяет корреляционную (фрактальную) размерность аттрактора, т.е. минимальную размерность фазового пространства. Величина m, с которой начинается насыщение D_c(m), задаст минимальную размерность пространства вложения. Для сильно зашумленного сигнала насыщения D_с(m) нс наблюдается.

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

Важно отметить следующие особенности практического применения рассмотренного алгоритма. Координаты временного ряда должны быть выбраны с шагом, равным величине τ. Причем величина шага между элементами ряда (2.10) и величина задержки τ могут и не совпадать. Кроме того, в практических задачах значение m, определяемое формулой (2.11), часто оказывается завышенным, и можно ограничиться пространством меньшей размерности, а именно, воспользоваться условием m = |d+1|¹⁴ и даже m < d¹⁵.

На рис. 2 в качестве примера приведены графики зависимости корреляционной размерности D_с от размерности вложения m для модельных сигналов, а именно, для гармонического сигнала (рис. 49, а), отображения Хенона (рис. 2, 6) и белого шума (рис. 2, в). Все кривые построены по наблюдаемым последовательностям.

Модельные сигналы заданы в виде:

x(i) = 88 • sin(π • i/100) + 600 гармонический сигнал;

x(i) = 220 • randn(N, 1) + 600 — белый шум с математическим ожиданием 600 и дисперсией – 220, где randn функция, генерирующая случайную величину с нормальным распределением, N — количество отсчетов дискретного сигнала.

Отображение Хенона — двумерный аналог логистического отображения, имеет следующий вид:

x(i + 1) = 1 – a• x(i)² + b • x(i – 1),

где а и b — внешние параметры, равные а = 1,4 и b = 0,3.

Как видно из рисунка, кривая зависимости корреляционной размерности D_c от размерности вложения для гармонического сигнала с увеличением m возрастает и выходит на плато на уровне 1. Для отображения Хенона кривая D_c переходит в насыщение на уровне дробной величины 1,23, что свидетельствует о хаотическом характере генерируемого процесса. Для случайного сигнала (белого шума) увеличение размерности вложения m на единицу приводит к примерно такому же увеличению показателя D_c, поэтому теоретически для его описания требуется бесконечное число переменных.

Рассчитанные значения корреляционной размерности для заданных модельных сигналов, а также значения, найденные в литературных источниках, практически совпадают (см. табл. 2). Небольшие различия в полученных значениях можно объяснить малой длиной анализируемой выборки, которая составляла 1000 отсчетов.

Таблица 2

На рис. 3 приведен вид временной последовательности Х(i), i = 1…..300 для отображения Хенона (рис. 3, а), а также его реконструированный фазовый портрет в пространстве размерности m = 2 (рис. 3, б).

По уравнению движения известно, что аттрактор Хенона является двумерным¹⁶¹⁷, но здесь восстановленное фазовое пространство размерности m = 2 и фазовый портрет получены расчетным путем.

Распознавание аритмий по ритмограмме

Рассмотренный метод реконструкции фазового пространства применялся для анализа биомедицинских сигналов (ЭКГ, ЭЭГ). Проведен анализ корреляционной размерности реальных данных — сигнала сердечного ритма при некоторых видах аритмий¹⁸. В первую очередь, это мерцательная аритмия (МА, atrial fibrillation, фибрилляция предсердий), надежное обнаружение которой является одной из главных задач кардиомониторного наблюдения. При этом нарушении выражена нерегулярность в появлении желудочковых сокращений, наблюдается хаотический разброс длительностей кардиоциклов (рис. 4, а). Анализ ритмограмм показал, что зависимость корреляционной размерности от размерности вложения, как и в случае модельных экспериментов (отображение Хенона), имеет насыщение на уровне дробной величины d = 5,5 (рис. 4, б). Это свидетельствует о хаотической природе сердечного ритма при МА.

Аналогичные зависимости были получены для нормального синусового ритма (HP) и частой экстрасистолии (ЧЭ). Оказалось, что для HP характерно значение d = 3,5, а в случае ЧЭ d = 1. Эго два наиболее распространенных вида ритма сердца, на фоне которых трудно распознавать МА, используя традиционные статистические и спектральные методы. Кривые были построены согласно предположению о независимости смежных отсчетов ритмограммы, что, конечно, позволяет получить лишь приближенную оценку фрактальной размерности d.

В практических задачах выбор размерности т на основе расчетных значений d зависит от типа решаемых задач. Для задач моделирования, прогнозирования временных рядов, исследования свойств аттрактора сердечного ритма размерность пространства вложения необходимо выбирать в соответствии с приведенными выше рекомендациями, ориентируясь на полученные значения d и условие (2.11). Классификация сигналов возможна на основе более простых описаний, например, представлении псевдофазового портрета в пространстве меньшей размерности. При решении таких задач важно выявить отличительные особенности фазового портрета и сформировать набор признаков, обеспечивающий надежное распознавание сигналов заданных классов.

Учитывая большие различия в полученных оценках минимальной размерности фазового пространства, для решения задачи обнаружения МА на фоне HP и ЧЭ выбрана размерность псевдофазового пространства равной 2. На рис. 5 представлены проекции псевдофазовго портрета на плоскость для динамической модели сердечного ритма при МА. HP, ЧЭ.

Анализ особенностей полученных изображений позволил сформировать следующий набор морфологических признаков:

• длина фазовой траектории;
• площадь портрета, ограниченная контуром;
• длина контура;
• угловой показатель.

Первые три признака довольно просто вычисляются и традиционно используются при обработке изображений¹⁹. Четвертый признак определяет относительно число переходов от короткого интервала к длинному и наоборот на всей выборке данных. Именно такие переходы наиболее характерны для сигналов, имеющих слабую регулярность, в частности, для последовательности кардиоинтервалов при МА. Фактически вычисляется число линий фазовой траектории, расположенных под углом 100°- 170° к горизонтальной оси, задаваемой координатой x(i).

Методом линейного дискриминантного анализа, основанном на использовании критерия Фишера, получено решающее правило, которое реализовано в алгоритме автоматического обнаружения МА. Результаты экспериментов, проведенных на выборке реальных записей ЭКГ показали, что ошибка классификации составляет 2% при длине последовательности, равной всего 50 кардиоинтервалам.

Важно отметить, что данный метод позволяет достаточно надежно обнаруживать МА на фоне HP и ЧЭ по коротким фрагментам ритмограммы, что отличает его от ранее рассмотренною метода, использующего оценки аппроксимированной энтропии.

1. Анищенко В. С. Знакомство с нелинейной динамикой. 3-е изд. – М.: Издательство ЛКИ, 2008. — 224 с.
2. Nonlinear biomedical signal processing / ed. by Metin Akay. Vol. 2. Dynamic analysis and modelling. — New York: IEEE, 2001. — 341 p.
3. Шустер Г. Детерминированный хаос: Введение. М.: Мир, 1988. 240 с.
4. Малинецкий Г. Г., Потапов А. Б., Подлазов А. В. Нелинейная динамика: Подходы, результаты, надежды. Изд. 3-е. — М.: Книжный дом «ЛИБРОКОМ», 2011. – 280 с.
5. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
6. Шустер Г. Детерминированный хаос: Введение. М.: Мир, 1988. 240 с.
7. Шустер Г. Детерминированный хаос: Введение. М.: Мир, 1988. 240 с.
8. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
9. Nonlinear biomedical signal processing / ed. by Metin Akay. Vol. 2. Dynamic analysis and modelling. — New York: IEEE, 2001. — 341 p.
10. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
11. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
12. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
13. Малинецкий Г. Г., Потапов А. Б., Подлазов А. В. Нелинейная динамика: Подходы, результаты, надежды. Изд. 3-е. — М.: Книжный дом «ЛИБРОКОМ», 2011. – 280 с.
14. Головко В. А. Нейросетевые методы обработки хаотических процессов // Научная сессия МИФИ-2005: Лекции по нейроинформатике. — М.: МИФИ. 2005. – С. 43-91.
15. Анищенко В. С. Знакомство с нелинейной динамикой. 3-е изд. – М.: Издательство ЛКИ, 2008. — 224 с.
16. Малинецкий Г. Г., Потапов А. Б., Подлазов А. В. Нелинейная динамика: Подходы, результаты, надежды. Изд. 3-е. — М.: Книжный дом «ЛИБРОКОМ», 2011. – 280 с.
17. Шустер Г. Детерминированный хаос: Введение. М.: Мир, 1988. 240 с.
18. Манило Л. А., Зозуля Е. П. Оценка корреляционной размерности временных рядов в задачах анализа аритмий // Биомедицинская радиоэлектроника. 2009, № 11. С. 32-39
19.  Гонсалес P., Вудс P. Цифровая обработка изображений I Пер. с англ. М.: Техносфера, 2005. — 1072 с.