Документация Engee
Notebook

Мониторинг ЭКГ с использованием радара непрерывного излучения и глубокого обучения (Часть 1)

Введение

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

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

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

О данных

В работе в качестве примера используется открытый набор данных, представленный Schellenberger et al. (Scientific Data, 2020)​. Этот датасет содержит ~24 часа синхронных записей радара и эталонных датчиков (ЭКГ, импедансная кардиография, непрерывное АД) от 30 здоровых добровольцев​.
Из датасета был извлечен одно изменерие, над которым проводилось исследование

Импорт библиотек

Импортируем необходимые для работы библиотеки из Julia и Python ( skimage, scipy, neurokit2), а также импортируем в скрипт файл utils.jl, который содержит все необходимые функции

In [ ]: 
include("$(@__DIR__)/install_packages.jl");

Для работы скрипта необходимо установить несколько питоновских пакетов, функции которыъ вызываются через PyCall. Для установки необходимо перейти в командную строку и нажать символ ";". Далее откроется shell, куда нужно прописать:

pip install -r requirements.txt

важно находится в рабочей директории.

In [ ]: 
# ENV["PYTHON"] = '/user/.user_venvs/py3.11/bin'
using PyCall
using Measures
using Peaks
using Plots
using DSP
using Statistics
using Wavelets
skimage = pyimport("skimage");
nk = pyimport("neurokit2")
scipy = pyimport("scipy")
include("$(@__DIR__)/utils.jl");

Загрузка данных

Загрузим данные из файла с помощью функции load_data, где i - синфазная составляющая радарного сигнала, q - квадратурная составляющая радарного сигнала, fs_radar - частота дискретизации радарного сигнала, ecg_data1, ecg_data2 - сигналы ЭКГ, записаные с использованием двух различных электродов, расположенных на теле испытуемого fs_ecg - частота дискретизации сигнала ЭКГ

In [ ]: 
Dirpath = "$(@__DIR__)"
filepath = raw"GDN0011_1_Resting.mat"
order_filter = 6;
full_path = joinpath(Dirpath, filepath);

i, q, fs_radar, ecg_data1, ecg_data2, fs_ecg = load_data(full_path);
cur_ecg = vec(ecg_data2); # сигнал ЭКГ, который будем анализировать

Визуализируем сырые радарные и ЭКГ сигналы

In [ ]: 
p = plot(
  layout        = (4,1),
  size          = (800,800),
  margin        = 5mm,
  bottom_margin = 15mm,     
)

plot!(p, i[1:10000],     subplot=1, label="I",   xlabel="Отсчет",        ylabel="Амплитуда, В")
plot!(p, q[1:10000],     subplot=2, label="Q",   xlabel="Отсчет",        ylabel="Амплитуда, В")
plot!(p, ecg_data1[1:10000], subplot=3, label="1 отведение", xlabel="Отсчет", ylabel="U, мВ")
plot!(p, ecg_data2[1:10000], subplot=4, label="2 отведение", xlabel="Отсчет", ylabel="U, мВ")
Out[0]:

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

In [ ]: 
complex_radar_Data = Complex.(i, q);
sigma = angle.(complex_radar_Data);
lambda_ = 12.5e-3;
range_ = unwrap(sigma) ./ (4 * pi) * lambda_;
plot(layout = (2, 1), size=(600,600),margin=5mm)

plot!(sigma[10000:20000] .* (180/π), subplot=1, label="Разность фаз", xlabel="Отсчет", ylabel="Разность фаз, градусы", legendfontsize = 10)
plot!(range_[10000:20000] .*1e3, subplot=2, label="Смещение грудной клетки", xlabel="Отсчет", ylabel="Расстояние, мм", legendfontsize = 10)
Out[0]:

Анализ сигнала с радара

Функция EllipseReconstruction(i, q) нормализует IQ-сигнал радара: аппроксимирует его эллипсом, центрирует, поворачивает и масштабирует в окружность, после чего извлекает фазу и переводит её в микросмещения цели. Это позволяет получить чистый сигнал перемещения, из которого затем выделяются дыхание и пульс c помощью функции getVitalSignals.

