Methodological Recommendations for the Creation of Sensor Measurement Systems for Respiratory Rate Monitoring Based on Photoplethysmographic Signal Processing

Full Text

1. ВВЕДЕНИЕ

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

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

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

Наряду с подсчетом вручную и использованием грудных электродов для измерения частоты дыхания используются пульсометры и бесконтактные методы. Известно об использовании в Европе пульсометра Masimo MightySat Rx (фирмы Masimo SET), в котором для вычисления ЧДД используется сигнал от венозной крови в отличие от традиционной пульсоксиметрии, которая предполагает, что на участке, где проводится измерение, пульсирует только артериальная кровь. Принцип работы прибора основан на обработке плетизмограммы, отражающей изменения пульсовой волны в реальном времени. Отмечается, что этот прибор обеспечивает высокую точность измерений, но конкретные показатели точности не указаны.

Из наиболее точных приборов следует отметить респираторный достаточно дорогой монитор частоты пульса и частоты дыхания Capnostream-35 и портативный монитор пациента Kernel KN601M, назначение которого состоит в определении сатурации артериальной крови, частоты пульса на периферических сосудах и получения фотоплетизмографических сигналов с точностью определения частоты дыхания ±2 вдоха/мин (Tiara Medical, 2013).

Из бесконтактных методов известно использование устройств (Гаранин и др., 2023) на основе видеомониторинга, тепловизионной камеры (тепловизора), допплеровского и ультрашироковолнового радаров, тензодатчиков. Наряду с преимуществами этих устройств – точностью измерений, возможностью длительного мониторинга, измерения частоты дыхания через одежду – известны и недостатки, связанные с зависимостью от температуры окружающей среды, чувствительностью к окружающим предметам в кадре, освещенности, особенностями пациента и сложностями обработки данных.

В ряде исследований подтверждается, что частота дыхания может быть достоверно определена с помощью фотоплетизмографического сигнала (ФПГ), отражающего наполнение капилляров кровью и состояние микроциркуляторного русла. В работе (Lázaro et al., 2011) представлен метод, основанный на вариабельности пульса и использовании скорости и дисперсии пульсовой волны. Частота дыхания оценивалась по месту нахождения наибольшего пика в скользящем среднем спектра мощности. Расчетная ошибка оценки частоты дыхания составляет 1,3±7,8%, а погрешность оценки частоты дыхания достигает 6,7%.

Устойчивая мера отличия субъектов с различной частотой дыхания предложена в (Garde et al., 2014), где алгоритм мониторинга основан на обобщенной спектральной корренотропии, содержащей информацию о статистической и временной структуре пульсовых волн. В работе (Shelley et al., 2006) изучается частотный спектр сигнала пульсового оксиметра, его изменение во времени. На основе изменения частоты, интенсивности и амплитуды, извлекаемых из ФПГ сигнала, предложен совместный частотно-временной анализ формы сигнала пульсового оксиметра, который используется для определения частоты дыхания пациентов, находящихся на искусственной вентиляции легких.

В статье (Johansson et al., 2003) показано, что применение нейронной сети (НС) позволяет повысить точность измерения частоты дыхания. Был использован простейший вариант НС, а сигналы ФПГ в эксперименте регистрировались со лба у 15 здоровых людей. Из этих сигналов были извлечены форма систолической волны, форма диастолической волны и амплитуда пульса. Общий уровень ошибок измерения ЧДД при этом составлял 9.5–9.6%.

В исследовании (DeKorte et al., 2018) рассмотрен частотно-временной подход для определения мгновенной частоты дыхания (IRR). При этом на первом этапе реализовано частотно-временное преобразование, на основе вейвлет-анализа, для извлечения сигналов интенсивности, амплитуды и частоты (синхросжатое преобразование (SST)), на втором – преобразование SST применялось для каждого извлеченного сигнала вариации, вызванной дыханием, для оценки соответствующей IRR, и на третьим этапе проводился расчет окончательного значения IRR. Этот метод обеспечил наиболее точные оценки амплитуды и явился одним из лучших методов извлечения информации о частоте дыхания из сигналов ФПГ. Однако его сложность реализации состоит в четырехкратном применении SST-преобразования, что увеличивает время получения оценки ЧДД.

