Тихонов Э.П.

Санкт-Петербургский Государственный электротехнический университет «ЛЭТИ», Россия

Адаптивные алгоритмы измерения ЭКГ с высоким

временным разрешением в условиях воздействия

помех

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

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

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

y(nτ) = x(nτ)+g(nτ)+ξ(nτ), n = 0,1,2. …, (1)

где y(nτ) – дискретные значения сигнала на выходе АЦП;

x(nτ) – реальный ЭКС;

g(nτ) – сетевая помеха.

ξ(nτ) – широкополосный шум или помеха общего вида;

τ – интервал дискретизации исходного или входного сигнала y(t).

Первостепенной задачей является измерение ЭКС x(nτ) с минимальной или заданной погрешностью в условиях воздействия помех. При этом адекватность математической модели (1) реальному сигналу достаточна для решения поставленной задачи измерения ЭКС. Недостатком известных способов фильтрации является то, что алгоритм фильтрации, описываемый в общем случае некоторым оператором H, воздействуя на сигнал y(nτ), преобразует все его составляющие, включая и ЭКС x(nτ), поэтому в результате применения оператора H получаем, что η(nτ) ) = H{y(nτ)} ¹ x(nτ). Поиск оператора фильтрации, при котором достигался бы эффект минимального искажения реального сигнала x(nτ) при максимальном подавлении суммарной помехи затрудняется из-за противоречивости сформулированного требования. В связи с этим необходимо синтезировать такой вид оператора H, для которого выполнялось бы условие

(2)

где K{ … } – критерий, характеризующий меру отклонения x(nτ) от η(nτ).

Синтез оператора H, удовлетворяющего (2), является достаточно сложной задачей. Решению этой задачи посвящено большое число работ (см., например, [2] и приведённую в нём библиографию). Если подходить к решению поставленной задачи с точки зрения измерения x(nτ), то условие (1) может быть представлено в виде

Р{K{x(nτ),η(nτ)} ≤ δ₀} = b, (3)

где Р{K{x(nτ),η(nτ)} ≤ δ₀} – вероятность того, что характеристика отклонения сигналов x(nτ) и η(nτ) в соответствии с функцией меры K{x(nτ),η(nτ)} не превысит априорно заданную величину δ₀, равна b.

Оценка погрешности осуществляется в виде систематической и случайной составляющих, которые оцениваются на основе оператора определения математического ожидания (МО). Средний квадрат объединяет обе составляющие погрешности. Требование (3) к измерению ЭКС оказывается достаточно жёстким, так как полная информация о виде сигнале x(nτ) отсутствует. Поэтому обычно говорят не об измерении ЭКС, а только об её оценке либо о регистрации ЭКГ. Основная задача статьи заключается в исследовании возможности осуществления (3) на базе современных информационных технологий, доступных для реализации на основе последних достижений в области микроэлектроники.

Основная часть. Обычно из всевозможных операторов в качестве H, потенциально удовлетворяющих условию (2), на практике выбирается класс операторов, обладающих свойством линейности, т. е.

η(nτ)=H{x(nτ)} +H{g(nτ)}+H{ξ(nτ)}, (4)

Равенство (4) позволяет видоизменить (2) и приблизится к условию (3) за счёт выполнения следующих требований

(5)

где ρ[ … ] – функция меры, характеризующая отклонение результата измерения ЭКС от его действительного значения x(nτ);

ε₀ – заданная или допустимая величина отклонения ЭКС от результата её измерения в соответствии с установленной функцией меры;

Θ{ … } – некоторый априори задаваемый критерий, в смысле которого осуществляется поиск минимума при воздействии оператором Н на соответствующую составляющую сигнала η(nτ).

Рассматривать оценку как результат измерения можно лишь в том случае, когда величина ε₀ не превышает допустимую для измерений величину. Очевидно, что найти или синтезировать оператор (3) со свойством (5) в отсутствии чёткого контроля ненаблюдаемого сигнала x(nτ) по наблюдаемому сигналу η(nτ) затруднительно, тем более, в условиях, когда для решения задачи (2) необходимо выполнить противоречивые требования (5). Истинное значение сигнала x(nτ) неизвестно и к тому же вариативно и проверка первого условия в (5), даже эмпирически, невозможно. В этом случае представляет интерес, как теоретически, так и практически, метод решения данной задачи, постулирующий синтез композиционного оператора H{y(nτ)} в виде

η(nτ) = Н₁{ Н₂{x(nτ)}}+Н₁{ Н₂{g(nτ)}}+Н₁{ Н₂{ξ(nτ)}}. (6)

В равенстве (6) теоретически возможно так выбрать операторы Н₁ и Н₂, чтобы они удовлетворяли требованию (2) и, как следствие, (5). При этом эти операторы должны быть достаточно эффективными для решения поставленной задачи и доступными для практической реализации. Наиболее подходящим для раздельного выполнения требований (5) является свойство линейности (4), представленного в (6), которому удовлетворяют, например, операторы интегрирования Н₁ и дифференцирования Н₂ или операторы, в силу дискретного представления сигнала, близкие к ним. Однако остаётся вопрос о возможности выполнения на основе указанных видов операторов Н₁ и Н₂ сформулированных требований (2) и (5), так как для его решения наличие свойства линейности необходимо, но, явно, недостаточно. Для того чтобы продолжить поиск решения поставленной задачи, обратим внимание на то, что дискретность представления сигналов создаёт предпосылки для перехода к синтезу дискретных операторов Н_д1 и Н_д2 , аппроксимирующих соответствующие непрерывные операторы, в виде суммы и конечных разностей. Для перехода к дискретному варианту синтеза искомых операторов, рассмотрим общий подход их синтеза на основе дифференциального уравнения первого порядка, устанавливающего равенство между скоростью изменения выходного сигнала и заданным преобразованием входного сигнала, в виде

(7)

где a – коэффициент пропорциональности;

F^*[ … ] – априорно задаваемое преобразование ЭКС, включая помехи.

Дифференциальное уравнение (7) вытекает из баланса, который наступает в текущий момент времени t между скоростью изменения выходного сигнала электрокардиографа и соответствующим преобразованием входного сигнала y(nτ). Коэффициент пропорциональности a может зависеть от времени t, а его вид, значение и размерность либо задаются априорно, либо подстраиваются на основе дополнительных алгоритмов. В рассматриваемом случае будем считать его постоянным, причём, его величина и размерность устанавливается априорно и согласованно с преобразованием F^*[ … ]. Не уточняя конкретный вид преобразования F^*[ … ], представим (7) в эквивалентной форме

dh(t) = aF^*[x(t),g(t),ξ(t)]dt (8)

Переход к дискретному варианту представления равенства (8) осуществляется интегрированием правой и левой частей уравнения (8) в пределах от nτ до (n+1)τ. В результате интегрирования получим условие аппроксимации баланса (7) на заданном интервале дискретизации τ сигнала y(nτ)

Зафиксировав для преобразования F^*[ … ] момент времени t = nτ, приходим к аппроксимации (8) в виде итерационного алгоритма, сформированного при фиксированном временном интервале дискретизации τ в следующем виде

h[(n +1)τ] = h(nτ) + aτF[x(nτ) + g(nτ) + ξ(nτ)], n = 0, 1, 2, 3, …

Выполняя итерации при h(0) = 0, получаем равенство для определения выходного сигнала в момент времени (n +1)τ в виде

(9)

Сравнивая (9) с (4) и (6) нетрудно убедиться в том, что правая часть (9) будет удовлетворять свойствам композиционного оператора Н, представленного в (6), если преобразование F[ … ] относительно суммы переменных будет обладать свойством линейности, так как сумма в (9) при надлежащем выборе интервала дискретизации t аппроксимирует интегральный оператор Н₁. Следовательно, синтез оператора Н₂ сводится к выбору конкретного вида преобразования F[ … ]. Свойство линейности будет выполняться, если преобразование F[ … ], по крайней мере, будет аппроксимировать оператор дифференцирования через конечные разности. В этом случае с учётом линейности композиционного оператора каждую составляющую исходного сигнала (3) можно рассматривать в отдельности на предмет выполнения для неё соответствующего условия (4). Аппроксимация оператора дифференцирования через конечные разности в условиях воздействия помех, как указано в работе [3], оптимизируется на основе минимума среднего квадрата для функции, заданной в дискретных отсчётах. С учётом вида указанной оптимизации введём следующую модификацию аппроксимации оператора дифференцирования по конечным разностям, предложенного в [3] в виде преобразования

(10)

где Т – временной параметр, характеризующий величину временного отклонения в операторе F[ … ] отчётов сигнала y(nτ) при i =1, 2, … m от точки дифференцирования kt;

m – число соседних точек (отсчётов сигнала) с обеих сторон относительно точки дифференцирования kt.

Представленное в (10) преобразование F[ … ] зависит от 3-х параметров τ, m и T. В общем случае параметр Т = рτ, так как он. устанавливается кратно интервалу дискретизации τ. При m = 1, 2,3, … и Т = τ, т. е. для р = 1, преобразование F[ … ] совпадает с аппроксимацией оператора дифференцирования, предложенного в [3]. Для дальнейшего решения задачи, необходимо провести анализ всех трёх составляющих сигнала y(nτ) с целью их сравнения между собой и установления для каждой из них соответствующей уточняющей модели реальной составляющей. Остановимся сначала на анализе ЭКС x(nτ). Этот сигнал имеет значительную специфику и его можно охарактеризовать как квазидетерминированный сигнал, так как наряду с априорно непредсказуемыми вариациями в достаточно широких пределах его параметров и характеристик он состоит из некоторого набора (множества) фрагментов. Факт наличия этих фрагментов повторяющихся с определённой периодичностью или, вернее квазипериодичностью, носит детерминированный характер. Обычно ЭКС разбивается на фрагменты в пределах интервала его квазипериодичности в соответствии с установленными параметрами и характеристиками, которые являются отличительными признаками для выделенного фрагмента. Из множества признаков для k-го фрагмента можно выделить существенный или главный признак. Например, для R зубца главным признаком является постоянство (в установленном смысле) модуля первой производной ЭКС на определённом временном интервале. В соответствии с квазипериодичностью ЭКС можно утверждать и об определённой вариации локальных динамических характеристик фрагментов ЭКС, например, скорости изменения или производной ЭКС в области qRs комплекса. В отличие от сигнала x(nτ) у остальных двух составляющих сигнала y(nτ) динамические свойства изменяются случайно, но эти изменения имеют стационарный характер, включая и детерминированность частоты следования сетевой помехи g(nτ). В принципе, на отличии характеристик и параметров составляющих сигнала y(nτ) построены все известные методы выделения ЭКС x(nτ) из помех. При этом априорно выбираются такие параметры фильтрующих операторов, лежащих в основе указанных методов, которые обеспечивали бы необходимый компромисс между уровнем подавления помех и искажением полезного сигнала на интегральном уровне без учёта наличия у него локальных и постоянно эволюционирующих фрагментов. Иначе говоря, при фильтрации в известных методах не учитывается локальная динамика сигнала x(nτ), обусловленная наличием у него различных по форме фрагментов, в то время как для этого сигнала для различных фрагментов условия подавления помех значительно меняются. Например, в области изменения комплекса мощность ЭКС значительно превышает мощность суммарной помехи. В связи с этим повышение эффективности измерения ЭКС x(nτ) в условиях воздействия помех целесообразно осуществить за счёт адаптации соответствующих параметров оператора H к фрагментам ЭКС. К одному из этих параметров, прежде всего, целесообразно отнести интервал дискретизации τ. Если адаптировать интервал дискретизации к фрагменту ЭКС в соответствии с характеристикой его формы, то можно с одной стороны уменьшить погрешность измерения полезного сигнала, а с другой стороны существенно повысить эффективность подавления помех. Для решения задачи адаптации интервала дискретизации к различным фрагментам ЭКС, как уже отмечалось, целесообразно воспользоваться вероятностным методом адаптивной дискретизации случайных процессов достаточно полно исследованным, например, в [4, 5]. Предварительно заметим, что в основу любого адаптивного алгоритма должен быть положен принцип установления баланса (равновесия) между теми или иными существенными для решаемой задачи параметрами или характеристиками объекта. В рассматриваемом случае устанавливается равенство между скоростями изменения интервала дискретизации в текущий момент времени и динамикой погрешности восстановления ЭКС, зависящей от его скорости изменения, на этом же интервале дискретизации в соответствии с дифференциальным уравнением

(11)

где t(t) – интервал временной дискретизации для текущего момента времени t;

