Перспективы применения преобразования Гильберта-Хуанга для автоматизированной обработки электрокардиосигналов Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 53.088.7, 612.172.4

Л. Ю. Кривоногов, канд. техн. наук, А. Ю. Тычков, аспирант,

ГОУ ВПО «Пензенский государственный университет»

Перспективы применения преобразования Гильберта—Хуанга для автоматизированной обработки электрокардиосигналов

Ключевые слова: электрокардиосигнал, преобразование Гильберта—Хуанга, подавление помех. Key words: elecrocardiography signal, Hilbert—Huang transform, denoising.

В статье обосновано применение адаптивных базисных функций для анализа электрокардиосигналов; рассмотрено преобразование Гильберта—Хуанга — новый метод анализа нелинейных и нестационарных процессов, описан алгоритм подавления помех в электрокардиосигналах, основанный на их разложении по эмпирическим модам, анализе и обработке отдельных мод и последующей реконструкции сигналов; приведены результаты подавления различных помех; разработаны виртуальные приборы в среде графического программирования ЬаЪУ1ЕЖ.

Введение

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

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

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

Обоснование применения адаптивного базиса для анализа ЭКС

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

Традиционно для анализа ЭКС применяются амплитудно-временные методы. На сегодняшний день выявлены и описаны связи амплитудных и временных параметров ЭКС с различными состояниями сердечно-сосудистой и других функциональных систем организма. Однако в последнее время специалисты стараются принимать решения о патологиях сердеч-

но-сосудистой системы с учетом результатов анализа ЭКС в частотной области.

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

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

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

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

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

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

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

• качественное подавление помех на основе разделения в пространстве признаков информации о помехах и сигнале с последующим использованием для восстановления только информации о сигнале;

• разделение различных процессов (например, ЭКС плода на фоне ЭКС матери);

• эффективное обнаружение и классификацию информативных участков сигнала (комплексов, зубцов, сегментов) и измерение их параметров с высокой точностью;

• реализацию удобного для дальнейшего анализа преобразования, в результате которого сигнал представляется в частотно-временной области (в виде ЗБ поверхности в системе координат «амплитуда — частота — время» или «энергия — частота — время»), в связи с чем появится возможность выявления скрытых модуляций и областей концентрации энергии, а следовательно, и новых диагностических признаков в ЭКС.

Преобразование Гильберта—Хуан га

Для формирования адаптивного базиса, отвечающего реальным изменениям сигнала во времени, наибольший практический интерес представляют точки экстремумов, разрывов, перегибов, нарушения монотонности. В 1995 году Н. Хуанг предложил метод анализа нелинейных и нестационарных сигналов, который впоследствии был назван преобразованием Гильберта—Хуанга [З]. Предложенное преобразование включает два основных этапа:

• разложение сигнала на компоненты, или декомпозиция на эмпирические моды (ДЭМ);

• формирование по полученным эмпирическим модам (ЭМ) спектра Гильберта.

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

• локальность, то есть возможность учета локальных особенностей сигнала;

• ортогональность, обеспечивающая восстановление сигнала с определенной точностью;

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

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

• общее число экстремумов равняется общему числу нулей с точностью до единицы;

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

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

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

большого объема вычислений и предполагает выполнение следующих действий [З]:

• определение локальных экстремумов (максимумов и минимумов) входного сигнала /(п);

• вычисление верхней огибающей и^(п) по найденным локальным максимумам и нижней огибающей ¿у(п) по локальным минимумам; огибающие строят чаще всего с помощью кубической сплайн-интерполяции;

• вычисление среднего значения огибающих (локального тренда): т^(п) = (и^(п) — й^{п))/2.

• получение приближения к ЭМ (вычитание локального тренда из результата, полученного на предыдущем шаге декомпозиции): gíj + ^(п) = gij(п) — - Щ(п);

• повторение предыдущих действий для выделенных приближений к ЭМ gij(п) (так называемая процедура отсеивания) и получение новых приближений и новых локальных трендов т^(п). Итерации продолжаются до тех пор, пока приближения не станут удовлетворять критерию останова (порогу по среднеквадратической разности между двумя результатами последовательных операций приближения). При выполнении критерия останова итерации прекращают и считают приближение gij(п) эмпирической модой Д(п) = gi^•(п);

• расчет остатка — разности между сигналом и ЭМ: г1 +\(п) = г^п) — Д(п). Первые пять действий приводят с остатком г + \(п).

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