In [ ]: 
I_comp, Q_comp, sigma, range_ = EllipseReconstruction(i, q);
radar_respiration, radar_pulse, radar_heartsound, radar_heartsound_denoise = getVitalSignals(range_, fs_radar, order_filter);

Поскольку сигналы достаточно большие, для визуальной оценки визуализируем извлеченыне пульс, дыхание, сердченые звуки на интервале [1, 10000]

In [ ]: 
plot(layout = (2, 2), size=(800,500), margin=5mm)

plot!(radar_respiration[1:10000], subplot=1, label="Дыхание", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_pulse[1:10000], subplot=2, label="Пульс", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_heartsound[1:10000], subplot=3, label="Звуки сердца", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_heartsound_denoise[1:10000], subplot=4, label="Звуки сердца после фильтра", xlabel="Отсчет", ylabel="Амплитуда, мВ")
Out[0]:

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

Поиск R пиков на радарном сигнале

Далее вычислим ускорение движения грудной клетки, отфильтруем его на полосе [4, 20] Гц

In [ ]: 
acc = diff(diff(radar_pulse)) * fs_radar^2;
radar_pulse_acc = Filter(acc, 4, 20, fs_radar, 6);
time_acc = range(0, step=1/fs_radar, length=length(radar_pulse_acc));

Визуализируем сигнал ускорения грудной клетки на интервале [10000, 40000] Гц

In [ ]: 
plot(time_acc[10000:40000], radar_pulse_acc[10000:40000], xlabel="Время (с)", ylabel="Ускорение, м/с^2", title="Ускорение грудной клетки", label="")
Out[0]:

Далее из радарного пульсового сигнала вычисляется амплитудная огибающая с помощью преобразования Гильберта, после чего в ней ищутся сердечные пики — только те, что достаточно высокие (выше 20% медианы) и разделены минимум на 0.76 секунды.

In [ ]: 
env = abs.(hilbert(radar_pulse_acc));
radar_peacks, _ = scipy.signal.find_peaks(env,
                      distance=Int(fs_radar * 0.76), height=median(env) * 0.2);

Отобразим фрагмент радарного сигнала: чёрная линия показывает огибающую, оранжевая — пульсовую компоненту (ускорение грудной клетки), а красные точки отмечают найденные R-пики в пределах выделенного временного окна.

In [ ]: 
time_start = 80000; time_end = 100000;
r_subset = filter(x -> time_start ≤ x ≤ time_end, radar_peacks)
radar_part = env[time_start:time_end]
time_part = time_acc[time_start:time_end]

time_peaks = time_acc[r_subset]

plot(time_part, radar_part;
    color = :black,
    lw = 1.5,
    label = "Огибающая",
    xlabel = "Время, с",
    ylabel = "Амплитуда, м/с^2",
    title = "Полный фрагмент"
)

plot!(time_part, radar_pulse_acc[time_start:time_end];
    label = "Ускорение грудной клетки"
)

scatter!(time_peaks, env[r_subset];
    marker = (:circle, 6),
    mc = :red,
    label = "R-пики"
)
Out[0]:

Анализ ЭКГ

Далее проанализируем имеющийся ЭКГ сигнал

Из ЭКГ-сигнала сначала удаляются шумы (ecg_clean), затем находятся положения R-пиков (ecg_peaks), после чего выполняется поиск пиков P, Q, S и T с помощью метода peak

In [ ]: 
ecg_cleaned = nk.ecg_clean(cur_ecg, sampling_rate=fs_ecg);
signals, rpeaks = nk.ecg_peaks(ecg_cleaned, sampling_rate=fs_ecg);
_, waves = nk.ecg_delineate(ecg_cleaned, rpeaks, sampling_rate=fs_ecg, method="peak");

Визуализируем участок ЭКГ с пиками P, Q, R, S, T волн с помощью кастомной функции, определенной в файле utils.py

In [ ]: 
plot_ecg_with_peaks(ecg_cleaned, waves, rpeaks; range_start=90000, range_end=100000)
Out[0]:
In [ ]: 
rr_intervals_ecg = diff(rpeaks["ECG_R_Peaks"]) / fs_ecg;
times_ecg = cumsum([0.0; rr_intervals_ecg]);

rr_intervals_radar = diff(radar_peacks) / fs_radar;
times_radar = cumsum([0.0; rr_intervals_radar]);

