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

УДК 519.246.8, 612.829

М. И. Богачёв, К. Е. Громова, Д. М. Клионский, А. С. Красичков, В. Н. Команцев Ю. А. Малых, О. А. Маркелов, С. А. Пыко, В. С. Рамаданов, А. А. Соколова, Ю. Д. Ульяницкий

Санкт-Петербургский государственный электротехнический

университет "ЛЭТИ"

| Флуктуационный анализ физиологических сигналов1

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

Физиологические сигналы, нестационарная динамика, медленные колебания, флуктуационный анализ, дифференциальная диагностика

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

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

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

1 Работа выполнена при поддержке Минобнауки РФ в рамках ФЦП "Научные и научно-педагогические кадры

инновационной России на 2009-2013 годы" (Государственный контракт № 14.740.11.0577 от 05.10.2010). © Богачёв М. И., Громова К. Е., Клионский Д. М., Красичков А. С., Команцев В. Н., 37

Малых Ю. А., Маркелов О. А., Пыко С. А., Рамаданов В. С., Соколова А. А., Ульяницкий Ю. Д., 2012

Известия вузов России. Радиоэлектроника. 2012. Вып. 6======================================

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

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

Исследование медленной регуляции физиологических ритмов и попытки описания их динамики активизировались на рубеже 1990-2000-х гг. (см., например, [3]). В ней и в ряде других работ показано, что для многих физиологических ритмов характерно отсутствие локализации регуляторных колебаний около некоторой фиксированной частоты. Напротив, реализация достаточно сложных нелинейных автономных регуляторных механизмов приво-

-х масштабов (и, соответственно, в широком диапазоне частот). C одной стороны, показана диагностическая информативность данных показателей, в первую очередь при диагностике серьезных функциональных и регуляторных нарушений. В частности, выявлена диагностическая ценность оценки снижения адаптационных возможностей на основе ослабления нелинейной составляющей зависимости между отсчетами сердечного ритма в широком временном диапазоне для оценок риска инфаркта миокарда и развития сердечной недостаточности (см., например, [4]). С другой стороны, показана относительная устойчивость масштабно-инвариантных свойств процесса и параметров его медленной регуляции с возрастом и при невыраженных сопутствующих нарушениях, их слабая зависимость от физической активности и от изменения симпатической и вагусной регуляций, например при кратковременной фармакологической стимуляции симпатической активности [5]. Также отмечена информативность свойств нелинейных составляющих медленных регуляторных колебаний при анализе зависимости динамики сердечного ритма от различных фаз сна, когда в фазе "быстрого" сна медленная регуляция была более выражена, чем в фазе "медленного" сна [6].

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

Метод флуктуационного анализа в базисе Хаара (метод мультифрактального анализа на основе вейвлет-преобразования (multi fractal wavelet transform analysis - MF WTA)) реализуется в соответствии со следующим алгоритмом (см. например, [7]). Доступная ис-

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

После сегментации в каждом из сегментов вычисляются локальные суммы:

к$

Ук= Ъ

/= к-1 5

На следующем этапе берутся первые (вторые и т. д.) разности (для методов MF \VTA1, МБ \УТА2 и т. д. соответственно): у'к ~Ук~Ук-Ь Ук ~ Ук ~Ук-Ъ •••> после чег0 вычисляются флуктуационные функции различных порядков q усреднением по всем к ок-

1 ч

, где • - операция взятия целой части аргумента; т -

нам: Fq s =

, 2 N/s

— z

2N ¿i

Укт

порядок используемой производной ряда yk, соответствующий порядку метода. Типичный диапазон длины окна приближенно определяется как 1 <s<N/8; коэффициент 2 в знаменателе и в верхнем пределе суммы используется только при реализации алгоритма с двойным проходом.

Флуктуационная функция F s характеризует линейную зависимость флуктуаци-онной составляющей динамического ряда, а флуктуационные функции порядков q ^ 2 отражают нелинейные взаимосвязи отсчетов динамического ряда. Существует взаимно однозначное соответствие между F s и автокорреляционной функцией динамического

ряда С ^ ~ , выраженное соотношением F2 s

Метод флуктуационного анализа с исключением тренда (detrended fluctuation analysis - DFA), впервые предложенный в работе [8], предполагает вычисление профиля (кумулятивной суммы) элементов анализируемого временного ряда xt : Ys = ^ х, - х, который за-

i=1

