Обзор алгоритмов обнаружения QRS-комплекса.

В литературе мне удалось найти 2 метода обнаружения QRS-комплекса: 1) методы, основанные на производной; 2) алгоритм Пана-Томпкинса.

Метод обнаружения QRS-комплекса, основанный на производной, имеет 2 различных решения.

Решение 1. Отметим, что для QRS-комплекса характерна наибольшая крутизна наклона сигнала (скорость изменения напряжения) в сердечном цикле, в связи с тем что желудочки обладают свойствами быстрого проведения и деполяризации. Поскольку скорость изменения задаётся оператором производной, то операция d/dt является наиболее логичной отправной точкой в попытке разработать алгоритм для обнаружения QRS-комплекса.

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

Рис. 4. Результат фильтрации сигнала ЭКГ с низкочастотной помехой с использованием оператора трехточечной центральной разности.

Baida и др. предложили основанный на производной алгоритм выделения QRS-комплекса, который позже был исследован и оценен Alhstrom и Tompkins, Friesen и др. и Tompkins. Алгоритм работает следующим образом. Аппроксимируется сглаженная трёхточечная первая производная у_о(п) от данного сигнала х(п):

у_о(п) = | х(п) - х(п- 2)|. (1)

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

y₁(п) = | х(п) — 2 х(п — 2) + х(п — 4)|. (2)

Эти два результата взвешиваются и комбинируются, что даёт

y₂(п) = 1,З у_о(п) + 1,1 y₁(п). (3)

Результирующий сигнал y₂(п) сканируется с использованием порога равного 1,0. Во всех случаях, когда этот порог пересекается с сигналом, последующие восемь отсчётов также проверяются с использованием того же самого порога. Если по меньшей мере шесть из восьми точек проходят тест с использованием порога, данный сегмент, состоящий из восьми отсчётов, считается частью QRS-комплекса. В результате этой процедуры получается импульс, имеющий ширину, пропорциональную ширине QRS-комплекса; однако этот метод весьма чувствителен к шуму.

Иллюстрация применения. Верхняя кривая на рис. 5 иллюстрирует два цикла Фильтрованного сигнала ЭКГ. Сигнал был пропущен через фильтр нижних частот Баттерворта восьмого порядка f_с = 90 Гц, прорежен в пять раз и профильтрован с использованием режекторного фильтра fо = 60 Гц. Эффективная частота дискретизации составляет 200 Гц. Сигнал был нормализован делением на его максимальную величину.

Рис. 5. Сверху вниз: два цикла фильтрованного сигнала ЭКГ; выходной сигнал уо(п) основанного на первой производной оператора, заданного уравнением (1); выходной сигнал у₁(п) основанного на второй производной оператора, заданного уравнением (2); комбинированный результирующий сигнал y₂(п), полученный после использования уравнения (3); результирующий сигнал y₃(п).

Вторая и третья кривая на рис. 5 показывают производные уо(п) и у₁(п) сответственно; четвёртая кривая иллюстрирует комбинированный результирующий сигнал y₂(п). Стоит обратить внимание на относительно высокие величины результирующего сигнала после применения производных в зоне расположения QRS-комплекса; выходной сигнал значительно ниже или пренебрежимо мал в местах расположения зубцов Р и Т, несмотря на то что в исходном сигнале имелись необычно крутые и высокие Т-зубцы. Также видно, что результирующий сигнал имеет многочисленные пики на всём протяжении QRS-комплекса, что связано с наличием у QRS-комплекса в данном примере трех основных склонов: Q-R, R-S и S-ST (в других сигналах ЭКГ может также присутствовать дополнительный склон PQ-Q).

Последняя кривая на рис. 5 показывает результат сглаживания уз(п), полученный после прохождения сигнала y₂(п) через фильтр 8-точечного скользящего среднего. В результате мы получили единственный импульс с амплитудой, превышающей 1,0 на всем интервале регистрации соответствующего QRS-комплекса. Для обнаружения QRS-комплекса может быть использован простой алгоритм поиска пиков. Для определения истинной локализации соответствующего QRS-комплекса необходимо вычесть из найденных точек локализации QRS-комплекса суммарную задержку, вносимую фильтрами.

