СПОСОБ ФОРМИРОВАНИЯ ДВУМЕРНОГО ИЗОБРАЖЕНИЯ БИОСИГНАЛА И ЕГО АНАЛИЗА

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

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

Известен способ представления ЭКС, посредством которого кардиосигнал расчленяют на RR-отрезки, которые затем накладывают последовательно один на другой, синхронизируя их по максимуму RR-зубца на кардиомониторе. В результате наложения RR-отрезков друг на друга изображение ЭКС "оживает", становятся заметными малейшие флюктуации временных интервалов, амплитуд и форм зубцов ЭКС. Здесь возможны два вида наложения: со стиранием предыдущих отрезков и без стирания. В первом случае оператору предоставляется возможность оперативного обнаружения отклонений любого из параметров ЭКС и качественной оценки диапазона отклонений. Во втором случае оператор располагает интегральной картиной изменения структуры ЭКС, позволяющей произвести количественную оценку диапазона обнаруженного отклонения параметра (патент РФ №2033076 МПК⁷ A61B 5/04 «Способ представления электрокардиосигнала» [Текст] / Бакаев В.М.; Бакаев М.В.).

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

Наиболее близким к заявленному способу является способ формирования двумерного изображения биосигнала и его анализа, заключающийся в том, что в каждой точке исходного временного сигнала y(t) численными методами оценивается первая производная dy/dt и вся последующая обработка сигнала выполняется на фазовой плоскости в координатах y(t) - dy/dt. Компьютерная обработка ЭКГ предусматривает разделение фазовой траектории на отдельные сердечные циклы, селекцию траекторий с одинаковой морфологией (отбраковка ненадежных траекторий, вызванных артефактами либо экстрасистолами), усреднение траекторий в фазовом пространстве с последующей оценкой "эталонного" цикла во временной области по усредненной фазовой траектории. Этот способ позволяет одновременно оценивать как амплитудные, так и скоростные параметры любых элементов электрокардиосигнала, что дает возможность с высокой точностью оценить форму ЭКГ и обнаружить в ней такие отклонения, которые обычно скрыты от врача при традиционном анализе ЭКГ во временной области (Файнзильберг Л.С. Компьютерный анализ и интерпретация электрокардиограмм в фазовом пространстве [Текст] / Системнi дослiдження та iнформацiйнi технологii, 2004, №1. С. 32-46).

Недостатки данного способа формирования двумерного изображения биосигнала и его анализа обуславливаются следующими факторами. Отличительной особенностью способа является использование дополнительной информации, содержащейся в скоростных характеристиках исследуемого процесса. Однако скоростные характеристики процесса, определяемые на основе дифференцирования цифрового сигнала, не являются однозначными и определяются как выбранными правилами (параметрами окна цифрового фильтра) дифференцирования, так и способом фильтрации исходного сигнала. При анализе посредством этого способа имеется возможность использовать статистические параметры изображения, например, для электрокардиосигнала это параметр σ (рассеивание точек фазовых траекторий), дающий интегральное представление о вариабельности морфологии, отдельных циклов ЭКГ; угол α ориентации усредненной фазовой траектории, главным образом характеризующий соотношение амплитуд зубцов комплекса QRS; параметры симметрии отдельных фрагментов усредненной фазовой траектории относительно оси , которые характеризуют соотношение скоростей на восходящем и нисходящем участках соответствующих волн, в частности волны Т. Однако эти параметры не являются самодостаточными и используются как дополнительные параметры при оценке ЭКГ в стандартных отведениях.

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

Поставленная задача достигается тем, что в исследуемой живой системе выбирают две подсистемы, процессы в которых наблюдаются в виде синхронизированных квазипериодических сигналов x₁(t) и x₂(t), осуществляют синхронную запись этих сигналов в течение времени Т, определяемом целями исследований и природой исследуемой системы, в координатах x₁(t)-x₂(t) строят фазовый портрет исследуемой динамической системы на интервале Т, затем определяют замкнутую кривую, оконтуривающую этот портрет, определяют дескрипторы Фурье полученной замкнутой кривой, после чего осуществляют классификацию исследуемой системы посредством обучаемого классификатора, построенного в пространстве дескрипторов Фурье.

Для обучения классификатора формируют обучающую выборку, состоящую из изображений замкнутых кривых, оконтуривающих фазовые портреты динамических систем известных классов, построенных в координатах x₁(t)-x₂(t), определяют дескрипторы Фурье всех кривых, входящих в обучающую выборку, определяют число дескрипторов K_i для каждой i-й кривой обучающей выборки, определяют кривую с максимальным числом дескрипторов К_max, к дескрипторам каждой i-й кривой в высокочастотной области добавляют К_max-K_i дескрипторов с нулевым значением амплитуды, все спектральные составляющие в спектре i-го контура умножают на величину K_max/K_i, задавшись ошибкой восстановления кривых по дескрипторам Фурье и оптимизирующим критерием, определяют число дескрипторов К, используемых для классификации динамической системы, и обучают нейронную сеть для классификации кривых с К нейронами в первом слое.

Для классификации исследуемой динамической системы по дескрипторам Фурье к ее дескрипторам Фурье в высокочастотную область добавляют К_max - дескрипторов Фурье с нулевой амплитудой, умножают все дескрипторы на величину К_max/ и из К_max полученных дескрипторов выбирают К дескрипторов с младшими номерами, которые подают на К входных нейронов классификатора.

На фиг.1 показана структурная схема устройства, реализующего данный способ.

На фиг.2 представлена схема алгоритма, реализующего представленный способ.

На фиг.3 представлены детализированные схемы алгоритмов, реализующие отдельные укрупненные блоки схемы алгоритма фиг.2.

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

На фиг.5 показан фазовый портрет динамической системы, построенный на основе двух сигналов, представленных на фиг.4.

На фиг.6 показана замкнутая кривая, полученная в результате оконтуривания фазового портрета фиг.5.

На фиг.7 показаны дескрипторы Фурье оконтуривающей замкнутой кривой, представленной на фиг.6.

На фиг.8 показаны примеры восстановления границы контура по дескрипторам Фурье: М - число коэффициентов Фурье, использованных при восстановлении.

На фиг.9 приведены примеры параметрических кривых исходной и восстановленной оконтуривающей кривой фазового портрета.

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

Способ осуществляется с помощью устройства, структурная схема которого показана на фиг.1. В качестве квазипериодических сигналов x₁(t) и x₂(t) в данном устройстве используются фотоплетизмосигнал и ЭКС. Устройство состоит из датчика пульса 1 - фотоэлектрический датчик (устанавливается на подушечку большого пальца левой руки), датчика ЭКС 2 (усилитель биопотенциалов с электродами), двухканального аналого-цифрового интерфейса 3, к двум входам которого подключены выходы датчиков 1 и 2, ЭВМ 4, к системной шине которой подключен аналоговый интерфейс 3, клавиатуры 5, подключенной к порту ЭВМ 4, и монитора 6, подключенного к выходу ЭВМ 4.

Способ осуществляется согласно схеме алгоритма, представленной на фиг.2. В блоке 1 осуществляется ввод в компьютер времени мониторинга биосигналов Т и шага дискретизации Δt. В блоке 2 осуществляется синхронный ввод в компьютер двух сигналов x₁(t) и x₂(t), отражающих динамическое состояние системы. Детализированная схема алгоритма синхронного ввода представлена на фиг.3а.

В блоке 3 осуществляется построение фазового портрета динамической системы в координатах x₁(t)-x₂(t). Детализированная схема алгоритма построения фазового портрета в этих координатах представлена на фиг.3б. В блоке 13 определяют размер портрета в пикселях, а в блоке 14 - реальные динамические диапазоны сигналов. В блоке 15 определяют масштабные коэффициенты, позволяющие выразить отсчеты сигналов в координатах пикселей. В блоках 16 и 17 непосредственно получаем фазовый портрет системы путем активации пикселей в соответствующих координатах. При этом вполне возможно, что некоторые пиксели будут активироваться неоднократно, что при рассматриваемом способе построения двумерного изображения не имеет значения.

Оконтуривание фазового портрета осуществляется путем использования морфологических операций, реализованных в блоках 4 и 5. Непосредственно оконтуривание фазового портрета осуществляют посредством морфологической операции bwperim - выделение границы бинарного объекта (блок 5 фиг.2). Так как морфологический оператор bwperim оконтуривает бинарный объект без дыр, то фазовый портрет не должен содержать внутренних дыр. Для выполнения этого условия в блоке 4 осуществляют морфологическую операцию дилатации.

Замкнутую кривую фиг.6 разлагают в ряд Фурье, определяют соответствующие дескрипторы Фурье, модули которых используют как информативные признаки на входе нейронной сети прямого распространения (блоки 6 и 7 фиг.2).

На фиг.7 показаны модули дескрипторов Фурье, полученные для контура, показанного на фиг.6.

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

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

На фиг.3в представлена схема алгоритма работы блока классификации 7 в режиме обучения. На первом этапе обучения классификатора задают общее число отсчетов в оконтуривающих кривых, которое должно быть одинаковым для всех контуров контрольных и обучающих выборок. Это число определяют по результатам статистических исследований. Для фазового портрета, отражающего динамику сердечно-сосудистой системы, число отсчетов (пикселей) в контурах колеблется от 500 до 30000. Для того чтобы иметь возможность сравнивать дискретные отсчеты частот, соответствующие разным контурам, необходимо, чтобы контуры имели одно и то же количество отсчетов (априорно полагается, что частота дискретизации у них одинакова и равна единице, то есть один пиксель). Чтобы выровнять число отсчетов в выборке анализируемых контуров, необходимо довести число отсчетов в каждом контуре до максимального в выборке К_max.

С этой целью высокочастотную часть спектральной полосы заполняют нулями и тем самым доводят число отсчетов в спектре каждого контура границы сегмента до максимального, например 30000.

Таким образом, в процессе обучения классификатора осуществляют следующую цепочку преобразований: Kmax→ дополнение нулями спектральных отсчетов в области высоких частот до К_max, которая осуществляется в блоке 19.

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

Дополнение нулями спектра соответствует интерполяции в пространстве сигналов, то есть появлению дополнительных отсчетов между исходными отсчетами. Так как кривая дискретизирована с предельной частотой дискретизации, то промежуточные отсчеты попадают в те же координаты, что и исходные. В том случае, если спектр определяется на контуре, то координаты отсчета определяют и величину отсчета. Тогда дополнительные отсчеты имеют значения, равные значениям в смежных узлах интерполяции. Это ведет к нарушению требований равенства Парсеваля. Чтобы сохранить энергетический эквивалент между пространством сигналов и пространством частот, необходимо увеличить амплитуду дескрипторов пропорционально числу промежуточных отсчетов, появившихся в пространстве сигналов в результате дополнения спектра нулями. Поэтому на втором этапе обучения классификатора все спектральные составляющие в спектрах i-й кривой умножают на величину K_max/K_i, где K_i - количество отсчетов в i-й кривой (блок 20).

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

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

Предположим, что вместо всех коэффициентов Фурье α(u) используются только первые М из них. Это значит, что при u>М-1 α(u)=0. Результатом восстановления окажется следующее приближение :

для

Так как при вычислении каждой компоненты используется лишь М членов, k по-прежнему пробегает весь диапазон от 0 до К_max-1, т.е. в приближенной кривой будет то же самое число точек, но для восстановления их координат используется меньшее число дескрипторов. Из рассмотрения преобразования Фурье следует, что высокочастотные составляющие описывают мелкие детали, тогда как низкочастотные компоненты определяют общую форму границы. Поэтому, чем меньше М, тем больше деталей границы теряется.

На фиг.8 представлена граница произвольного контура, состоящая из К=130 точек, а также результаты ее восстановления с помощью уравнения (1) при разных значениях М. Следует отметить, что, начиная со значений М около 33, у восстановленной границы угловые точки начинают «выдаваться» из последовательности. При М=98 восстанавливается почти точная копия оригинала. Отсюда следует, что нескольких коэффициентов низшего порядка достаточно для описания общей формы границы, однако для точного восстановления резких деталей, например углов и прямолинейных участков, требуется значительно большее число членов высокого порядка.

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

x_k=f₁(k) (2)

y_k=f₂(k). (3)

После перехода от К дескрипторов к М дескрипторам (М<К) получают также параметрические кривые, отражающие геометрию границы фазового портрета в виде

и

Примеры исходных и восстановленных по М дескрипторам параметрических кривых (2), (3) и (4), (5) представлены на фиг.9.

Ошибка восстановления определяется из следующего выражения:

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

Для оптимизации числа дескрипторов Кэттель (Халафян, А.А. STATISTICA 6. Статистический анализ данных. [Текст] / А.А.Халафян. - М.: ООО Бином-Пресс, 2007. - 512 с.) предложил найти такое место на графиках типа фиг.9, где убывание критерия Λ слева направо максимально замедляется. Предполагается, что справа от этой точки находится только «информативная осыпь» - «осыпь» является геологическим термином, обозначающим обломки горных пород, скапливающиеся в нижней части скалистого склона. В соответствии с этим критерием можно оставить в примере фиг.10 для тестового контура 1-10 дескрипторов, а для тестового контура 2-40 дескрипторов.

Процесс определения оптимального К осуществляется в блоках 20-25.

Таким образом, перед обучением нейросетевого классификатора по обучающей выборке определяются максимальное число в обучающей выборки К_max и число нейронов во входном слое классификатора К. Затем осуществляется настройка нейронной сети (блок 27), выполняемая по известным алгоритмам, например по алгоритму обратного распространения ошибки (Осовский С. Нейронные сети для обработки информации [Текст] // С.Осовский / Пер. с польского И.Д.Рудинского. - М.: Финансы и статистика, 2004. - 344 с).

Схема алгоритма работы блока 7 в режиме классификации оконтуривающей кривой представлена на фиг.3г.

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

Конкретные примеры

Пример 1. Пациентка Ф., 57 лет, обратилась в клинику с жалобами на повышение уровня АД (150/95 мм рт.ст.), сердцебиение, головные боли, иногда приливы жара, повышенную потливость, раздражительность, чувство беспокойства, сухость кожи, лишний вес.

Из анамнеза известно, что в течение последних двух лет отмечаются эпизоды повышения АД, в связи с чем обращалась к кардиологу, но практически не лечилась. В течение последнего года присоединились эпизоды сердцебиения, раздражительность, чувство беспокойства, сухость кожи. Менопауза с 56 лет. При обследовании по месту жительства однократно уровень глюкозы крови составлял 5,8 ммоль/л. Случаев раннего развития ССЗ у ближайших родственников не выявлено.

При физикальном обследовании: состояние удовлетворительное. Рост 166 см, вес 88 кг (ИМТ 33,09 кг/м², индекс ОТ/ОБ 0,93). Кожные покровы чистые, отеков нет. Частота дыхания 18 в 1 мин. Аускультативная картина в легких и сердце - без особенностей. ЧСС 82 в мин, АД 130/80 мм рт.ст., живот при пальпации мягкий, безболезненный, печень не увеличена.

По результатам клинического анализа крови и общего анализа мочи патологических изменений не выявлено. ЭКГ: ритм синусовый с ЧСС 80 уд. в 1 мин. Отклонение электрической оси сердца влево.

По данным ЭхоКГ: существенных отклонений от нормы не выявлено.

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

Пациентке были даны подробные рекомендации по изменению образа жизни, диетические рекомендации, разъяснена важность контроля АД, уровня глюкозы крови.

После проведенного обследования пациентке был назначен моксонидин (Физиотенз, Solvay Pharma) в суточной дозе 0,4 мг, который она принимала в течение 12 недель.

При анализе полученных результатов можно отметить, что помимо достижения целевого уровня АД по данным СМАД уменьшилась выраженность вазомоторных и психоэмоциональных нарушений. Пациентка субъективно отметила улучшение общего состояния, снизился вес тела на 6 кг, уменьшились головные боли. Кроме того, снизился индекс ОТ/ОБ на 4,1%, а также ИМТ на 6,77% от исходного. Уровни общего холестерина и триглицеридов снизились на 5,92 и 11,4% соответственно. Уровень тощаковой глюкозы на фоне лечения моксонидином 0,4 мг/сут также снизился на 12,7% от исходных цифр, кроме того, улучшились показатели перорального глюкозотолерантного теста.

Пример 2. Пациент А., 67 лет. Диагноз: АГ, II стадии, третьей степени, кризовое течение. Обследование осуществляется с помощью выявления факторов риска и лабораторного исследования показателей холестерина, триглицеридов и β-липопротеидов.

Исследования фазового портрета кардиосигналов посредством нейросетевого классификатора, настроенного на диагностику артериальной гипертензии, показали высоком риске развития АГ (четвертый класс риска).

На момент обследования у пациента действительно был уже установлен диагноз АГ.

Данный способ является простым, экономичным для оценки показателя риска развития АГ на основе выбранных факторов риска. Он может использоваться при диспансеризации и профилактических осмотрах лиц от 20 до 76 лет в амбулаторных и стационарных условиях.

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