В статье (Chang et al., 2018) представлен эффективный подход к извлечению частоты дыхания из ФПГ на запястье, основанный на спектральном анализе Голо-Гильберта (Holo-Hilbert spectral analysis – HHSA)¹.

В 2016 году Хуанг (Huang et al., 2016) предложил спектральный анализ Голо-Гильберта для выявления совместной частотной и амплитудной модуляции. По сравнению с традиционным преобразованием Гильберта–Хуанга, который использует однослойный вариант эмпирической модовой декомпозиции (EMD) для учета изменений внутри волновой частотной модуляции, HHSA использует вложенный подход EMD для организации межмасштабного умножения сигналов амплитудной модуляции. HHSA позволяет расширить возможности существующего спектрального анализа и создает полное спектральное представление о нелинейных и нестационарных данных со всеми возможными режимами амплитудной и частотной модуляциями, как аддитивными, так и мультипликативными. Этот метод обеспечивает получение реализаций высокоразмерного спектра.

В статье (Chang et al., 2018) отмечается высокая точность метода, но вместе с тем указано на проблему использования HHSA, связанную с тяжелыми вычислениями при выполнении вложенного EMD. Его реализация для компактных и быстродействующих сенсорных устройств будет представлять значительные трудности и увеличивает стоимость устройств.

По этой причине в статье предложено для оценки частоты дыхания использовать преобразование Гильберта–Хуанга (HHT), которое хорошо описывает нелинейные и нестационарные характеристики сигнала, а реализуемо проще, чем HHSA.

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

• в дополнение к модовой декомпозиции преобразования Гильберта–Хуанга разработан и добавлен алгоритм определения моды, которая в наибольшей степени связана с искомой частотой дыхания;
• разработана методика выбора приоритетных регионов для оценки частоты дыхания;
• предложен критерий определения требуемого натяжения браслета при измерениях.

Оставшаяся часть статьи построена следующим образом: в разделе II представлена методика оценки частоты дыхания на основе HHT с выбором приоритетной моды, в разделе III представлены результаты исследований, в разделе IV дается обсуждение результатов, в разделе V – выводы. В приложении рассмотрены методика выбора приоритетных регионов для оценки частоты дыхания и критерий определения требуемого натяжения браслета при измерениях.

2. МЕТОДИКА ОЦЕНКИ ЧАСТОТЫ ДЫХАНИЯ НА ОСНОВЕ ПРЕОБРАЗОВАНИЯ ГИЛЬБЕРТА–ХУАНГА С ОПРЕДЕЛЕНИЕМ КОНТРОЛЬНОЙ МОДЫ

Методика предполагает:

• предварительную цифровую обработку фотоплетизмографических сигналов, – фильтрацию и аппроксимацию;
• модовую декомпозицию сигналов на основе преобразования Гильберта–Хуанга;
• определение из набора вычисленных мод наиболее связанной с частотой дыхания.
• вычисление функции коррентропии для указанной приоритетной (контрольной) моды;
• вычисление спектральной плотности функции коррентропии для оценки частоты дыхания.

Необходимость в выполнении фильтрации, как в случае других биологических сигналов (Кубланов и др., 2020), связана с требованием уменьшения воздействия шумов (артефактов и помех) на вычисляемые параметры сигналов. Наличие помех при обработке фотоплетизмографических сигналов является следствием влияния движений рукой, на которой закреплен датчик (в том числе плохого натяжения браслета с датчиком на запястье), мышечных сокращений, изменений температуры в помещении при измерениях, воздействия прямых солнечных лучей и яркого освещения.

На точность измерений могут оказывать влияние принятые перед измерениями физические процедуры, прием лекарственных средств, методические ошибки и погрешности измерительной аппаратуры. При устранении артефактов при цифровой обработке сигналов применяются известные методы фильтрации (Рангайян, 2007; Айфичер и др., 2008). Интерполяция проводится путем повышения частоты дискретизации сигнала, исходя из предположения об ограничении исходного спектра, а также для необходимости сопряжения и совместного анализа различных источников сигналов. Так как указанные статистические процедуры широко известны, они применительно к данной методике не рассматриваются.

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