тем подразделяется на N/s сегментов равной длительности s отсчетов. Далее полиномиальной аппроксимацией порядка p каждого из сегментов с помощью метода наименьших квадратов из экспериментальных данных удаляется тренд, после чего последовательным изменением размеров окна 5 в диапазоне р + 2 < s < N/4 для всех возможных вре х локализаций v вычисляется флуктуационная функция:

о 1 s 2

F s,\ =-У Y v — 1 s+k -Pv к ,

2 N/s tl V

где Pv к - значения аппроксимирующего полинома для v-ro окна длины

Порядком аппроксимирующего полинома определяется порядок метода (метод DFA0 предполагает вычитание только постоянного значения, DFA1 - линейного тренда и т. д.). Результирующая флуктуационная функция вычисляется усреднением по всем вре-

м положениям v: F s =

1 2 N/. ^

I F2 s, V

. Аналогично предыдущему

.2 Ы/х У=1

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

Наклон а графика функции Е s , представленной в двойном логарифмическом масштабе, отражает скорость убывания корреляции при удалении точек. Указанный наклон соотносится с показателем наклона корреляционной функции, как а = 1-у/2. В отсутствие корреляции или при наличии только кратковременной корреляции 0 < а < 0.5 соответствует отрицательной долговременной корреляции, 0.5 < а < 1 - положительной.

Для характеристики нелинейной долговременной зависимости (анализа флуктуаци-онных функций порядков с/ > 2), а также для дробных моментов предложен мультифрак-тальный DFA-метод (МР DFA) [9], в рамках которого вычисляется флуктуационная функ-

IV ч

F \ 1 Y Г/г*

ция порядка q\ Fq s 1 v

12 7V/s V=1

q/2

В указанной работе установлено, что для широкого класса случайных процессов, называемых фрактальными, при больших s значение Fq s асимптотически пропорционально s q, где Hq - обобщенный показатель Хёрста динамического ряда x. Для процессов с исключительно линейной зависимостью Hr. =Н = const.

ч

Графическое отображение результатов исследования удобно представлять функцией Fq s/.У2, поскольку при этом отсутствие зависимости между отсчетами динамического

ряда соответствует горизонтальной линии Fq s j.У2 = const, Vs, q. При изучении процессов с отрицательными значениями Hq для некоторых q обычно прибегают к двойному интегрированию вместо одинарного, т. е. повторяют операцию вычисления кумулятивной суммы дважды, и в дальнейшем визуализируют функцию Fq s j's3'2. Характерные искажения флуктуационных функций в присутствии трендов сложной формы, регулярных (квазипериодических) и случайных составляющих рассмотрены в работе [10].

На рис. 1 представлены результаты применения рассмотренных методов для анализа электроэнцефалографических (ЭЭГ) и электромиографических данных. На рис. 1, а приведены флуктуационные функции для ЭЭГ в норме, на рис. 1, б - при патологии (лейко-энцефаломиелит без эпилептических приступов), полученные с помощью метода MF WTA1. На рис. 1, а в структуре флуктуаций в зоне I (колебания с постоянной времени 0.01 < 5 < 0.1 с) доминирует альфа-ритм с характерной постоянной времени 5 « 0.05 с (что на графиках отображается характерным "горбом" кривых). В зоне II наблюдается устойчивый рост Fq s с Hq «0.8...0.9, что соответствует достаточно выраженной линейной

зависимости в стационарном режиме. Необходимо отметить, что в зонах I и II Fq s

0.5 10"

1 1 1 д = 5 2 1 1 1 1 1 1

— 1 ^Л 1 2

/ I 1 1 ! II 1 1 ! ш ! 1 1 1 0.2

10"

10"

10°

10"

10"

10"

100

5, с

0.2 10

1 <7 = 5 Ууг

' тт

ш |

/ I 1 II ' и

-з

10

-2

10"

100

10"

10"

10"

100

Рис. 1

практически не зависит от q, а это свидетельствует об отсутствии выраженных нелинейных эффектов в исследуемом сигнале и о доминирующей линейной зависимости. Напротив, в зоне III (характеризующей колебания с постоянными времени 5 = 1.. .10 с) наблюдается расхождение функций ^ 5 для различных значений q, что является признаком

присутствия нелинейной зависимости.

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

отсчетами реализации, индикатором которой является спад функции ^ ^ .

На рис. 1, в приведен пример флуктуационных функций электромиографического сигнала при удержании массы 0.5 кг, на рис. 1, г - массы 1.5 кг, полученных методом MF DFA2

м масштабам ^ < 0.03 с,

вид флуктуационной функции идентичен для различных уровней нагрузки и отмечается квазипериодический процесс с характерной постоянной времени 5«0.03 с. В зоне II 0.03 < л < 0.1 с также наблюдается положительная корреляция, при этом различия в наклоне функций незначительны. Наиболее существенные различия наблюдаются в зоне III 5 > 0.1 с , где характер корреляции флуктуационной составляющей изменяется с положительного (см. рис. 1, в) при меньшей нагрузке на отрицательный (см. рис. 1, г) при большей нагрузке. При этом при малой нагрузке (см. рис. 1, в) также отмечается нелинейная зависимость (отражающаяся в различном наклоне флуктуационных функций для раз-

41

5

з

2

3

2

5, с

б

а

2

з

2

5, с

я, с

в

г

Известия вузов России. Радиоэлектроника. 2012. Вып. 6======================================

личных моментов q). При увеличенной нагрузке (см. рис. 1, г) указанные различия практически незаметны, однако наблюдается рост Fq s в области .v > 0.5 с, связанный, по-

видимому, с проявлением тренда.

Альтернативой рассмотренным ранее методам может служить подход, основанный на декомпозиции на эмпирические моды (ДЭМ), впервые предложенный в [11]. Эмпирическая мода (ЭМ) (intrinsic mode function - IMF) - функция, сформированная на основе значений сигнала, заданная непрерывно на интервале существования сигнала или дискретно в виде вектора отсчетов, имеющая в общем случае произвольную форму и произвольную аналитическую запись (если таковая существует), но всегда удовлетворяющая двум необходимым условиям:

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

2. Локальное (мгновенное) среднее значение функции, равное полусумме верхней огибающей, полученной интерполяцией найденных локальных максимумов, и нижней огибающей, полученной интерполяцией найденных локальных минимумов, не должно превышать заранее определенного порогового значения г|, зависящего от машинной точности в, и погрешностей, связанных с получением, преобразованием и передачей сигнальной информации. В качестве средства интерполяции чаще всего используются кубические сплайны. Добиться точного равенства нулю локального среднего значения в каждый момент времени невозможно по ряду объективных причин: вычислительных погрешностей, плохой обусловленности систем уравнений для расчета коэффициентов сплайнов, а также особенностей самого сигнала (например, краевых эффектов - сильных осцилляций интерполирующей функции на краях, благодаря которым вблизи краев локальное среднее практически всегда отличается от нуля).

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

Максимально возможное число ЭМ определяется как М «log2 А''±1, где N - число

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

Рис. 2 иллюстрирует применение метода ДЭМ к обработке поверхностной элек-тромиограммы, снятой с мышцы предплечья ^xtensor digitorum) испытуемого при выполнении им статической нагрузки в виде удержания груза массой 0.5 кг. На рис. 2 представлены исходный сигнал х, первые четыре эмпирические моды IMF^-.-IMF\ и остаток Ах, полученный вычитанием из сигнала сформированных ЭМ.

х-104 2 0 - 2 - 4

0.1

тщ- ю'

0.2

0.3

0.4

0.5

г, с

- 2 - 4

I I I I ■ I

ЦЦ

0

0.1

0.2

0.3

0.4

0.5

г, с

тщ- ю'

1.5 0

- 1.5

- 3

4

0

0.1

0.2

0.3

0.4

0.5

г, с

тщ- ю5

5

- 5 - 10 - 15

ллмлЛшкллАА 11<||||Ь1|а ЛалиишА^ьУк!^ 1ыЛиш/ии|1ьм11||| Аш^Л л и 1/1 м|Л

1ШЛ-10 2.1 0

- 2.1

- 4

Ах-10 1.5 0

- 1.5

о

5

0.1

0.2

0.3

0.4

0.5

г, с

0.1

0.2

0.3

0.4

0.5

г, с

,5

0

0.1

0.2

0.3

0.4

0.5

г, с

Рис. 2

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

основных ритмов электроэнцефалограммы - дельта-ритма тета-ритма 1\), альфа-рит-

2 „ ма и бета-ритма /р . Как визуальное сопоставление, так и результаты статистической

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

0

0

0

2 Дельта-ритм сосредоточен в области частот 0.5...4 Гц, тета-ритм занимает область частот 4...8 Гц. Альфа-ритм сосредоточен в области 8.13 Гц, а бета-ритм занимает область 13.30 Гц.

ст, мкВ

10

1 2 3 4 5 6 7 8 9 10 11 12 13 n

2

ft, мкВ

150

100

50

Ра, мкВ'

I- LjJ.