Отметим, что процедура поиска пика не может быть выполнена непосредственно по сигналу ЭКГ: QRS-комплекс может быть не самой высокой по амплитуде волной в кардиоцикле, а имеющиеся артефакты могут легко сделать процедуру поиска неэффективной. Обратите внимание, что сигнал ЭКГ в данном примере перед вычислением производных был профильтрован для ограничения полосы частот в полосе до 90 Гц и что он был свободен от дрейфа базовой линии^[3].

Решение 2. Murthy и Rangaraj предложили алгоритм обнаружения QRS- комплекса, основанный на взвешенном и возведённом в квадрат операторе первой производной и на фильтре скользящего среднего. В этом случае основанный на производной оператор определяется следующим образом:

g₁(n) = , (4)

где х(п) — сигнал ЭКГ, N — ширина окна, в пределах которого разность первого порядка вычисляется, возводится в квадрат и взвешивается с использованием коэффициента (N — i + 1). Взвешивающий коэффициент линейно спадает, начиная от текущей разности до разности, вычисленной на N отсчётов ранее по времени, что обеспечивает сглаживающий эффект. Дальнейшее сглаживание результата выполняется с использованием фильтра скользящего среднего по М точкам:

g(n) = (5)

При частоте дискретизации 100 Гц ширина окна фильтра устанавливается как M= N = 8. Данный алгоритм даёт единственный пик для каждого QRS-комплекса и подавляет зубцы Р и Т.

Поиск пика в обработанном сигнале g(n) может быть выполнен с помощью простого алгоритма поиска пиков следующим образом.

1. Сканирование фрагмента сигнала g(n), на котором ожидается наличие пика, и определение максимальной величины g_max. Максимальное значение g(n) на всём протяжении анализируемого участка может быть также принято равным g_max.

2. Определение порога как некоторой доли от максимума, например, Th = 0,5g_max.

3. Для всех g(n) > Th выбираются те отсчёты, для которых соответствующие величины g(n) больше, чем определённое заданное число М предыдущих или последующих отсчётов g(n), т. е.

{p} = [n | g(n) > Th] И

[g(n) > g(n - i), i = 1,2, ...,M] И (6)

[g(n) > g(n + i), i = 1,2, ...,M].

Определённый таким образом набор {р} содержит индексы всех найденных в сигнале g(n) пиков.

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

Иллюстрация применения. Верхняя кривая на рис. 6 иллюстрирует два цикла фильтрованного сигнала ЭКГ. Фильтрация сигнала выполнялась с использованием фильтра нижних частот Баттерворта 8-го порядка с f_с = 40 Гц. Далее сигнал был прорежен в 10 раз, в результате чего эффективная частота дискретизации составила 100 Гц, что соответствует параметрам, использованным Murthy и Rangaraj. Сигнал был нормализован делением на его максимальную величину.

Рис. 6. Сверху вниз: два цикла фильтрованного сигнала ЭКГ; выходной сигнал g₁(n) оператора взвешивания и возведения в квадрат первой производной, заданного уравнением (4); выходной сигнал g(n) сглаживающего фильтра, заданного уравнением (5).

Вторая и третья кривые на рис. 6 показывают выходные сигналы основанного на производной оператора и сглаживающего фильтра. Обратите внимание, что окончательный выходной сигнал содержит единственный гладкий пик для каждого QRS-комплекса и что зубцы Р и Т не дают никаких существенных значений выходного сигнала. Для детектирования и сегментации каждого комплекса может быть использован простой алгоритм поиска пиков^[3].

Алгоритм обнаружения QRS-комплекса Пана Томпкинса. Пан и Томпкинс (Pan и Tompkins) предложили алгоритм обнаружения QRS-комплексов, рассчитанный на работу в реальном времени и основанный на анализе наклона, амплитуды и ширины QRS-комплексов. Этот алгоритм состоит из следующей последовательности фильтров и методов: фильтр нижних частот, фильтр верхних частот, оператор производной, возведение в квадрат, интегрирование, адаптивная пороговая процедура и процедура поиска. На рис. 7 в схематической форме показаны основные этапы этого алгоритма.