Использование в методике преобразования Гильберта–Хуанга (HHT) объясняется тем, что пульсовые сигналы содержат нелинейные или несинусоидальные характеристики, которые сложно выделить и количественно оценить. Разложение на эмпирические моды предлагает потенциальное решение, разделяющих сигнал на набор эмпирических мод (IMF), которые допускают физически интерпретируемое преобразование Гильберта (Huang et al., 1998) и последующий анализ.

Преобразования Гильберта–Хуанга (Huang et al., 1998; Huang et al., 2003) является эффективным методом обработки нестационарных сигналов. Оно позволяет выявить скрытые в шумах амплитудные и частотные модуляции, устранить перекрывающиеся колебания сигнала для увеличения значимости мгновенных частот и неравномерность амплитуды, когда соседние по спектру колебания сильно отличаются друг от друга. В соответствии с этим преобразованием обработка данных проводится итеративным алгоритмом в два этапа.

На первом этапе исходные данные с помощью эмпирической модовой декомпозиции разлагаются в ряд отдельных компонентов (процесс “просеивания”), называемых эмпирическими модами (IMD) (Huang et al., 1998). Эмпирическая модовая декомпозиция (EMD) проводит разложение сигнала на несколько модовых функций (IMF), которые раскрывают нестационарные и стационарные компоненты, а также фиксируют респираторные особенности.

Задача процесса “просеивания”, во временной области, состоит в изолировании самых быстрых колебаний во временном ряду путем итеративного отсеивания более медленных. Медленные колебания удаляются путем вычитания среднего значения огибающей верхней и нижней амплитуды сигнала до тех пор, пока среднее не станет достаточно близким к нулю. Этот изолированный компонент сигнала известен как эмпирическая мода; он вычитается из исходного сигнала, и процесс просеивания повторяется для определения, следующей модовой функции, которая будет содержать более медленную составляющую. Этот процесс повторяется до тех пор, пока в сигнале останется только тренд.

Необходимо отметить преимущество преобразования, состоящее в том, что разложение данных происходит по базису, который получается из самих данных, а не выбирается из заранее известного набора. После работы алгоритма EMD сигнал x(t) может быть определен как

$x\left(t\right)=\sum _{i=1}^{n_c-1}{c}_i\left(t\right)+r_{n_c}\left(t\right)$,

где n_c – общее количество IMF, c_i(t) – эмпирические моды, $r_{n_c}$(t) – остаток разложения.

Критерий остановки процесса “просеивания” (Ш.Ч. Кан и др., 2010) может состоять в ограничении нормированной величины L, которая вычисляется из двух последовательных результатов “просеивания”, при этом остаток 𝑟_𝑛 становится монотонной функцией

$L=\sum _{t=1}^T\left[\frac{\left|h_{1\left(k-1\right)}\left(t\right)-h_{1\left(k\right)}\left(t\right)\right|}{h_{1\left(k-1\right)}^2\left(t\right)}\right]\in \left[\mathrm{0.2,} 0.3\right]$,

где h_1(k)(t) – разность при “просеивании” на k-м шаге.

Далее методика предполагает применение к эмпирическим модам преобразования Гильберта для сохранения временных локальных особенностей анализируемых временных рядов. Это позволяет оценить нестационарность сигнала за счет определения мгновенной частоты (Huang et al., 2009) и мгновенной амплитуды:

$H\left[c_i\left(t\right)\right]=\frac{1}{\pi }\underset{-\infty }{\overset{\infty }{}}\frac{c_i\left(t\right)}{t-\tau }d\tau$.

Для аналитического сигнала $z_i\left(t\right)=c_i\left(t\right)+jH\left[c_i\left(t\right)\right]=a_i\left(t\right)$, для заданного момента времени, мгновенная частота может быть определена в виде

$f_i\left(t\right)=\frac{1}{2\pi }\frac{d{\theta }_i\left(t\right)}{dt}$.

Следует заметить, что извлечение мгновенных частот зависит не от свертки (как в модели Фурье), а от производных по времени. Мгновенная частота представляет собой однозначную функция времени, – каждому заданному моменту времени соответствует только одно значение мгновенной частоты.