Г-1

Г (л) = X П (л)+ 'V

1=1

При этом каждая функция является монокомпонентной, имеет свой характерный временной масштаб осцилляций, который убывает с увеличе-

остаток ги(п)

Рис. 1\ Упрощенная схема алгоритма декомпозиции сигнала

а) и, мВ

в) и, мВ 200 150

д) и, мВ

'-) и, мВ 150

к) и, мВ з -,

2 1

0 -1

-2 -3 -4

л) и, мВ 0

-5 -10

-15 -20 -25

п

Рис. 2

Пример декомпозиции фрагмента ЭКС на эмпирические моды: а — исходный ЭКС; б—к — эмпирические моды; л — глобальный тренд

нием ее номера ¿, а процесс вычисления сводится к устранению локального тренда, соответствующего данному масштабу. На рис. 2 приведен фрагмент ЭКС и результат его декомпозиции на ЭМ.

На втором этапе преобразования Гильберта—Ху-анга по полученным ЭМ формируется частотно-временное представление сигнала — спектр Гильберта. Для каждой ЭМ рассчитывается мгновенная частота с помощью преобразования Гильберта, затем мгновенные частоты нормируются и наносятся на ЗБ-диаграмму. Информация об амплитуде или энергии отображается соответствующим цветом. Спектр Гильберта позволяет определить наличие и характер амплитудных и частотных модуляций в сигнале, идентифицирует временные и частотные диапазоны концентрации энергии. Используя выражение, связывающее мгновенную амплитуду, мгновенную частоту и время (номер отсчета), можно построить трехмерную поверхность в системе координат « амплитуда — частота — время» или «энергия — частота — время». Представление сигнала в виде ЗБ-поверх-ности в системе координат «энергия — частота — время» характеризует распределение мгновенной энергии в каждой точке частотно-временной плоскости. На рис. З показана поверхность в системе координат «энергия — частота — время» для фрагмента ЭКС (см. рис. 2, а).

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

Применение преобразования Гильберта— Хуанга для подавления помех в ЭКС

Рассмотрим подробнее вопросы применения технологии преобразования Гильберта—Хуанга для подавления помех на этапе предварительной обработки ЭКС. Известно, что основы качественного автоматизированного анализа ЭКС закладываются на этапах регистрации и предварительной обработки [4]. При этом главной причиной ошибок автоматической диагностики в электрокардиологии являются погрешности, допущенные на этапе измерения амплитудно-временных параметров ЭКС [5]. Сами погрешности измерений напрямую связаны с наличием помех в регистрируемом сигнале и искажений, полученных в результате некачественной регистрации и предварительной обработки.

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

Традиционно различают следующие виды проявления помех в электрокардиографии:

• дрейф изолинии (низкочастотные колебания с частотой менее 0,5 Гц);

• сетевая помеха (суперпозиция гармоник разных фаз с частотами, кратными 50 Гц);

Рис. 3 | Поверхность в системе координат «энергия — частота — время»

Рис. 4\ Структура, реализующая процедуру подавления помех в ЭКС на основе ДЭМ

• мышечный тремор (хаотические колебания с частотой З0—100 Гц);

• артефакты движения (одиночные или циклические волны с частотой от единиц до З0—40 Гц);

• высокочастотные шумы электродов и усилителей;

• импульсные помехи.

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

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

Для эффективного подавления помех в ЭКС необходимо провести анализ каждой ЭМ в целях их классификации на три группы:

• сигнальные, содержащие только полезную информацию;

• помеховые, содержащие только помеховую составляющую;

• смешанные, в которых присутствует полезная информация вместе с частотными компонентами помех.

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

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

процедур. Наиболее распространенные методы удаления помех из ЭМ — нелинейная пороговая обработка, фильтрация Савицкого—Галея и линейная частотная. Исследовав алгоритмы [6—8], мы предлагаем алгоритм подавления помех в ЭКС.

Алгоритм подавления помех в ЭКС

Подавление высокочастотных помех (сетевой помехи и мышечного тремора). Исследуемый ЭКС раскладывается на ЭМ, затем несколько высокочастотных ЭМ подвергаются пороговой обработке в соответствии с выражением:

Л(f (п)) -

_ isign(f (п)) х(|f (п)| -Pj ), if If.(п)| > Pj 1 "(о, if I f (п )| < P J'