Рис. 7. Блок-диаграмма алгоритма обнаружения QRS-комплекса Пана-Томпкинса.

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

(7)

Выходной сигнал у(п) связан с входным сигналом х(п) разностным уравнением:

у(п) = 2 у(n - 1) - у(п -2) + 1/32[x(n) - 2x(n - 6) + х(п -12)]. (8)

При частоте дискретизации 200 Гц этот фильтр имеет довольно низкую частоту среза f_с=11 Гц и вносит задержку, составляющую 5 отсчётов или 25 мс. Этот фильтр имеет ослабление более 35 дБ на частоте 60 Гц и эффективно подавляет сетевую наводку, если она существует.

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

(9)

отношение между входом и выходом задаётся разностным уравнением:

у(п) = у(п — 1) + х(п) — х(п — 32). (10)

Передаточная функция H_hp(z) фильтра верхних частот задаётся следующим образом:

H_hp(z) = (11)

Аналогично, выходной сигнал р(п) фильтра верхних частот задаётся разностным равнением:

р(п) = р(п - 1) – 1/32x(n) + x(n-16) – x(n-17) + 1/32x(n-32), (12)

где х(п) и у(п) связаны между собой в соответствии с уравнением (10). Этот фильтр верхних частот имеет частоту среза 5 Гц и вносит задержку 80 мс.

Оператор производной. Операция дифференцирования, использованная Паном Томпкинсом, задаётся следующим образом:

у(п)=1/8 [2х(п) + х(п — 1) — х(п — 3) — 2х(п — 4)], (13)

и при частотах до 30 Гц аппроксимирует идеальный оператор d/dt. Эта процедура взятия производной подавляет низкочастотные компоненты зубцов Р и Т и имеет высокий коэффициент усиления для высокочастотных компонент, появляющихся из-за крутых склонов QRS-комплекса.

Возведение в квадрат. Операция возведения в квадрат делает результат положительным и усиливает большие разности, возникающие из-за QRS-комплексов; маленькие разности, возникающие на Р- и Т-зубцах, при этом подавляются. Высокочастотные компоненты в сигнале, связанные с QRS-комплексом, ещё более усиливаются.

Интегрирование. Как можно было наблюдать в предыдущем разделе, выходной сигнал операции, основанной на производной, демонстрирует многочисленные пики в пределах длительности отдельного QRS-комплекса. Алгоритм Пана-Томпкинса выполняет сглаживание выходного сигнала предыдущих операций с использованием интегрирующего фильтра типа скользящего окна, задаваемого уравнением

у(п) =1/N [х(п - (N - 1)) + х(п — (N — 2)) + ... + х(п)]. (14)

Выбор ширины окна N должен выполняться с использованием следующих соображений: при выборе слишком большой ширины выходные сигналы, связанные с QRS-комплексом и Т-зубцом будут сливаться, в то время, как слишком маленькая ширина приведёт к нескольким пикам для единственного QRS-комплекса. Для частоты дискретизации f_s = 200 Гц была выбрана величина N = 30. На рис. 8 показано влияние ширины окна на выходной сигнал интегратора и его связь с шириной QRS-комплекса.

Рис. 8. Взаимосвязь QRS-комплекса и выходного сигнала интегратора, основанного на скользящем окне: схематичное изображение сигнала ЭКГ (выше); выходной сигнал интегратора, основанного на скользящем окне (ниже). QS — ширина QRS-комплекса, W — ширина окна интегратора, заданная как N/f_s с.