Ключевым отличием от ранее известных вариантов применения преобразования Гильберта–Хуанга и особенностью настоящей методики является использование разработанного алгоритма нахождения моды c_k(t), которая наиболее связана с частотой дыхания. Назовем эту моду контрольной или приоритетной. Для ее определения должны быть вычислены мгновенные частоты {f_k(t)}_D в диапазоне дыхания $D\in \left(\mathrm{0.2,}\dots \mathrm{,0.3}\right)$ Гц, который свойственен взрослому человеку в покое (в том числе и во сне), что соответствует (12–18) вд./мин (Гуцол и др., 2014). Контрольная мода является наиболее связанной с частотой дыхания в том смысле, что она обеспечивает бóльшую площадь под кривой плотности вероятности для соответствующих мгновенных частот в D-диапазоне.

Оценим распределение плотности вероятности мгновенных частот f_k(t), соответствующих контрольной моде. Для построения гистограммы найдем границы частных непересекающихся интервалов $\Delta =f_{k_j}-f_{k_{j-1}}$ мгновенных частот $f_{k_1},f_{k_2},\dots ,f_{k_n}$ удовлетворяющие следующему уравнению: $P\left(f\le f_j\right)=\frac{j}{k}$, где $k=\left\{\frac{N}{5}, при N<100\right.$, где N – объем выборки мгновенных частот.

Далее вычислим частоту попадания множества мгновенных частот в каждый из интервалов. Тогда гистограмма определится в соответствии с известной классической процедурой:

где $s_i=\frac{n_i}{\Delta N}$ – количество значений, попавших в i-й интервал функция мгновенных частот, F_k – соответствующая контрольной моде c_k(t).

Тогда площадь под кривой гистограммы в диапазоне D частот дыхания имеет вид

$S_{f_k}=F_k\left(a\le f_k\le b\right)=\sum _{f_k=a}^{f_k=b}{F}_k\left(f_k\right)$,

где a = 0.2, b = 0.33.

Для определения контрольной моды необходимо вычислить функцию плотности вероятности для каждой моды из всех найденных мод (обычно их 7–10). Контрольная мода определяется путем сравнения площадей $S_{f_j}$ под кривой плотности вероятности, для всех мод в указанном диапазоне частот. При этом контрольная мода определяется после вычисления площади $S_{f_k}=\left\{\text{max}\left(S_{f_j}\right)\right\}$, где $j=\mathrm{1,}\dots ,N$, N – общее количество рассматриваемых мод. Результат выполнения этого процесса наглядно представлен в разделе III.

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

Коррентропия вычисляется в виде формулы с учетом разложения в ряд Тейлора гауссова ядра k(x(n), x(n – m)) и содержит все моменты четного порядка случайной величины (x(n), x(n – m)) (Weifeng et al., 2007; Santamaria et al., 2006):

$V\left(m\right)=\frac{1}{\sqrt{2\pi \sigma }}\sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k}{2{\sigma }^{2k}k!}E\left[{\left(x\left(n\right)-x\left(n-m\right)\right)}^{2k}\right]$,

где

$k\left(x\left(n\right)-x\left(n-m\right)\right)=\frac{1}{\sqrt{2\pi \sigma }}\text{exp}\left[-\frac{{\left(x\left(n\right)-x\left(n-m\right)\right)}^2}{2{\sigma }^{2k}}\right]$

– функция гауссова ядра, σ − размер ядра, определяемый правилом плотности Сильвермана (Silverman. 1986; Herawati et al., 2017). Различные функции ядра дают различные спектральные разрешения, но так как ядра нелинейные, то они обеспечивают статистические оценки более высокого порядка для информации о случайном процессе. Если объем исходных данных мал, то размер ядра должен выбираться как компромисс между эффективностью оценки (малая дисперсия) и наличием выбросов. Необходимо задавать размер ядра в соответствии с сигналом ошибки каждой итерации (Weifeng et al., 2007; Silverman, 1986), тогда ${\sigma }_{SM}=1.06\min \left\{{\sigma }_E,\frac{R}{1.34}\right\}{N}^{-\frac{1}{5}}$, где σ_E – мощность ошибки, а R – квартильный размах ошибки. Величина σ_SM изменятся в пределах 1.8…3, что находится в окрестности ее наилучших значений (Weifeng et al., 2007).

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

Для стационарного случайного процесса допустимо записать оценку коррентропии в виде (Weifeng et al., 2007)

$V\left(m\right)=\frac{1}{N-m+1}\sum _{n=m}^N\left[k\left(x\left(n\right)-x\left(n-m\right)\right)\right]$.

Как указано ранее, частота дыхания определяется в диапазоне частот (0.2–0.3) Гц по коррентропии контрольной моды. Для этого вычисляется спектральная плотность функции коррентропии на основе представления дискретным преобразованием Фурье с корреляционным окном W(m):

$P_{Sp}\left(f\right)=\sum _{m=-L}^{L}W\left(m\right)V\left(m\right)\exp \left(-j2\pi fmT\right)$,

где частота $-\frac{1}{2T}\le f\le \frac{1}{2T}$, L ≈ N/10.

3. РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ

При проведении исследований использовались следующие наборы фотоплетизмограмм:

• выборка из 15 записей сигналов ФПГ с частотой дискретизации 50 Гц для здоровых людей, средний возраст которых 32 года. Измерения амплитудных значений пульсовых волн выполнено сенсором Biofy SFH7072. Для получения контрольных значений частоты давления использован прибор Biopac MP100;
• ограниченная выборка из баз данных (Real-World PPG dataset, 2019; PPG-BP Database, 2022; Elgendi et al., 2022). В первом наборе из них были сигналы ФПГ от 35 здоровых людей с частотой дискретизации 50 Гц. Данные из базы PPG-BP записей данных содержали сигналы для 19 испытуемых и были записаны с помощью прибора Biopac RES.

Обработка сигналов ФПГ проводилась в соответствии с методикой, представленной в разде ле II. Предварительный этап обработки сигналов для уменьшения влияния артефактов и помех, основан на известных методах цифровой обработки сигналов. Для улучшения разрешения вокруг пика частоты дыхания, пики, соответствующие частоте сердечных сокращений (ЧСС), отфильтровываются с использованием фильтра нижних частот в диапазоне (0.5–3) Гц. Сигналы ФПГ сегментировались с помощью временных окон (Марпл.-мл., 1980), длительностью (7–10) с с перекрытием 50%. Оценка погрешности вычисления частоты проводится путем сравнения с контрольным значением частот дыхания.

В ходе проведения исследований проводился анализ сигналов, снятых с различных регионов, например, со лба испытуемого или с кисти руки. Определение наиболее приемлемого региона, который обеспечивает сигнал с бóльшей амплитудой, было основано на использовании методики выбора приоритетных регионов для оценки ЧДД (см. приложение А).

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

В качестве примера оценки ЧДД рассмотрим результаты обработки ФПГ сигналов согласно методики раздела II. Аналогичный алгоритм обработки был применён и к другим сигналам ФПГ. На рис. 1, в качестве примера, показаны графики сигналов двух пульсовых волн, после предварительной обработки, снятых с кисти руки с помощью сенсора Biofy SFH7072. Рассмотрим результаты обработки этих сигналов.

Рис. 1. Графики двух пульсовых волн на 30 с интервале (пунктиром показаны тренды пульсовых волн)

Результаты модовой декомпозиции разложения Гильберта–Хуанга для 1-й волны имеют вид, изображенный на рис. 2.

Рис. 2. Моды (IMD) разложения Гильберта–Хуанга (вверху графика показана кривая пульсовой волны)

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

На рис. 3 показаны графики плотностей вероятности для двух пульсовых волн, в виде огибающих гистограмм мгновенных частот. Вычислены площади под кривыми плотностей вероятности, соответствующие всем модам в диапазоне D частот (0.2…0.3) Гц (на графиках показаны кривые только для трех мод, так как для других мод кривые проходят в стороне от указанного диапазона частот дыхания).

Рис. 3. Графики плотностей вероятности для двух пульсовых волн, показаны 2-я, 3-я и 4-я моды с отмеченным диапазоном дыхания S_k, 3-й моде соответствует кривая, отмеченная синим цветом. Графики (a) соответствуют 1-й пульсовой волне, графики (b) – второй волне

Из графиков видно, что для 1-й волны площадь S_k под кривой, которая соответствует 3-й моде (PSD[3]) в диапазоне частот дыхания, больше, чем соответствующая площадь для 2-й моды, и тем более чем для 4-й, что говорит о том, что третья мода в наибольшей степени связанна с частотой дыхания. По этим причинам 3-я мода является контрольной и выбрана для проведения дальнейшего расчета частоты дыхания. Для второй пульсовой волны 4-я мода (PSD[4]) является контрольной, именно она определяет частоту дыхания.

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

Рис. 4. График коррентропии для контрольной моды 1-й волны

Спектральная плотность функции коррентропии вычисляется на основе дискретного преобразования Фурье. На рис. 5 показаны графики спектров для двух пульсовых волн, которые изображены на рис. 1.

Как видно из рис. 5, частота дыхания для первой пульсовой волны равна 0.205 Гц, что соответствует 12.30 вд./мин, а для второй 0.208 Гц –12.48 вд./мин. Абсолютная погрешность относительно контрольного значения частоты составила для первой волны 1.2 вд./мин, а для второй – 0.8 вд./мин.

Рис. 5. Графики спектров функции коррентропии двух пульсовых волн (слева спектр для 1-й волны, справа – для 2-й) с отмеченным значением измеренной частоты дыхания (синим цветом показана кривая сглаженных значений, красным цветом закрашен диапазон частот дыхания D)

В ходе исследования проведена обработка большого набора пульсовых волн, соответствующих мужчинам и женщинам с группировкой по возрастам. При этом получена величина среднего значения абсолютной погрешности², которая для 58.8% всей обработанной выборки не превышало порога в 1 вд./мин., а 95% доверительный интервал для средней абсолютной погрешности всей выборки составил [0.72–2.2] вдохов/мин.

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

Таблица 1. Оценка точности вычисления частоты дыхания от длительности интервала обработки

4. ОБСУЖДЕНИЕ

Пульсовые сигналы могут быть зарегистрированы только на ограниченных интервалах и являются нелинейными и нестационарными. В связи с этим следует отметить значительные преимущества преобразования Гильберта–Хуанга (HHT) для их обработки в сравнении с классическими спектральными методами, основанными на анализе Фурье. В этих условиях спектральный анализ Фурье имеет ограниченное применение.

Значимые преимущества HHT детально обсуждены в статье (Huang et al., 1996). “Основными концептуальными нововведениями являются введение “собственных модовых функций”, основанных на локальных свойствах сигнала, что делает мгновенную частоту значимой; введение мгновенных частот для сложных наборов данных устраняет необходимость в побочных гармониках для представления нелинейных и нестационарных сигналов”.

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

В методе HHT декомпозиция является адаптивным и высокоэффективным процессом, это разложение применимо к нелинейным и нестационарным процессам. На основе преобразования Гильберта вычисляют мгновенные частоты как функции времени, что обеспечивает представление результатов в виде спектра Гильберта, как распределение энергии, частоты, времени.

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

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

Критичным параметром методики является величина длительности интервала обработки пульсовой волны. Анализ результатов влияния длительности интервала пульсовой волны на точность вычисления частоты дыхания показывает, что длительность интервала обработки пульсовой волны не должна быть меньше периода колебания тренда пульсовой волны. Для информации, приведенной в табл. 1, величина периода составляла 30 с. Точность оценки может быть улучшена при осреднении результатов нескольких измерений частоты дыхания на различных интервалах пульсовой волны.

Выбора приоритетных регионов для оценки частоты дыхания основан на вычислении спектрального эксцесса. Известно, что сигнал, имеющий гауссово распределение имеет эксцесс, равный нулю в то время, как при наличии нестационарной составляющей в сигнале, эксцесс становится больше нуля. Достоинством использования спектрального эксцесса являются его высокая чувствительность к изменению спектра сигнала под действием скрытых нестационарностей. Спектральный эксцесс принимает бóльшие значения в полосе частот, где на интервале сигнала присутствуют нестационарности и, равен нулю в полосе частот, где в спектре нестационарности отсутствуют. Однако следует учитывать: если изменения спектра становятся значительными, эксцесс может уменьшается. Чтобы этого избежать, рекомендуется вычислять среднеквадратическое значение эксцесса на исследуемых интервалах и следить за изменениями этой величины.

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

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

5. ВЫВОДЫ

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

Отличительной особенностью методики является применение алгоритма анализа распределения плотности вероятности мгновенных частот в области диапазона дыхания. По результатам применения преобразования Гильберта–Хуанга показано как определять контрольную эмпирическую моду, которая в наибольшей степени связана с искомой частотой дыхания и использовать ее для оценки частоты дыхания.

Экспериментальная проверка разработанной методики на большом объеме данных показала, что среднее значение абсолютной погрешности для 58.8% выборки не превышало порога в 1 вд./мин, а 95% доверительный интервал для средней абсолютной погрешности всей выборки составил [0.72–2.2] вд./мин.

БЛАГОДАРНОСТИ

Автор выражает благодарность сотрудникам лаборатории Российского исследовательского центра LG Electronics Ltd. (г. Москва) за организацию и проведение измерений частоты дыхания.

КОНФЛИКТ ИНТЕРЕСОВ

Автор статьи подтвердил отсутствие конфликта интересов, о котором необходимо сообщить.

ПРИЛОЖЕНИЯ

A. Методика выбора приоритетных регионов для оценки ЧДД

В статье для оценки частоты дыхания использовалась пульсовая волна, снятая с кисти руки (рис. 1). Но для анализа ЧДД могут быть использованы волны, снятые не только с кисти, но и с других участков, например со лба пациента. В этом случае должны быть определены рекомендации по выбору приоритетного региона для снятия пульсовой волны.

В результате исследований этого вопроса предлагается выбирать приоритетные регионы на основе вычисления функции спектрального эксцесса (СЭ) (Nita, 2007). Он является нормированным кумулянтом четвёртого порядка преобразования Фурье и используется для отображения отклонения от гауссовости спектральных составляющих сигнала. СЭ позволяет обнаружить наличие в сигнале скрытых нестационарностей и указать, в каких полосах частот они расположены. Идея предлагаемого метода состоит в сравнении скрытых нестационарностей для нескольких пульсовых волн.

Предлагается следующий алгоритм.

1. Вычислить для каждой k-пульсовой волны значение несмещенной оценки СЭ (Nita, 2007; Vrabie, 2003):

$S{K}_Y\left(f_k\right)=\frac{M}{M-1}\left(\frac{\left(M+1\right){{\displaystyle \sum }}_{i=1}^M{\left(P_{ki}\left(f_k\right)\right)}^4}{{\left({\displaystyle {\sum }_{i=1}^M{P}_{ki}{\left(f_k\right)}^2}\right)}^2}-2\right)$,

где $P_k\left(f_k\right)=\frac{2}{N\sum w_n^2}{\left|X_k\right|}^2$ – оценка спектральной плотность мощности (СПМ) сигнала пульсовой волны, $f_k=\frac{2k}{N}{f}_c$, $f_c=1/2\Delta f$ – частота Найквиста, ∆t – частота дискретизации сигнала, при вычислении СПМ $k=\mathrm{0,}\dots ,N\text{/}2$.

$X_k=\sum _{n=0}^{N-1}{w}_n{x}_n{e}^{-\frac{2\pi ikn}{N}}$ – комплексные коэффициенты дискретного преобразования Фурье (ДПФ), {x_n} – дискретные значения пульсовой волны, {w_n} – набор из N действительных коэффициентов временного окна n (Хеннинга или Хэмминга), при вычислении ДПФ, $k=\mathrm{0,..,}N-\mathrm{1,} i=\mathrm{0,}\dots ,M$, M – общее число спектральных окон.

2. Вычислить площади под кривой СЭ в диапазоне частот дыхания (0.2–0.3) Гц.
3. Вычислить значения площадей (п. 2) для всех пульсовых волн сессии. Найти регион с наименьшими значениями площади. Этот регион и будет приоритетным для вычисления оценки частоты дыхания.
4. При вычислении частоты дыхания применяется усреднение по перекрывающимся сегментам.

B. Критерий определения допустимого натяжения браслета при измерениях

Точность оценки частоты дыхания зависит от натяжения браслета на кисти руки. Влияние величины натяжения браслета исследовалось применительно к двум задачам: оценки ЧДД и исследованию сатурации. Для анализа натяжения браслета предложено использовать значения корреляции Пирсона (Pir) между волнами IR и Red (инфракрасного и красного) излучателей датчика биомониторинга, а также значения корреляции Пирсона между значениями соседних 15 с скользящих окон волн IR и Red (Pir1 и Pir2). При этом критерий основан на оценке не превышения функциями корреляции порогов, соответствующих допустимому натяжению браслета (в данном случае показаны пороги: 0.9, 0.8, 0.8) на рис. 6.

Рис. 6. Графики функций корреляций Пирсона волн IR и Red (n – номер скользящего окна)

В качестве критерия принято отношение суммы значений функции Pir, которые расположены ниже порога (Pir_n) к сумме значений, которые выше порога $Pi{r}_m:Pi{r}_{sumratio}=Pi{r}_n/Pi{r}_m$.

На рис. 7 показаны графики коэффициента натяжения браслета и значения параметра критерия R относительно вычисленных величин сатурации S. Номера у точек от 1 до 4 на рис. 7а соответствуют увеличению натяжения браслета, причем коэффициент натяжения обратно пропорционален величине натяжения. Увеличению натяжения браслета соответствует возрастание значений сатурации.

Рис. 7. Графики коэффициента натяжения браслета и коэффициента корреляции относительно референсного значения (точка 1 соответствует слабому натяжению, точка 4 – сильному)

Достоверные значения R соответствуют вычисленному значению сатурации, которое равно референсному значению сатурации (Ref, %), для рассматриваемого варианта (на графике красная пунктирная линия). На рис. 7б показано поведение функции Pir_sumratio. С увеличением натяжения браслета значения параметра критерия уменьшаются. Выброс точек 1 и 2, находящихся на рисунке в окружностях на рис. 7б, и соответствующие им Pir_{sum_ratio} не соответствуют требуемому натяжению браслета.

В данном случае порогом достоверных значений является величина Pir_sumratio = 0.775 (рис. 7б). Таким образом, для определения допустимого натяжения браслета надо предварительно вычислить корреляции Пирсона для волн IR и Red и значения порогов.

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

² Абсолютная погрешность – разность между вычисленными значениями и контрольными значениями ЧДД, получены с использованием прибора Biopac MP100.

About the authors

P. B. Petrenko

Author for correspondence.
Email: prof.petrenko54@gmail.com
Russian Federation, Saint Petersburg

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Graphs of two pulse waves at a 30 s interval (dotted line shows pulse wave trends)

Download (156KB)

Indexing metadata

3. Fig. 2. Modes (IMD) of the Hilbert-Huang decomposition (pulse wave curve is shown at the top of the graph)

Download (215KB)

Indexing metadata

4. Fig. 3. Probability density plots for two pulse waves, the 2nd, 3rd and 4th modes are shown with the Sk breathing range marked, the 3rd mode corresponds to the curve marked in blue. Graphs (a) correspond to the 1st pulse wave, graphs (b) - to the second wave

Download (276KB)

Indexing metadata

5. Fig. 4. Correntropy plot for the control mode of the 1st wave

Download (102KB)

Indexing metadata

6. Fig. 5. Plots of the spectra of the correntropy function of two pulse waves (left spectrum for the 1st wave, right spectrum for the 2nd wave) with the marked value of the measured respiration rate (blue colour shows the curve of smoothed values, red colour shades the range of respiration rates D)

Download (278KB)

Indexing metadata

7. Fig. 6. Plots of Pearson correlation functions of IR and Red waves (n - number of sliding window)

Download (173KB)

Indexing metadata

8. Fig. 7. Graphs of bracelet tension coefficient and correlation coefficient relative to the reference value (point 1 corresponds to weak tension, point 4 - to strong tension)

Download (200KB)

Indexing metadata