где fj(n) — отсчеты конкретной j-й ЭМ; pj — порог, рассчитанный на основе медианной оценки уровня помех в соответствующей ЭМ, pj = Rj[(median | fj(n) -- median (fj(n))|)/0,6745], Rj — индивидуальный коэффициент для каждой ЭМ.

Результаты подавления тремора мышц предложенным алгоритмом приведены на рис. 5.

а) и,мВ 500

L t HI AJ 1

.Л; /181 kr^ejl^ 1/э01 Ьэ1 1081 уЦ|71 i*621

Рис. 5

Результат подавления тремора мышц в ЭКС: а — ЭКС с помехой; б — ЭКС после процедуры помехо-подавления

Биомедицинская инженерия

а) и,мВ 500 п

Рис. 6

Результат подавления дрейфа изолинии в ЭКС: а — ЭКС с помехой; б — ЭКС после процедуры помехоподавления

Подавление низкочастотных помех (дрейфа изолинии и артефактов движения). Исследуемый ЭКС подвергается декомпозиции на эмпирические моды, глобальный тренд удаляется, несколько самых низкочастотных ЭМ фильтруются с помощью фильтра Баттерворта второго порядка (частота среза — 0,2 Гц) для подавления низкочастотных составляющих этих мод. На рис. 6 приведены результаты подавления дрейфа изолинии предложенным алгоритмом.

Результаты моделирования

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

В =

0,5

Е (ъ

Е Ъ

¿=1

100%.

Проведенные исследования показали, что значения среднеквадратической ошибки для большинства реальных помех не превышают 2,5-3,0 % и лишь для некоторых комбинированных помех достигают 4,0-5,0 % (результаты получены без учета краевых участков сигнала, на которых искажения больше).

Алгоритм реализован в среде графического программирования LabVIEW. На рис. 7 показана блок-диаграмма виртуального прибора LabVIEW, реализующего процедуру подавления высокочастотных помех в ЭКС.

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

Блок ввода данных

Рис. 7\ Реализация процедуры подавления высокочастотных помех в ЭКС (виртуальный прибор LаЪVIEW)

Вывод

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

Визуализация ЭКС в виде поверхности энергетической плотности частотных составляющих (ЗБ-по-верхности в системе координат «энергия — частота — время») даст возможность врачам выявить новые диагностические признаки, невидимые при амплитудно-временном представлении сигнала.

I Литература I

1. Дроздов Д. В., Леванов В. М. Автоматический анализ ЭКГ: проблемы и перспективы // Здравоохранение и медицинская техника. 2004. № 1. Режим доступа: Ьйр://

www.altonika.ru/article.php?id=187. Дата обращения: 05.03.2010.

2. Ватанабе С. Разложение Карунена—Лоэва и факторный анализ. Теория и приложения // Автоматический анализ сложных изображений / Под ред. Э. М. Браверман-на. М.: Мир, 1969. С. 163-181.

3. Huang N. Е., Shen Z., Long S. R. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis // Proc. Royal Soc. Lond. 1998. Vol. 454. P. 903-995.

4. Злочевский M. С. Критерии и особенности качества средств автоматического анализа электрокардиосигналов // Приборы, средства автоматизации и системы управления. ТС-10 «Медицинские приборы, оборудование и инструменты». М.: Информприбор, 1988. Вып. 7.

5. Чирейкин Л. В., Шурыгин Д. Я., Лабутин В. К. Автоматический анализ электрокардиограмм. Л.: Наука, 1977. 248 с.

6. Pat. № 20080269628 USA. А61В5/0402; А61В5/0402. De-

noising and artifact rejection for cardiac signal in a sensis system / D. W. Koertige, H. Zhang, В. Pelzek. Assignee: Siemens Medical Solutions USA. Inc. N 11/831143; filing date 31.07.2007; publ. date 30.10.2008.

7. Blanco-Velasco M., Weng В., Barner К. Е. A New ECG Enhancement Algorithm for Stress ECG Tests // Computers in Cardiology. 2006. Vol. 33. P. 917-920.

8. Na Pan, Vai Mang, Mai Peng Un, Pun Sio hang. Accurate removal of baseline wander in ECG using empirical mode decomposition // Noninvasive Functional Source Imaging of the Brain and Heart and the International Conference on Functional Biomedical Imaging, 2007. NFSI-ICFBI 2007. Joint Meeting of the 6th International Symposium on. Hangzhou, China, October 12-14. Hangzhou, 2007. P. 177-180.