Адаптивная пороговая процедура. Пороговая процедура в алгоритме Пана-Томпкинса адаптируется к изменениям в сигнале ЭКГ путём вычисления скользящих оценок пиков, связанных с сигналом и шумом. Считается, что пик обнаружен каждый раз, когда выходной сигнал изменяет своё направление в пределах заданного интервала. В ходе дальнейшего рассмотрения SPKI будет означать уровень пика, который на этапе обучения алгоритма был интерпретирован как пик, соответствующий QRS, a NPKI будет обозначать уровень пика, связанного с событиями, не являющимися QRS-комплексами (шум, ЭМГ и т. д.). Величины THRESHOLD II и THRESHOLD I2 — это два значения порогов, которые используются для отнесения обнаруженных пиков к одной из двух категорий: сигналу (QRS-комплекс) или шуму.

Каждый новый обнаруженный пик относится либо к категории «пик сигнала», либо к категории «пик шума». Если пик превышает THRESHOLD I1 во время первого шага анализа, он классифицируется как пик QRS-комплекса (сигнала). Если используется процедура обратного поиска (описанная в следующем параграфе), то данный пик, для того чтобы быть названным QRS-комплексом, должен быть выше THRESHOLD I2. После того как обнаружен каждый пик, уровни пиков и порогов обновляются и классифицируются следующим образом:

SPKI = 0,125РEAKI + 0,875SPKI, если PEAKI-пик сигнала;

NPKI = 0,125PEAKI + 0,875NPKI, если PEAKI-пик шума; (15)

THRESHOLD I1 = NPKI + 0,25(SPKI - NPKI);

THRESHOLD I2 = 0,5THRESHOLD I1. (16)

Для случая процедуры обратного поиска с использованием величины THRESHOLD I2 формула обновления величины SPKI в случае обнаружения QRS-комплекса изменяется на следующую:

SPKI = 0,25РЕАКI + 0,75 SPKI. (17)

Процедура обратного поиска. Алгоритм Пана-Томпкинса постоянно обновляет две средние величины RR-интервала: RR AVERAGE1 — это среднее для восьми последних комплексов и RR AVERAGE2 — это среднее для восьми последних комплексов, имеющих RR-интервалы в диапазоне, заданном пределами RR LOW LIMIT = 0,92 х RR AVERAGE2 и RR HIGH LIMIT = 1,16 x RR AVERAGE2. Каждый раз, когда QRS-комплекс не обнаружен на определённом интервале, заданном как RR MISSED LIMIT = 1,66 х RR AVERAGE2, за QRS-комплекс принимается пик, амплитуда которого лежит в диапазоне между порогами, установленными для процедуры обратного поиска.

Данный алгоритм продемонстрировал очень низкую степень ошибки, составляющую 0,68%, или 32 комплекса в час на базе данных, содержащей примерно 116 000 QRS-комплексов, полученных из 24-часовых записей ЭКГ от 48 пациентов.

Иллюстрация применения. Верхняя кривая на рис. 9 иллюстрирует тот же самый сигнал, что и на рис. 8. На том же самом рисунке показаны выходные сигналы различных этапов алгоритма. Подробное рассмотрение этого рисунка приводит примерно к тем же выводам, что и для метода, основанного на производной. Оператор производной подавляет Р и Т зубцы и даёт большой выходной сигнал на участках сигнала, где находятся QRS-комплексы. Операция возведения в квадрат преимущественно увеличивает большие значения и усиливает высокочастотные компоненты.

Рис. 9. Результаты работы алгоритма Пана-Томпкинса. Сверху вниз: два цикла фильтрованного сигнала ЭКГ (тот же самый сигнал, что показан на рис. 8); выходной сигнал полосового фильтра (комбинация низкочастотного и высокочастотного фильтров); выход оператора, основанного на производной; результат возведения в квадрат; результат интегрирования, умноженный на 100.

Результирующий сигнал содержит многочисленные пики для каждого QRS-комплекса и, следовательно, должен быть сглажен. Окончательным выходным сигналом интегратора является единственный гладкий импульс для каждого QRS-комплекса. Обратим внимание на сдвиг между истинным положением QRS-комплекса и выходным импульсом, вызванный суммарными задержками различных фильтров^[3].

Не нашли, что искали? Воспользуйтесь поиском по сайту: