Алгоритм оценки дрейфа изоэлектрической линии кардиосигнала при анализе длительных мониторограмм Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

СПИСОК ЛИТЕРАТУРЫ

1. Амплитудно-фазовая конверсия / Г. М. Крылов, В. З. Пруслин, Е. А. Богатырев и др. М.: Связь, 1979. 256 с.

2. Зенькович А. В. Снижение влияния изменений амплитуды сигнала в фазовых системах // Радиоэлектроника. Изв. вузов. 1972. Т. XV, № 11. С. 1402-1403.

3. А. с. СССР 420951 А1 МКИ6 G01K25/00. Измеритель разности фаз и фазовой модуляции / А. В. Зенькович (СССР). Опубл. 25.03.1974. БИ № 11.

4. Зенькович А. В., Типашов В. И. Определение коэффициента преобразования паразитной амплитудной модуляции в фазовую (частотную) в смесителях // Радиоэлектроника. Изв. вузов. 1974. Т. XVII, № 1. С. 80-85.

A. V. Zenkovich Nizhniy Novgorod state technical university n. a. R. E. Alexeev

Method for reducing the influence of signals amplitude in location and navigation phase systems

New compensating method for reducing the amplitude-to-phase conversion caused by varying of input signals amplitude is presented. Theoretic base and features are given, experimental results and usage guidelines are discussed.

Signals amplitude changes, amplitude-to-phase conversion, phase systems

Статья поступила в редакцию 22 марта 2013 г.

УДК 621.37

А. С. Красичков

Санкт-Петербургский государственный электротехнический университет "ЛЭТИ" им. В. И. Ульянова (Ленина)

Алгоритм оценки дрейфа изоэлектрической линии кардиосигнала при анализе длительных мониторограмм1

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

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

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

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

За сутки при использовании холтеровского мо-ниторирования накапливается более 80 тысяч сердечных циклов, что приводит к огромным временным затратам на расшифровку мониторограмм. Ситуация усугубляется тем, что при съеме элек-

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

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

1 Работа выполнена в рамках договора № 1258/ФКТИ2 от 01.10.2012 г. между ОАО «Концерн "Океанприбор"» и СПбГЭТУ "ЛЭТИ". © Красичков А. С., 2013

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

Наиболее простой подход основан на устранении дрейфа с помощью фильтра верхних частот [2]. Недостатком указанного способа является то, что вместе с дрейфом изолинии из смеси элек-трокардиосигнала (ЭКС) и помехи удаляется часть спектральных составляющих ЭКС. Подавление составляющих в низкочастотной части спектра ЭКС приводит, в свою очередь, к искажению формы ЭКС на выходе фильтра верхних частот, что может повлиять на результаты анализа.

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

Недостатками данного подхода являются следующие:

- не всегда в ЭКС присутствует явно выраженный и стабильно лежащий на изолинии сегмент РР или ТР (например, желудочковые экстрасистолы не имеют РР-сегмента, ТР-сегмент на высокой частоте сердечных сокращений (ЧСС) может вообще отсутствовать);

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

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

- наблюдаются искажения сигнала на сегменте БТ;

- не позволяет полностью ликвидировать дрейф изоэлектрической линии.

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

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

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

При построении статистической модели сигнала приняты следующие допущения:

1. Миографическая помеха представляет собой аддитивный нормальный "белый" шум с нулевым средним значением и спектральной плотностью мощности 0.5^0. Нормальное приближение шумов электрической активности мышц следует из центральной предельной теоремы как результат сложения большого числа электрических импульсов, среди которых нет превалирующих. Таким образом, в пределах анализируемого кардиокомплек-са, состоящего из N отсчетов (случай дискретного сигнала), отсчеты шума п(г), / = 1, 2, ..., N являются независимыми нормальными случайными величинами с нулевым средним значением и с

~ 2 /

дисперсией о (в пределах анализируемого кар-диокомплекса дисперсия помехи считается постоянной величиной, однако может изменяться от кардиокомплекса к кардиокомплексу).