ρ[ y(t),j(t,t(t)] – функция меры, характеризующая отклонение результата измерения сигнала y(t) от его восстановленного значения j(t,τ(t)) на интервале дискретизации τ(t) в момент времени t;

y{ … } – производная для одномерного случая или градиент для многомерного варианта от функции меры, характеризующей отклонение в момент времени t функции меры ρ[ y(t),j(t,t(t)] от её априорно задаваемого значения d₀;

a – коэффициент пропорциональности.

Для того чтобы практически использовать уравнение (11) для решения задачи адаптации на основе реальной цифровой системы по выходному сигналу АЦП необходимо перейти от дифференциального уравнения к его аппроксимации в конечных разностях в итерационной форме

t[(n+1)Dt] =t(nDt)+ b{y{d₀ – r[y (nDt ),j( nDt, t(nDt))]}}, (12)

где b = aDt, причём для сходимости алгоритма требуется, чтобы aDt < 1;

t[(n+1)Dt] и t(nDt) – интервалы временной дискретизации для n + 1-го и n-го тактов итерации

Обычно желательно контролировать адаптивный алгоритм формирования интервала дискретизации (12) не по абсолютной величине d₀, а по её относительному аналогу, задаваемому в виде n₀r[y(nDt)] при n₀ < 1.

На основании изложенного функционирование адаптивной система для измерения ЭКС в общем случае описывается математически нелинейной системой вида

(13)

аналитическое решение которой достаточно проблематично даже при задании в явном виде всех преобразований, учитывая, что сигнал y(t), с учётом всех его составляющих в пределах интервала периодичности x(t), не стационарен.

Прежде всего, отметим, что в силу дискретности представления интервала дискретизации в компьютерных ЭКГ в алгоритме (12) постоянный параметр b равен единице, а преобразование y{ … } соответствует знаковой функции. В этом случае динамика изменения интервала дискретизации может быть описана графом, приведённым на рис. 1.

Рис. 1 Граф, характеризующий последовательность изменения интервала дискретизации за один временной такт итерации

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

Учитывая, что математический анализ свойств алгоритма достаточно громоздкий и в силу отсутствия адекватной математической модели описания ЭКС, дальнейшее исследование необходимо провести на реальном сигнале с применением методов имитационного моделирования. Такое исследование для предложенных алгоритмов было выполнено для реальных ЭКС сигнала в полосе частот 0.05 - 250 Гц с различными интервалами дискретизации Dt, полученных без предварительной фильтрации сетевой помехи. На рис. 2 приведены результаты фильтрации ЭКС сигналов в соответствии с (13) и интервалами дискретизации Dt, равными 20 мкс и 40 мкс. На рис. 3 показаны некоторые интересные результаты моделирования для различных пациентов, которые дают интегральную информацию о ЭКС, связывающую эволюцию RR интервалов и фрагменты ЭКС.

Рис. 2 Результат фильтрации ЭКС с указанием адаптивного изменения интервала дискретизации в соответствии с алгоритмом (13) и исходным после АЦП интервалом дискретизации Dt, равным 40 мкс для разных значений параметров однотипного алгоритма

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

Литература

1. Тихонов Э. П. Алгоритм поиска интервала квантования случайного процесса / В кн.: Автоматический контроль и методы электрических измерений // Тр. VIII Всесоюзн. конф., т.1, Новосибирск,: Наука, СО, 1971, С.46 – 51.

2. Рангайян Р. М. Анализ биомедицинских сигналов. Практический подход. / Пер. с англ. под ред. А. П. Немирко. – М. ФИЗМАТЛИТ. 2007. – 440 с., ил.

3. К. Ланцош Практические методы прикладного анализа. Спавочное руководство / Пер. с англ. М. З. Кайнера под ред. А. М. Лопшица. – М. ФИЗМАТЛИТ. 1961. – 524 с., ил.

4. Тихонов Э. П. Вероятностные адаптивные алгоритмы дискретного представления аналоговых сигналов. Ч. 1: Исследование свойств // Информационно-управляющие системы. 2011. №2. С. 8 – 15.

5. Тихонов Э. П. Вероятностные адаптивные алгоритмы дискретного представления аналоговых сигналов. Ч. 2: сравнительный анализ и численные данные // Информационно-управляющие системы. 2011. №3. С. 9 – 14.