CC BY
132
УДК 53.088.7, 612.172.4
Л. Ю. Кривоногое, Д. В. Папшев
ПОВЫШЕНИЕ ЭФФЕКТИВНОСТИ ПОДАВЛЕНИЯ ВЫСОКОЧАСТОТНЫХ ПОМЕХ В ЭЛЕКТРОКАРДИОСИГНАААХ
L. Yu. Krivonogov, D. V. Papshev
EFFECTIVE ENHANCEMENT OF HIGH-FREQUENCY INTERFERENCE SUPPRESSION IN ECG SIGNALS
Аннотация. Рассмотрены вопросы подавления помех в электрокардиосигналах (ЭКС) на основе эмпирической модовой декомпозиции (EMD). Обосновано применение адаптивного базиса для анализа ЭКС. Описан алгоритм эмпирической модовой декомпозиции. Предложен усеченный алгоритм EMD. Рассмотрены задачи нелинейной пороговой обработки эмпирических мод ЭКС. Ключевой параметр нелинейной пороговой обработки - определение значения порога дискриминации, от него зависит эффективность подавления помех в ЭКС. Рассмотрены известные робастные оценки масштаба. Выбрана Q-оценка для вычисления параметра уровня шума в отдельных эмпирических модах. Разработан алгоритм определения порога дискриминации при подавлении высокочастотных помех в ЭКС на основе усеченной EMD.
Abstract. In this article questions of ECG signals denoising based on empirical modes decomposition are considered. The application of adaptive basis for the analysis of ECG. The algorithm of empirical mode decomposition algorithm are considered. The truncated algorithm of empirical modes decomposition algorithm proposed. The tasks of nonlinear thresholding empirical modes shown. The key parameter in the thresholding process is the choice of the threshold value. The known robust scale estimators are considered. Q-estimator is selected for the calculation of noise parameter in some empirical modes. An algorithm for determining the threshold value in suppressing high-frequency interference in ECG signals based on truncated EMD.
Ключевые слова: электрокардиосигналы, подавление помех, разложение на эмпирические моды, нелинейная пороговая обработка, робастная оценка масштаба.
Key words: ECG signals, interference suppression, empirical mode decomposition, nonlinear thresholding, robust scale estimators.
Введение
По данным Всемирной организации здравоохранения, сердечно-сосудистые заболевания являются основной причиной смерти населения во многих экономически развитых странах, в том числе и в России. Повышение эффективности лечения и возвращение пациентов к активной жизни связаны, прежде всего, со своевременной диагностикой, основой которой является электрокардиография. С начала 1960-х гг. в электрокардиографию начинают активно внедряться методы автоматического анализа электрокардиосигналов. В настоящее время автома-
тический анализ ЭКС дает возможность внедрения новейших достижений мировой кардиологии в широкую медицинскую практику, что достигается тиражированием программных продуктов ведущих мировых исследовательских центров для компьютерных кардиоанализаторов. Несмотря на очевидные успехи, достигнутые за полувековую историю автоматизации электрокардиографических исследований, эффективность компьютерного анализа ЭКС еще недостаточно высока. Ошибки и неточности автоматической диагностики имеют целый ряд причин, одной из которых является искажение полезного сигнала помехами или процедурами первичной обработки. Если в результате подавления помех происходят искажения формы информативных участков ЭКС, то это может привести к ошибочным или неточным диагностическим заключениям. В этой связи актуальна разработка способов и алгоритмов фильтрации ЭКС, обеспечивающих значительное подавление помех при минимальном искажении формы полезного сигнала [1].
ЭКС представляют собой нестационарный сложноструктурированный сигнал с циклически повторяющимися локально сосредоточенными информативными участками. При формировании, регистрации и обработке ЭКС подвергаются воздействию различных помех, часть из которых представляет собой случайные аддитивные, высокочастотные (по отношению к ЭКС) шумы, вызванные сокращением скелетных мышц, работой электродов и усилителей и некоторыми другими причинами.
Обоснование применения адаптивного базиса для анализа ЭКС
До последнего времени в цифровой обработке сигналов в основном использовались методы линейной частотной фильтрации, что связано с наличием подходящего математического аппарата, простотой интерпретации и расчета линейных фильтров. В то же время применение линейных фильтров не обеспечивает высокое качество подавления помех в ЭКС. Так как спектры полезного сигнала и помех в общем случае перекрываются, применение линейных фильтров приводит к нежелательному искажению ЭКС. При использовании линейных фильтров нижних частот (ФНЧ) для устранения высокочастотных (ВЧ) помех острые зубцы (<2, Я, 5) и другие высокочастотные компоненты ЭКС сглаживаются. Такое искажение формы недопустимо, так как высокочастотные компоненты ЭКС несут диагностическую информацию. При расширении полосы пропускания ФНЧ частотный спектр части помех оказывается в полосе пропускания фильтра, и, соответственно, такие помехи не будут подавлены. Таким образом, применение линейных ФНЧ для подавления высокочастотных помех в ЭКС недостаточно эффективно.
Для повышения эффективности подавления помех в нестационарных сигналах в последние годы широкое применение нашли методы разложения сигналов на узкополосные составляющие по локально сосредоточенным базисам. При этом полагают, что при правильно подобранном базисе узкополосная помеха может занять всю частотную составляющую и тогда ее можно будет исключить. Если же помеха широкополосная, то она проявляется сразу в нескольких субполосах малыми шумами, которые удаляют с помощью специальной обработки.
В качестве локально-сосредоточенных базисов наибольшее распространение получили вейвлеты. Нестационарный сигнал анализируется путем разложения по базисным функциям, полученным из некоторого прототипа на основе его масштабных растяжений и сдвигов вдоль временной оси. Функция-прототип называется базисным вейвлетом. Фильтры на основе вейвлет-преобразования (ВП) стремительно завоевывают популярность для анализа нестационарных сигналов благодаря хорошо разработанному математическому аппарату, наличию множества специально разработанных вейвлетов, быстрым вычислительным алгоритмам, возможности решать широкий класс задач. Основная проблема, сдерживающая эффективное практическое применение ВП, в частности для анализа биомедицинских сигналов - большое многообразие и неочевидность выбора базисного вейвлета для решения конкретной задачи. Для корректного анализа биомедицинских сигналов необходимо формирование не априорно заданного, а адаптивного базиса, функционально зависимого от содержания самого сигнала (т.е. учитывающего структуру сигнала, изменчивость его параметров, наличие помех различного вида и интенсивности).
Преобразование Гильберта-Хуанга и эмпирическая модовая декомпозиция
Для адаптивного анализа сигналов наибольший практический интерес представляют точки экстремумов, разрывов, перегибов, нарушения монотонности, отвечающие реальным изменениям сигнала во времени. Именно такие точки (локальные экстремумы) были использованы Норденом Хуангом для формирования адаптивного базиса нового метода анализа нелинейных и нестационарных сигналов, который впоследствии был назван преобразованием Гильберта-Хуанга (Hilbert-Huang Transform, HHT) [2]. Предложенное преобразование включает два основных этапа:
- разложение сигнала на частотные компоненты - эмпирическая модовая декомпозиция (Empirical Mode Decomposition, EMD);
- формирование по полученным эмпирическим модам (ЭМ) спектра Гильберта (Hilbert Spectral Analysis, HSA).
Для задачи подавления помех в сигналах используется EMD. Основным преимуществом EMD является высокая адаптивность, связанная с тем, что базисные функции, используемые для разложения сигнала, конструируются непосредственно из самого исследуемого сигнала, что позволяет учесть все его локальные особенности, внутреннюю структуру, присутствие различных помех. Кроме адаптивности, разложение обладает и другими важными для практических приложений свойствами [3]: локальностью (возможностью учета локальных особенностей сигнала); ортогональностью, обеспечивающей восстановление сигнала с определенной точностью; полнотой, гарантирующей конечность числа базисных функций при конечной длительности сигнала.
Эмпирические моды (ЭМ) - это монокомпонентные составляющие сигнала, которые вместо постоянной амплитуды и частоты, как в простой гармонике, имеют меняющуюся во времени амплитуду и частоту. ЭМ не имеют строгого аналитического описания, но должны удовлетворять условиям, гарантирующим определенную симметрию и узкополосность базисных функций.
Процедура EMD содержит следующие шаги [2, 3]:
1. Определение локальных экстремумов (максимумов и минимумов) входного сигнала x(t).
2. Вычисление верхней огибающей Wy(t) по найденным локальным максимумам и нижней огибающей dij(t) по локальным минимумам. Огибающие строят чаще всего с помощью кубической сплайн-интерполяции.
3. Вычисление среднего значения огибающих (локального тренда):
j )=^ММ.
4. Вычитая локальный тренд из входного сигнала, получаем первое приближение к первой ЭМ:
j (t ) = x (t)-mj (t).
5. Продолжение шагов 1-4 для выделенных приближений к ЭМ gy(t) (так называемая процедура отсеивания) и получение новых приближений и новых локальных трендов miJ-(n). Итерации продолжаются до тех пор, пока приближения не станут удовлетворять критерию останова. При выполнении критерия останова итерации прекращают, и приближение gij(t) считают эмпирической модой f (t) = gij(t).
6. Расчет остатка - разности между сигналом и ЭМ:
rM (t) — Г (t)-f■ (t).
Шаги 1-5 приводят с остатком ri+1(t). Всю процедуру прекращают, когда в реализации ri+1(t) становится мало экстремумов.
В результате декомпозиции из исходного сигнала x(t) извлекаются конечное число ЭМ fi(t) и результирующий остаток r^t), который представляет собой либо постоянную величину, либо медленно меняющийся тренд, не подлежащий дальнейшему разложению (поскольку у него отсутствуют экстремумы, использующиеся при конструировании огибающих):
V -1
x(t) = X fr (t) + V (t).
i—1
При этом каждая функция f(t) является монокомпонентной, имеет свой характерный временной масштаб осцилляций, который убывает с ростом ее номера i , а процесс вычисления сводится к устранению локального тренда, соответствующего данному масштабу.
В работах [4, 5] показана эффективность применения эмпирической модовой декомпозиции для подавления помех в ЭКС.
Усеченный алгоритм EMD
Базовый алгоритм подавления ВЧ помех на основе EMD предусматривает декомпозицию сигнала на ЭМ, нелинейную пороговую обработку (НПО) некоторых из них и последующую реконструкцию сигнала. Пороговая обработка (Thresholding, Thr) ЭМ представляет собой специальную нелинейную процедуру дискриминации отсчетов на некотором пороговом уровне в соответствии с определенными правилами:
Thr (f (t )) = ¥(p),
где f - эмпирическая мода i-го уровня разложения; t - дискретное время; у - некоторая нелинейная пороговая функция;p - пороговое значение.
Наиболее серьезным недостатком вышеописанного базового алгоритма подавления ВЧ помех в ЭКС является вычислительная сложность алгоритма EMD, затрудняющая реализацию данного способа в реальном времени. С целью уменьшения времени разложения сигнала на монокомпонентные частотные составляющие авторами предложено исключить процедуру отсеивания из алгоритма EMD. При этом первые приближения к соответствующим ЭМ будут считаться частотными компонентами входного сигнала, и именно к ним необходимо применить НПО.
Частотные компоненты, полученные в результате такого усеченного EMD, не полностью удовлетворяют свойствам ЭМ (сформулированным Хуангом [2, 3]), что может приводить к появлению в них отрицательных мгновенных частот. В данной ситуации это не имеет значения, так как задача получения спектра Гильберта и соответственно построения 3D поверхности перед нами не стоит. А вот то, что базис является законченным и сходящимся (сумма всех частотных компонент и остатка равна исходному сигналу) не вызывает сомнений. И, что самое главное, базис остается адаптивным, так как получен непосредственно из анализируемых данных эмпирическим методом.
При реализации классического EMD выполняется 6-8 итераций для качественного отсеивания ЭМ, а при разложении участка ЭКС в 2500-3000 отсчетов получается 8-12 ЭМ. Таким образом, при усечении EMD количество итераций сокращается в 48-96 раз. Даже если для устранения ВЧ помех использовать только три первые частотные компоненты, то количество операций сократится примерно в 18-24 раза. Таким образом, алгоритм фильтрации на основе усеченной EMD декомпозиции позволяет значительно уменьшить вычислительные затраты при реализации алгоритма и позволяет реализовать фильтрацию ЭКС в реальном времени. Проведенные исследования показали, что эффективность подавления помех в большей степени зависит от типа пороговой функции и уровня порога, чем от качества отсеивания ЭМ. Следовательно, применение усеченного EMD для подавления помех в ЭКС вполне оправдано.
Задачи нелинейной пороговой обработки эмпирических мод ЭКС
Стратегия НПО определяет эффективность помехоподавления и включает решение следующих задач:
- выбор количества ЭМ, к которым необходимо применить пороговую обработку;
- выбор типа пороговой функции и способа ее применения;
- определение порогового значения p.
Первая задача решается исходя из следующих соображений. Общее количество ЭМ, на которое раскладывается участок сигнала длительностью в N отсчетов, приближенно оценивается как
M = (log2 N) ± 1.
Двоичный логарифм в формуле соответствует бинарной структуре алгоритма эмпирической модовой декомпозиции. Высокочастотные помехи после декомпозиции сосредоточиваются в первых нескольких ЭМ, их количество назовем «глубиной шума», которая приблизительно равна DN = M/4 и обычно не превышает 5. Например, участок ЭКС в 2500 отсчетов (5 с при частоте дискретизации 500 отсчетов в секунду) раскладывается на 10-12 ЭМ в зависимости от сложности сигнала, вида и интенсивности помех. При этом все ВЧ помехи сосредоточиваются в первых трех ЭМ, поэтому НПО целесообразно применять именно к трем первым (наиболее высокочастотным ЭМ).
Вторая задача связана с выбором решающего правила пороговой обработки и способом его применения. На сегодняшний момент известны несколько таких правил.
Простейшим решающим правилом пороговой обработки является так называемый Hard Thresholding (HThr), сохраняющий неизменными все отсчеты ЭМ f(t), большие или равные порогу p, и обращающий отсчеты, не удовлетворяющие данному условию, в нуль:
Решающее правило Hard Thresholding имеет два недостатка, которые снижают его ценность для задачи помехоподавления. Первый из них состоит в том, что сохранение отсчетов ЭМ, превышающих некоторое заданное значение порога, подразумевает сохранение также помех, присутствующих в них. Второй недостаток связан с возникновением в результирующем сигнале паразитных гармоник за счет обращенных в нуль отсчетов [4].
Решающее правило Soft Thresholding (SThr) наряду с обращением в нуль отсчетов ЭМ f(t), содержащих на практике лишь помеховую составляющую, производит уменьшение отсчетов на величину p, что соответствует шумоподавлению также в информативных отсчетах [5]:
Функция Soft Thresholding имеет тенденцию большего искажения сигнала из-за уменьшения отсчетов, соответствующих QRS комплексам ЭКС.
Для устранения недостатков описанных пороговых правил за последние годы был предложен целый ряд других функций (гиперболическая функция Видаковича, устойчивая функция Брюса и Гао, функция Nonnegative Garrote Бреймана, функция на основе гиперболического тангенса Пурначандры и Кумаравеля и т.д.). Некоторые из них позволяют уменьшить искажения ЭКС при помехоподавлении.
Для подавления помех в ЭКС представляет интерес функция Nonnegative Garrote:
введенная Брейманом в 1995 г. [6] для устранения недостатков функций ИТЬг и 8ТЬг. Применение этой функции (при оптимальном пороге) позволяет на 1-2 % уменьшить среднеквадратичную ошибку и примерно на 5 % - максимальную абсолютную ошибку реконструкции ЭКС.
Что касается способа применения пороговой функции, то проведенные исследования показали:
- для каждой ЭМ значения порогов определяются индивидуально в соответствии с уровнем помех в ней (спектр реальных помех отличается от спектра белого шума, поэтому интенсивность помех в разных ЭМ будет различной);
- значения порогов должны корректироваться (в идеале на каждом кардиоцикле) из-за нестационарного поведения помехи.
Третья задача - определение порогового значения р для каждой ЭМ, которое связано с уровнем шума в этой моде. Известные подходы к определению значения порогов, реализованные в вейвлет-шумоподавлении (методы sqrt-log, Берга-Массара, Штайна) и основанные на
' f (t)-p,if| f (t )|> p; SThr (f (t )) = J 0, if | f (t ))< p; f (t) + p, iff (t)j<-p.
2
предположении о нормальности распределения помехи, имеют ограниченное практическое применение для подавления помех на основе БМБ и НПО в связи со следующим обстоятельством. Широко распространенный подход к определению значения порогов на основе предположения о гауссовом (нормальном) распределении шума на разных уровнях разложения не соответствует действительности. С одной стороны, помехи в ЭКС не всегда имеют гауссову плотность распределения вероятностей (ПРВ), с другой стороны, в отличие от вейвлет-преобразования линейность оператора, соответствующего процедуре БМБ, строго не доказана, поэтому даже для гауссового шума распределение шумовых составляющих на разных уровнях декомпозиции не сохраняется гауссовым.
Обычно значение порога стараются установить на уровне максимальных отсчетов шума или даже на уровне возможных выбросов при данном распределении и количестве отсчетов. Значение порога, рассчитанное по формуле Крамера
составляет около 3 а для 100 отсчетов, около 4а для 2500 отсчетов, почти 5 а для 100 000 отсчетов и является завышенным практически для всех реальных ситуаций. Несомненно, такое значение порога обеспечивает гарантированное удаление шумов, но приводит к значительному искажению сигнала. На рис. 1 приведены графики зависимости процентной среднеквадратичной ошибки - PRD (а) - и максимальной абсолютной ошибки - MaxAE (б) - реконструкции ЭКС от значения порога при использовании пороговой функции Nonnegative Garrote. Из графиков видно, что необоснованное завышение порога p приводит к увеличению ошибки реконструкции (особенно для максимальной абсолютной ошибки).
Для оценки уровня шума в ЭМ ЭКС целесообразно использовать робастные оценки масштаба [7]. Под робастностью в статистике понимают нечувствительность к различным отклонениям и неоднородностям в выборке. Робастные методы позволяют оценить масштаб шума (эквивалент среднеквадратического отклонения в классической статистике) в конкретной ЭМ, уменьшив или даже исключив влияние выбросов, соответствующих информативным участкам ЭКС.
p = aV 2ln n,
PRD: %
МО P: U f
MaxAE. е.м.р.
Выбор робастной оценки масштаба для определения порога дискриминации
Одной из наиболее эффективных и известных робастных оценок масштаба является абсолютное медианное отклонение (Median Absolute Deviation, MAD), которое описывается следующим выражением [8]:
где b - константа, зависящая от ПРВ выборки (для гауссового распределения b = 1,4826).
MAD имеет наилучшую возможную устойчивость (робастность), равную 50 %. Робаст-ность оценки определяется количеством отсчетов выборки, которые можно изменить, чтобы это не повлияло на конечный результат. Оценка MAD приспособлена для симметричных распределений и имеет низкую эффективность для гауссовых распределений.
В работе [8] предложены другие робастные оценки, в частности S- и ^-оценки. S-оценка определяется выражением
где с - константа, предназначенная для адаптации оценки к виду ПРВ выборки (для гауссова распределения с = 1,1926). Эта оценка вычисляется через расстояния между наблюдениями, поэтому применима как для симметричных, так и для несимметричных распределений.
^-оценка также применима для несимметричных ПРВ и достаточно эффективна для гауссовых распределений. Вычисляется ^-оценка следующим образом:
где нижний индекс 0,25 означает взятие 0,25-го квантиля полученной выборки разностей. Константа й также используется для согласования с видом ПРВ выборки (для гауссовых распределений й = 2,2219).
Наибольший интерес для практического применения представляет ^-оценка. Это связано с тем, что распределение отсчетов ЭМ в общем случае негауссово и несимметрично из-за наличия информативных участков ЭКС.
Алгоритм обработки эмпирических мод при подавлении высокочастотных помех в ЭКС
Исследования, проведенные авторами, позволили разработать алгоритм нелинейной обработки эмпирических мод при подавлении высокочастотных помех в ЭКС на основе усеченной EMD, который заключается в следующем:
- для каждой ЭМ, подлежащей НПО, рассчитывается робастная ^-оценка;
- выбирается значение порога рi = 2Qi , при этом удаляется большая часть выбросов шума (проведенные эксперименты показали, что в интервале ±2Qi лежит более 95 % выбросов шума в эмпирических модах реальных ЭКС);
- проводится пороговая обработка каждой ЭМ в соответствии с правилом ОТЪг(/^(0);
- оставшиеся редкие выбросы (практически всегда одиночные) устраняются, для этого целесообразно использовать трехточечный медианный фильтр.
Предложенный алгоритм нелинейной обработки ЭМ позволяет устранить высокочастотные помехи в ЭКС с минимальным искажением полезного сигнала. При этом масштаб шума оценивается без сегментации ЭМ, что позволяет использовать этот алгоритм даже для патологических ЭКС, в которых отсутствуют QRS комплексы.
Применение для пороговой обработки функции Nonnegative Garrote и отказ от необоснованного завышения порога приводят к уменьшению искажений ЭКС при подавлении высокочастотных помех.
Использование усеченной EMD в алгоритме фильтрации позволяет значительно повысить скорость обработки ЭКС, при этом появляется возможность проводить подавление помех в ЭКС в реальном времени.
MAD n = b х med, х, - med ,-x
n
j j I'
Sn = cхmed, jmedy|x, -xjЦ,
Заключение
Предложенные алгоритмы обеспечивают повышение эффективности подавления высокочастотных помех в ЭКС.
Список литературы
1. Дроздов, Д. В. Влияние фильтрации на диагностические свойства биосигналов / Д. В. Дроздов // Функциональная диагностика. - 2011. - № 3. - С. 75-78.
2. Huang, N. E. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis / N. E. Huang, Z. Shen, S. R. Long // Proc. R. Soc. Lond. A. -1998. - Vol. 454. - P. 903-995.
3. Huang, N. E. The Hilbert-Huang Transform and Its Applications / N. E. Huang, S. S. P. Shen // World Scientific Publishing Co. Pte. Ltd. - Singapore, 2005. - 526 р.
4. Кривоногов, Л. Ю. Подавление помех в электрокардиосигналах на основе разложения по эмпирическим модам / Л. Ю. Кривоногов, А. Ю. Тычков // Известия ЮФУ. Технические науки. Тематический выпуск «Медицинские информационные системы». -2010. - № 8 (109). - С. 127-132.
5. Кривоногов, Л. Ю. Перспективы применения преобразования Гильберта-Хуанга для автоматизированной обработки электрокардиосигналов / Л. Ю. Кривоногов, А. Ю. Тычков // Биотехносфера. - 2011. - № 5-6 (17-18). - С. 41-48.
6. Breiman, L. Better Subset Regression Using the Nonnegative Garrote / L. Breiman // Tech-nometrics. - 1995. - № 37 (4). - P. 373-384.
7. Хьюбер, П. Робастность в статистике / П. Хьюбер. - М. : Мир, 1984. - 303 с.
8. Rousseeuw, P. J. Alternatives to the Median Absolute Deviation / P. J. Rousseeuw, C. Croux // Journal of the American Statistical Association. - 1993. - Vol. 88, № 424. - Р. 1273-1283.
Кривоногов Леонид Юрьевич
кандидат технических наук, доцент
кафедры медицинских информационных систем
и технологий,
Пензенский государственный университет E-mail: leonidkrivonogov@yandex.ru
Папшев Дмитрий Викторович
студент,
Пензенский государственный университет E-mail: roVer_d@mail.ru
Krivonogov Leonid Yur'evich
candidate of technical sciences, associate professor, sub-department of medical information systems and technology, Penza State University
Papshev Dmitriy Viktorovich
student,
Penza State University
УДК 53.088.7, 612.172.4 Кривоногов, Л. Ю.
Повышение эффективности подавления высокочастотных помех в электрокардиосигналах /
Л. Ю. Кривоногов, Д. В. Папшев // Измерение. Мониторинг. Управление. Контроль. - 2014. - № 2 (8). -С. 17-24.