2. Дрейф изоэлектрической линии в пределах кардиокомплекса к (/) описывается полиномом п-го

порядка к (г) = апгп + ... + с^г + а$. Справедливость данной аппроксимации подтверждается тем, что во многих работах (см., напр., [1], [3]) наиболее часто и достаточно эффективно использовалось кусочно-полиномиальное представление дрейфа изоэлек-трической линии. Ввиду небольшой продолжительности кардиокомплекса для аппроксимации дрейфа изоэлектрической линии в пределах анализируемого временного отрезка можно ограничиться использованием полинома 2-го порядка.

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

4. Индивидуальные временные зависимости длительностей зубцов и интервалов кардиоком-плекса от длительности интервала Я—Я установлены (например, на этапе предварительного обследования) и детерминированы [4], [5].

5. Форма всех зубцов кардиокомплекса установлена на стадии предварительного обследования, причем в силу неизменности механизма формирования эта форма сохраняется в интервале обследования (с изменением ритма изменяется только скорость проведения) [6], [7]. Математически условие неизменности формы состоит в том, что существует положительное не равное нулю число аг, такое, что sг (/) = аг/г (г), I = 1, 2, к, N, где sг (г) - зубец кардиокомплекса; /г (г) - эталонный сигнал для соответствующего зубца кардиокомплекса, причем г - индекс зубца кардиокомплекса.

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

Буквенное обозначение интервалов и индексов неудобно, поэтому индексу зубца присвоены следующие значения для типичного кардиокомплекса: г = 1 - интервал существования зубца Р; г = 2 - интервал существования зубца Q; г = 3 - интервал существования зубца Я; г = 4 - интервал существования зубца Б; г = 5 - интервал существования зубца Т. Неизменность формы зубцов дает возможность формально представить кардиокомплекс 5" (г) при нормальном состоянии пациента в виде суммы

5 (г ) = ^ аг/г (г ) =

г

= а1 /1 (г) + к + а5 /5 (г), (1)

где /г (г), г = 1, 2, .5 - функция, описывающая форму г-го зубца ЭКС (см. рисунок); а1, а2, к, а5 - весовые коэффициенты соответствующих зубцов модели.

Наблюдаемая реализация сигнала у (г) на кардиокомплексе представляет собой аддитивную смесь незашумленного кардиосигнала 5 (г), дрейфа изоэлектрической линии к (г) и миографической помехи п (г), причем к (г) и п (г) - совместно нормальные независимые случайные величины с одинаковой дисперсией и нулевым средним. Тогда реализация сигнала записывается следующим образом: у (г ) = 5 (г) + к (г) + п (г), где

к (г) = а2?'2 + сьг + сд.

В этом случае функция правдоподобия имеет вид

1

W [у| а, к (/)] =

(Т2я-а)

N

хехр

N

2о2

И у

И аг/г 0') + к 0)

г =1

(2)

где у = {у(1), у(2), к, у(N)} - вектор наблюдаемых отсчетов кардиокомплекса; а = {а1, а2, к, а5 } - вектор весовых коэффициентов.

У\ !

1 т ! У\ \ ! V 1 1

-V Г /2 —1- 1 1 1 1

/з 1 1 1 1 1 1 1 1 1

-VI/ . ^

/4

х

Г

\

2

1

/

V

0

0

Представим функцию правдоподобия в следующем виде:

1

W [у| а, к (г)] =

(л/2Р-о)

N

X

N

|2

I 1

X ехр {--- XIУ (г)_ к (г)]

I 2о г=1

I 1 N 5

Xехр { — у (г) - к (г)] X аг/г (г) -I о2 г=1

г=1

N

1 " X

2о г=1

5

X аг/г (г)

г=1

и введем обозначения: К1 = 1/ ^д/2Й-о) ; К2 =

= ехр {-¿^ У (г)-к (г )]2 |

Так как произведение Л( ) Л у ( ) = 0 при г Ф у,