scatter(
    times_radar[2:end], rr_intervals_radar;
    label   = "RR радар",
    xlabel  = "Отсчет",
    ylabel  = "RR",
    title   = "RR радар & ЭКГ",
    color   = :green
)

scatter!(
    times_ecg[2:end], rr_intervals_ecg;
    label = "RR ЭКГ",
    color   = :red
)
Out[0]:
In [ ]: 
print("Кол-во RR интервалов на радарном сигнале: ", length(rr_intervals_radar), "\n")
print("Кол-во RR интервалов на ЭКГ: ", length(rr_intervals_ecg))

Кол-во RR интервалов на радарном сигнале: 621
Кол-во RR интервалов на ЭКГ: 659

На ритмограмме красные точки (RR по ЭКГ) и зелёные (RR по радару) в основном лежат вплотную друг к другу вокруг, что говорит об общей хорошей сходимости интервалов. При этом у радара встречаются рассеянные выбросы вверх (до ~1.9 с), а у ЭКГ — вниз (до ~0.6 с), что указывает на редкие пропуски или ложные срабатывания при детекции. В средней части (примерно от 100 до 400 отсчёта) кластеры особенно плотные и согласованы. В целом между RR_radar и RR_ECG наблюдается высокая согласованность

Далее рассчитаем мгновенную частоту сердечных сокращений для радарных и ЭКГ-данных, а затем из полученных HR вычисляются базовые статистики — среднее, медиана, минимальное и максимальное значения, что позволяет сравнить усреднённый пульс и разброс ЧСС в обеих методиках.

In [ ]: 
HR_radar = 60 .÷ rr_intervals_radar;
HR_ecg = 60 .÷ rr_intervals_ecg;
mean_HRRadar   = mean(HR_ecg)
median_HRRadar = median(HR_ecg)
min_HRRadar    = minimum(HR_ecg)
max_HRRadar    = maximum(HR_ecg)

mean_HRECG   = mean(HR_radar)
median_HRECG = median(HR_radar)
min_HRECG    = minimum(HR_radar)
max_HRECG    = maximum(HR_radar);
In [ ]: 
println("Статистика HR (уд/мин)")
println("──────────────────────────────")
println("             |     ECG     |   Radar")
println("─────────────┼─────────────┼────────────")
println("Mean HR      | ", round(mean_HRECG, digits=2), "       | ", round(mean_HRRadar, digits=2))
println("Median HR    | ", round(median_HRECG, digits=2), "        | ", round(median_HRRadar, digits=2))
println("Min HR       | ", round(min_HRECG, digits=2), "        | ", round(min_HRRadar, digits=2))
println("Max HR       | ", round(max_HRECG, digits=2), "        | ", round(max_HRRadar, digits=2))

Статистика HR (уд/мин)
──────────────────────────────
             |     ECG     |   Radar
─────────────┼─────────────┼────────────
Mean HR      | 58.18       | 62.1
Median HR    | 58.0        | 60.0
Min HR       | 29.0        | 42.0
Max HR       | 77.0        | 112.0

В целом согласованность неплохая — средние и медианные значения ЧСС близки (≈58 vs ≈62 уд/мин), но у радарного метода заметно больше разброс: минимумы и максимумы значительно выше, что говорит о ложных пиках или пропусках в детекции.

Анализ вариабельности сердечного ритма (HRV)

Можно получить полный набор показателей HRV — от простых статистических метрик до продвинутых нелинейных и мульти–фрактальных индексов с помощью функции hrv из библиотеки neuroit2.