I MI И.....

1 2 3 4 5 6 7 8 9 10 11 12 13 n

40

20

-iji« I I -

1 2 3 4 5 6 7 8 9 10 11 12 13 n

Рис. 3

100 80 60 40 20 0 - 20

2

Pn, mkB

60

30

2

Pp, мкВ

20

10

I I I I I I

I I I L i.

1 2 3 4 5 6 7 8 9 10 11 12 13 n

■ l .........

1 2 3 4 5 6 7 8 9 10 11 12 13 n

]ц i l -I

1 2 3 4 5 6 7 8 9 10 11 12 13 n

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

Список литературы

1. Hildebrandt G., Moser M., Lehofer M. Chronobiologie und Chronomedizin. Stuttgart: Hippokrates Verlag in MVS Medizinverlage, 2002. 145 s.

2. Математические методы выявления регулярных и статистических закономерностей в биомедицинских и экологических данных большого объема / М. И. Богачёв, А. Р. Каюмов, А. С. Красичков, О. А. Мар-келов. СПб: Изд-во СПбГЭТУ "ЛЭТИ", 2012. 176 с.

3. Scaling behaviour of heartbeat intervals obtained by wavelet-based time-series analysis / P. C. Ivanov, M. G. Rosenblum, C.-K. Peng et al. // Nature. 1996. Vol. 383. P. 323-327.

4. Cardiac interbeat interval dynamics from childhood to senescence: comparison of conventional and new measures based on fractals and chaos theory / S. M. Pikkujamsa, Т. H. Makikallio, L. В. Sourander et al. // Circulation. 1999. Vol. 100. P. 393-411.

5

0

0

0

0

0

5. Schmitt D. Т., Ivanov P. C. Fractal scale-invariant and nonlinear properties of cardiac dynamics remain stable with advanced age: a new mechanism picture of cardiac control in healthy elderly // Am. J. physiol. regul. integr. comp. physiol. 2007. Vol. 293. P. 1923-1937.

6. Comparison of detrended fluctuation analysis and spectral analysis for heart rate variability in sleep and sleep apnea // Т. Penzel, J. W. Kantelhardt, L. Grote et al. // IEEE trans. biomed. eng. 2003. Vol. BE-50. P. 1143-1156.

7. Feder J. Fractals. New York: Springer, 1988. 310 p.

8. Mosaic organization of DNA nucleotides / C.-K. Peng, S. V. Buldyrev, S. Havlin et al. // Phys. rev. E. 1994. Vol. 49. P. 1685-1689.

9. Multifractal detrended fluctuation analysis of nonstationary time series / J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde et al. // Physica (Amsterdam). 2002. Vol. 316. P. 87-114.

10. On spurious and corrupted multifractality: the effects of additive noise, short-term memory and periodic trends / J. Ludescher, M. I. Bogachev, J. W. Kantelhardt et al. // Physica (Amsterdam). 2011. Vol. 390. P. 2480-2490.

11. The empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis / N. E. Huang, Z. Shen, S. R. Long et al. // Proc. royal soc. A. London. 1998. Vol. 454. P. 903 -995.

12. Sokolova A., Bogachev M. I., Bunde A. Clustering of ventricular arrhythmic complexes in heart rhythm // Phys. rev. E. 2011. Vol. 83. P. 021918 (1-7).

M. I. Bogachev, K. E. Gromova, D. M. Klionsky, A.S. Krasichkov, V. N. Komantsev Y. A. Malych, O. A. Markelov, S. A. Pyko, V. S. Ramadanov, A. A. Sokolova, Y. D. Uljanitski Saint-Petersburg state electrotechnical university "LETI"

Fluctuation analysis of physiological signals

Several fluctuation analysis methods with application to physiological signals are considered that are capable of detecting and analyzing fluctuations caused by slow regulatory feedback laps in physiological systems in the presence of non-stationary components (trends), regular quasi-periodic and noisy components are considered. The efficiency of these methods is demonstrated using examples of electroencephalographical and electromyographical signals analysis.

Physiological signals, non-stationary dynamics, slow fluctuations, fluctuation analysis, differential diagnostics

Статья поступила в редакцию 3 декабря 2012 г.

УДК 621.396.4

Е. Б. Соловьёва, А. В. Зубарев

Санкт-Петербургский государственный электротехнический

университет "ЛЭТИ"

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

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

Импульсный шум, речевой сигнал, нелинейный фильтр, метод расщепления сигналов

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

© Соловьёва Е. Б., Зубарев А. В., 2012 45