N

имеем X г=1 Тогда

X агЛ (г)

Г=1

N 5

= XX а2 л2 (г).

г=1 г=1

W [у| а, к (г)] = К1К2 X

I 1 N 5

Xехр{-2X[У (г) - к(г)] X агЛг(г) -

.о2 г=1

Г=1

, N 5

1 ^ ^ ,22 Л,

XXа2Я (гН =

2о г=1 г=1

= К1К2 П ехр +

5

П

=1

2о

2

N

аг X[У (') - к (г)] Л (г)

г=1

о

(3)

N

где Ег = XЛ (г), г = 1, ..., 5 - энергия

эталон-

г=1

ного зубца кардиокомплекса.

С помощью метода максимального правдоподобия на основе полученной условной плотности вероятности, которая является статистической моделью сигнала, учитывающей дрейф изоэлек-трической линии и миографическую помеху, можно найти совместную оценку весовых коэффициентов амплитуд аь к, а5 и коэффициентов полинома а0 , а1, а2 , описывающего дрейф изоэлектрической линии. 50

При таком подходе возникает проблема, связанная с решением большого числа совместных уравнений: решение системы, состоящей более чем из пяти уравнений, приводит к катастрофическому росту объема вычислений и не позволяет получить решение задачи аналитически. Однако предварительные исследования [5], [8] показали, что при длительном наблюдении неизвестные весовые коэффициенты а могут интерпретироваться как независимые нормальные случайные величины с равным единице средним значением т{аг} = 1 и дисперсией Б {а г}. При этом считается, что математическое ожидание и дисперсия амплитуд зависят только от индекса зубца и не зависят от номера в последовательности наблюдаемых кардиоком-плексов. В силу независимости имеем: W(а) = W(а1, к, а5) = W(а1 )W(а2(а5),

где

X ехр

W (аг) = (1/^2яБ{аг }) X -(аг -1)7(2Б{аг})], г = 1, к, 5. (4)

При таком подходе для нахождения безусловной относительно вектора весовых коэффициентов плотности вероятности усредним функцию правдоподобия по распределению W (а) случайных параметров:

W [У№ (г)] =

¥ ¥

= | ... | W[уу |(а1, к, а5), ку (г)]X

—¥ —¥

X W (а1 (а5) й а1 — й а5.

(5)

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

В результате, подставив выражения (3) и (4) в (5), получим:

¥ ¥ 5

W [у| к (г)] = К1К2К3 | к | П ехр

а г ^Е г

N

агX[У(г)-к(г)]Лг (г)

г=1

-¥ г=1

а - 2 а

2

о

2Б {а г}

2о

й а1 — й а5,

X

2

2

+

+

+

N

где K3 = П(1Д/2—D{ar})exp(2D{ar})J . где zr = Sy(i) fr (i) - скалярное произведение

r=1 Так как

í exp (-px2 - qx) dx = — exp ——, p > 0 J \ p 4 p

[9], имеем:

W [У h ( i )] =

= K1K2K3 П I exp (-pra2 - qrar ) dar =

r =1 -¥

= KK2II exp [—r2/( 4 Pr r=1

где

k = K1K3 Пл/VPr ; r =1

' N N

qr =—7 S у (i) fr (i)-S h (i) f( i) .i=1 i=1

pr = Erl ( 2o2 ) +1/ ( 2D {ar }). С учетом (1) запишем:

w[y h (i)] =

D{ar }:

= K exp

N

S[ y (i )- h (i)]2 i=1

2o2

5 q2 11 4 pr

r=1

= K exp

N

Sy2 (' ) i=1

2o2

x

xexp

f

5

1

N [N ]2

-2zrSh(i) fr (i) + Sh (i) fr (i)

i=1 _i=1 ]

1 o4

N

2Sh(i) fr (i) i=1

o2D{ar }

x exp

r=1^pr

f „2

2z„

xexp

v o4 o2D{ar } D2 {ar } y

NN -2¿ y (i ) h (i) + ¿ h2 (i )

i=1_i=1

2o2

x

i=1

входной реализации и эталонного сигнала.

Из полученного выражения следует, что для нахождения оценок коэффициентов методом максимального правдоподобия необходимо найти такую совокупность значений âo, ^1, â2, при котором достигается минимум функции

I(âb â2)-

->min,

где

l = -

N n 2

2¿ y (i ) h (i ) + ¿ [h (i )] i=1

i=1

5

+S

r

=14 pr

2zr ¿ h (i ) fr (i )-

i=1

S h (i ) fr (i) i=1

o

+

2¿ h ( i ) fr (i ) i=1

D {ar }

(6)

Введем выражение

N

1 ¿[h (i )]2

2

i=1

222

= â2 A + â1 B1 + âo G1 + â2 â[Ci + â2 âo J1 + â1âo^1,

где

N

, N , N

A = 2 ¿ i4; в = 2 S i2; G = N2; Q = ¿ i3; i=1 i=1 i=1

NN

J=S i2; ^=S i,

и получим:

5

S-

=14 pr o

i=1

' N

S h (i ) fr (i ) i=1

i=1

= â2 A2 + c2 B2 +

+ âG + â2â1C2 + â2âo J2 + â^o^2,

где

5

A = S

r=1 5

4 pr o2

B2 = S

r=1

4 pr o2

'N

Si2 fr (i) i=1

' N

S ifr (i) ,i=1

00

+

2

2

1

+

1

2

Л

x

Л

2

1

1

+

2

1

^2 = X

г=1

5

^2 = X

г=1 5

4 Рг о2

N

N

X Л (г)

. г=1

N

1 N N

ЧXг2Л (г)ЪЛг (г)

2 Рг о г=1

О = X

г=1

N

1 N N

Ц" X г2Л (г )X Л (г)

г=1 [ 2 Рг о г=1 г =1

5 Г 1 N N |

и2 = X{т-LтXгЛг (г)XЛ (г)к

г=1 12Рго г=1 г=1 С учетом введенных обозначений выражение (6) запишем в следующем виде:

2 2 2 / = ?2С2 + ?1С1 + ?оСо + С2 + С1 В + Со +

+ С2 а^ + С2 Со+ С1Со^1 - С2^2 - С2В2 -

- Со ^2 - С2С1С2 - С2 Со О2 - С1Со^2 =

= с2 ( А - А2) + С2 (В - В2 ) + С2 ( - ^2 ) +

+ ?2С2 + + (оСо + С2С (С - С2 ) +

+ а2ао (О1 - О2 ) + а1ао (и1 - и2 ),

где

N

(2 =-X г2 У('

г=1

5

+ X

г=1

5

+X

г=1

5

+X

г=1

1 г

2 Рг

1

о

б{аг}) и

N

X г2Л (г)

N

{1 —X'' У(г) + =1

2 Рг

1

о

})

N

X г Лг (г)

N

(о =-X У(г) + =1

2 Рг

1

о

2

Б {аг } ) ~1

N

XЛ (г)

С учетом замены переменных А = А - А2;

В = В1 - В2; 6 = О1 - О2; С = С1 - С2; 0 = О1 - 02; и = и - и2 окончательно выражение (6) представим в виде

2 2 2 / = С2 А + С1 В + Со О + ?2С2 + +

+ ?оао + С2С1С + С2Со О + а\аои.

Для нахождения оценок коэффициентов, при которых функция (6) минимальна, необходимо решить следующую систему уравнений:

Э// Эа2 = 2с2 А + ?2 + сС + со О = о; Э// Эс1 = 2а^ + ?1 + С2С + Сои = о; Э// Эао = 2аоО + ?о + С2 О + с^и = о.

С помощью метода Крамера получены оценки коэффициентов:

С2 = ^2^1 + ^1^2 + ^3; С = ^2 + ^1^4 +

Со = ^2^3 + +

(7)

где

W1 =(и2 - 4ВО)/Д; W2 =(2СО - Ои)/Д; W3 =( 2ОВ - Си УД; W4 =(О2 - 4 АО )/д; W5 =( 2 Аи - ОС)/Д; W6 = (с 2 - 4 АВ )/д,

причем Д = 2(4АВО - Аи2 - С2О + СОи - О2В).

На основании (7) оценка дрейфа изоэлектри-ческой линии окончательно определяется как

к (г) = сс2 г2 + С^ + Со. Таким образом, найдена оценка дрейфа изоэлектрической линии (определены коэффициенты, входящие в полином аппроксимации). В выражения для оценок коэффициентов полинома входит дисперсия амплитуд зубцов Б{аг }.

Результаты предварительных исследований показали, что средний уровень амплитуды зубцов значительно меняется от человека к человеку, однако относительный разброс остается практически неизменным [5]. На основании данных исследований для практического применения можно использовать Б {аг } = о.о22.

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

1

X ехр

W [у| а, к (г)] =

1 N I

-Ь X{У ,г)-

гФ Р

N

П(^о)

г=1 гФ Р

[ 5

X агЛг (г) + к (г)

г=1

)

где Р - область отсчетов неизвестных фрагментов кардиокомплекса.

При этом дальнейшая последовательность действий для нахождения оценки дрейфа изоэлектри-ческой линии (описанная в статье) не изменяется.

2

1

г

г

1 г

г

(

1

г

г

(

\

СПИСОК ЛИТЕРАТУРЫ

1. Кардиомониторы. Аппаратура непрерывного контроля ЭКГ / А. Л. Барановский, А. Н. Калиниченко, Л. А. Манило и др.; под ред. А. Л. Барановского и А. П. Не-мирко. М.: Радио и связь, 1993. 248 с.

2. Пат RU 2436502 С2. МПК A61 B5/04, A61 B5/0402 (2006.01). Способ устранения дрейфа изоэлектриче-ской линии электрокардиосигнала и устройство для его осуществления / О. В. Мельник, А. А. Михеев, Н. С. Штрунова. Опубл. 20.12.2011. Бюл. № 35.

3. Siddiah N., Srikanth T., Kumar Y. S. Nonlinear filtering in ECG signal enhancement // Int. J. of computer science and communication networks. 2012. Vol. 2, № 1. P. 134-139.

4. Красичков А. С., Киреенков И. С., Нифонтов Е. М. Определение индивидуальной зависимости между временными параметрами электрокардиограммы // Вестн. аритмологии. 2004. № 35. C. 33.

A. S. Krasichkov Saint-Petersburg state electrotechnical university "LETI"

5. Красичков А. С. Анализ статистических закономерностей ЭКС // Биомедицинская радиоэлектроника. 2011. № 5. С. 18-23.

6. Красичков А. С., Соколова А. А. Оценка точности воспроизведения кардиосигнала в процессе синхронного накопления // Изв. вузов России. Радиоэлектроника. 2010. Вып. 6. С. 48-53.

7. Красичков А. С., Нифонтов Е. М., Иванов В. С. Алгоритм сортировки кардиокомплексов для анализа длительных записей электрокардиосигнала // Биомедицинская радиоэлектроника. 2011. № 11. С. 24-28.

8. Информационное обеспечение оперативной диагностики функционального состояния сердечно-сосудистой системы человека / О. М. Андреева, М. И. Бога-чев, А. С. Красичков и др. СПб.: Изд-во СПбГЭТУ "ЛЭТИ", 2009. 120 с.

9. Прудников А. П., Брычков Ю. А., Маричев О. И. Интегралы и ряды. М.: Наука, 1981. 797 с.

Algorithm for ECG base line drift estimation in long-term monitoring records

Algorithm to estimate baseline drift in long-term ECG monitoring records under hindering noise conditions is considered. Electro cardio signal, myographic noise, base line drift, maximum likelihood methods, long-term monitoring Статья поступила в редакцию 5 июля 2013 г.