In [ ]: 
hrv_indices = nk.hrv(rpeaks, sampling_rate=fs_ecg, show=false)
Out[0]:
HRV_MeanNN HRV_SDNN HRV_SDANN1 HRV_SDNNI1 HRV_SDANN2 HRV_SDNNI2 HRV_SDANN5 HRV_SDNNI5 HRV_RMSSD HRV_SDSD HRV_CVNN HRV_CVSD HRV_MedianNN HRV_MadNN HRV_MCVNN HRV_IQRNN HRV_SDRMSSD HRV_Prc20NN HRV_Prc80NN HRV_pNN50 HRV_pNN20 HRV_MinNN HRV_MaxNN HRV_HTI HRV_TINN HRV_ULF HRV_VLF HRV_LF HRV_HF HRV_VHF HRV_TP HRV_LFHF HRV_LFn HRV_HFn HRV_LnHF HRV_SD1 HRV_SD2 HRV_SD1SD2 HRV_S HRV_CSI ... HRV_C1a HRV_SD1d HRV_SD1a HRV_C2d HRV_C2a HRV_SD2d HRV_SD2a HRV_Cd HRV_Ca HRV_SDNNd HRV_SDNNa HRV_DFA_alpha1 HRV_MFDFA_alpha1_Width HRV_MFDFA_alpha1_Peak HRV_MFDFA_alpha1_Mean HRV_MFDFA_alpha1_Max HRV_MFDFA_alpha1_Delta HRV_MFDFA_alpha1_Asymmetry HRV_MFDFA_alpha1_Fluctuation HRV_MFDFA_alpha1_Increment HRV_DFA_alpha2 HRV_MFDFA_alpha2_Width HRV_MFDFA_alpha2_Peak HRV_MFDFA_alpha2_Mean HRV_MFDFA_alpha2_Max HRV_MFDFA_alpha2_Delta HRV_MFDFA_alpha2_Asymmetry HRV_MFDFA_alpha2_Fluctuation HRV_MFDFA_alpha2_Increment HRV_ApEn HRV_SampEn HRV_ShanEn HRV_FuzzyEn HRV_MSEn HRV_CMSEn HRV_RCMSEn HRV_CD HRV_HFD HRV_KFD HRV_LZC
0 982.672231 138.257783 17.203166 134.240934 14.416785 136.780236 NaN NaN 174.439493 174.572008 0.140696 0.177515 999.5 105.2646 0.105317 140.25 0.792583 909.5 1080.0 63.429439 84.218513 532.0 1420.0 26.36 445.3125 NaN 0.011477 0.036597 0.070714 0.007943 0.126731 0.517539 0.288778 0.557983 -2.649114 123.441051 151.645602 0.81401 58808.391731 1.228486 ... 0.402532 95.415128 78.31781 0.29099 0.70901 81.802871 127.689776 0.413133 0.586867 88.870007 105.920627 0.731215 2.380941 1.043536 1.665046 -1.553554 -2.035018 -0.238965 0.003628 0.44613 0.526557 0.592637 0.639657 0.637235 0.591616 0.056814 -0.504086 0.000063 0.020011 1.361697 1.581122 8.666462 1.230945 0.872857 1.439033 1.934975 1.740408 1.928947 3.763396 0.952044

1 rows × 91 columns

Визуализируем несколько полученных статистик

In [ ]: 
time_keys  = ["HRV_SDNN", "HRV_RMSSD", "HRV_SDSD", "HRV_pNN50", "HRV_pNN20"]
freq_keys  = ["HRV_VLF", "HRV_LF", "HRV_HF", "HRV_LFHF"]
poinc_keys = ["HRV_SD1", "HRV_SD2"]
geom_keys  = ["HRV_TINN", "HRV_HTI"]
ent_keys   = ["HRV_SampEn", "HRV_ApEn"]


getval(key) = Float64(PyCall.getindex(hrv_indices[key], 1))


time_vals  = [getval(k) for k in time_keys]
freq_vals  = [getval(k) for k in freq_keys]
poinc_vals = [getval(k) for k in poinc_keys]
geom_vals  = [getval(k) for k in geom_keys]
ent_vals   = [getval(k) for k in ent_keys]


p = plot(layout = (1, 2), size=(500, 500), theme=:default, framestyle=:box, grid=true)


bar!(
    p[1, 1],
    1:length(time_keys), time_vals;
    xticks=(1:length(time_keys), time_keys),
    rotation=45,
    xlabel="Временной домен",
    title="Временной HRV",
    legend=false
)


bar!(
    p[1, 2],
    1:length(freq_keys), freq_vals;
    xticks=(1:length(freq_keys), freq_keys),
    rotation=45,
    xlabel="Частотный домен", 
    title="Частотный HRV",
    legend=false)
Out[0]:

Выводы

В этом демо-примере был разобран способ анализа ЭКГ сигналов и предобработка радарных данных