使用连续辐射雷达和深度学习的心电监测（第二部分）
心电图用于监测心脏：电极放置在皮肤上并记录心脏的活动。 然而，由于记录的持续时间短、贴纸的不适感和在移动过程中出现伪影的风险，这种方法并不总是方便的。
雷达系统提供非接触式监测：它们从心跳和呼吸中检测胸部的微运动，而不限制患者的运动，也不需要皮肤接触。 使用神经网络，可以从雷达信号的机械振动重建心电图，将胸部的运动转变为通常的P-QRS-T波。
作为演示示例的一部分，我们将重点介绍用于心脏监测的生物雷达以及使用AI算法将雷达信号转换为ECG的等效信号。
要了解该项目的第一部分，专门用于ECG信号特征和反射回波雷达的分析和比较，请转到ссылке.
工作准备
导入必要的库
我们从Julia和Python（skimage，scipy，neurokit2）导入必要的库，并导入utils。jl文件到脚本，其中包含所有必要的功能。
In [ ]:
include("$(@__DIR__)/install_packages.jl");
为了使脚本工作，您需要安装几个Python包，其功能通过PyCall调用。 要安装，请转到命令提示符并单击";"符号。 接下来，shell打开，您需要在其中注册:
pip安装-r要求.txt脧脗脭脴
重要的一个位于工作目录中。
In [ ]:
using Glob
using DSP
using CUDA
using Flux
using Statistics
using cuDNN
using BSON: @save, @load
In [ ]:
using Measures
由于使用的函数具有相当多的代码行数，因此它们已被移动到单独的文件中。下面导入的jls
的模型。jl文件实现所用模型的结构。 的数据集。jl文件用于实现数据集模块的逻辑，就绪数据。jl包含为训练准备数据的函数和utils。jl文件包含项目工作所需的所有其他功能。
In [ ]:
Dirpath = "$(@__DIR__)"
include(joinpath(Dirpath, "model.jl"))
include(joinpath(Dirpath, "dataset.jl"))
include(joinpath(Dirpath, "prepare_data.jl"))
include(joinpath(Dirpath, "utils.jl"));
数据准备
该论文使用Schellenberger等人提出的开放数据集。 （科学数据，2020）。 该数据集包含来自30名健康志愿者的雷达和参考传感器（ECG，阻抗心电图，连续BP）的~24小时同步记录。 从正在调查的数据集中提取了一个变异
创建一个文件夹，其中信号处理的结果将被保存。
In [ ]:
prepareRadarData = joinpath(Dirpath, "prepareRadarData")
isdir(prepareRadarData) || (mkdir(prepareRadarData); println("Папка prepareRadarData создана"));
接下来，我们将对演示示例中执行的原始ECG和雷达数据执行所有转换，以分析这些数据，即:
- 椭圆的重建
- 重要信号的提取-呼吸，脉搏
- 计算胸部加速度和随后的过滤
- 心电图过滤
该功能的实现位于就绪。jl文件。
In [ ]:
root = joinpath(Dirpath, "data")
pat = r"^GDN\d{4}_.+\.mat$"
order_filter = 6
files = [p for p in glob("**/*.mat", root) if occursin(pat, basename(p))]
save_path = joinpath(Dirpath, "prepareRadarData")
isdir(save_path) || mkdir(save_path)
for file in files
prepareData(file, save_path, order_filter)
end
使用上面的块，我们处理了数据并将其保存在一个单独的文件夹中。
创建数据加载器
接下来，我们将创建一个数据加载器，它将批量加载数据到模型中。 在下面的块中，我们将初始化数据集，将数据划分为训练集、验证集和测试集，并创建数据加载器。
数据集逻辑写在数据集中。jl文件
数据集将ECG和雷达信号切片到复盖的窗口上。 切片允许您减少由于使用长序列而对ANN的负载，并且使用重叠允许您考虑时间依赖性。
In [ ]:
ds = RadarECG("prepareRadarData");
N = length(ds)
train_len = floor(Int, 0.7 * N)
val_len = floor(Int, 0.15 * N)
test_len = N - train_len - val_len
train_idx = 1:train_len
val_idx = (train_len+1):(train_len+val_len)
test_idx = (train_len+val_len+1):N
X_train = [ ds[i][1] for i in train_idx ]
Y_train = [ ds[i][2] for i in train_idx ]
X_val = [ ds[i][1] for i in val_idx ]
Y_val = [ ds[i][2] for i in val_idx ]
X_test = [ ds[i][1] for i in test_idx ]
Y_test = [ ds[i][2] for i in test_idx ]
batch_size = 64
train_loader = Flux.DataLoader((X_train, Y_train);
batchsize = batch_size,
shuffle = false,
partial = true,
parallel = true,
collate = mycollate
)
val_loader = Flux.DataLoader((X_val, Y_val); batchsize=batch_size, shuffle=true, parallel = true, collate = mycollate)
test_loader= Flux.DataLoader((X_test, Y_test); batchsize=batch_size, shuffle=true, parallel = true, collate = mycollate);
让我们看看心电图和雷达信号是什么样子，切成重叠的窗口。
In [ ]:
radar_batch, ecg_batch = first(train_loader)
radar_batch = cpu(radar_batch)
ecg_batch = cpu(ecg_batch)
win, B = size(radar_batch)
p = plot(size=(500,500), layout = (2,1), margin=15mm)
for (arr, name) in zip((radar_batch, ecg_batch), ("radar", "ecg"))
target = name == "radar" ? 1 : 2
for j in 0:3
x = (1:win) .+ j*ds.step
y = arr[:, j+1]
if target == "radar"
plot!(p, x, y; subplot=target, lw=2, showlegend=false, alpha = 0.7, xlabel="Отсчет", ylabel="Ускорение, м/с^2")
else
plot!(p, x, y; subplot=target, lw=2, showlegend=false, alpha = 0.7, xlabel="Отсчет", ylabel="Амплитуда, мВ")
end
end
end
title!(p, "Радар", subplot=1)
title!(p, "ЭКГ", subplot=2)
display(p)
之后，我们将确定模型将被训练的设备。 如果有对GPU的访问，增加将加速几个数量级。
In [ ]:
device = CUDA.functional() ? gpu : cpu;
device == gpu ? train_loader = gpu.(train_loader) : nothing
device == gpu ? val_loader = gpu.(val_loader) : nothing
device == gpu ? test_loader = gpu.(test_loader) : nothing;
模型的架构呈现在模型中。jl文件。 初始化模型并将其传输到正在使用的设备
模型、优化器、损失函数和其他参数的初始化
In [ ]:
model = Radar2ECG_UNet_SE(32);
device == gpu ? model = gpu(model) : nothing;
让我们定义一些参数：学习率，epochs的数量和模型的名称需要在工作区中自动创建一个目录，其中模型的权重和各种度量的图表将在学习过程中存储。
我们还将初始化优化器，loss函数，它在utils中实现。jl文件，并使用辅助函数get_training_dir，它将创建必要的目录来保存学习结果。
In [ ]:
fs_ecg = 2000;
In [ ]:
learning_rate = 1e-4
num_epochs = 100
model_name = "trainedModel"
lossFunction = loss;
# opt = Flux.Adam(learning_rate)
opt = Flux.Optimiser(WeightDecay(1e-6), Adam(learning_rate))
epoch_number = 0
train_dir, name_new_dir = get_training_dir(Dirpath, model_name)
print(train_dir, " ", name_new_dir)
让我们定义计算指标的数据将存储的列表。 在这个例子中，我们根据tau水平计算MAE，RMSE，R2，精度。
我们还初始化了params变量，它存储了训练周期的一些参数。
In [ ]:
tau = 0.05
best_val_loss = Inf
train_losses = Float64[]
valid_losses = Float64[]
MAEList = Float64[]
RMSEList = Float64[]
R2List = Float64[]
AccList = Float64[]
params = Dict("epochs"=> num_epochs, "batch_size"=>batch_size, "learning_rate"=> learning_rate, "tau"=>tau, "model"=>"Radar2ECG_UNet", "info"=>"Train model without scheduler." );
模型训练
接下来，我们将开始模型的训练周期。 Train和validate函数位于utils中。jl文件。 完成培训后，结果将保存到创建的文件夹中
In [ ]:
no_improve_epochs = 0
best_model = nothing
try
for epoch in 1:num_epochs
println("-"^50 * "\n")
println("EPOCH $(epoch):")
opt, train_loss = train!(
model, train_loader, opt,
lossFunction, device, 1, num_epochs
)
val_loss, metrics = validate(model, val_loader, lossFunction, gpu)
@show "R2" metrics.R2
if val_loss < best_val_loss
best_val_loss = val_loss
no_improve_epochs = 0
else
no_improve_epochs += 1
if no_improve_epochs >= 10
old_lr = opt[2].eta
opt[2].eta *= 0.9
@info "LR reduced: $old_lr → $(opt[2].eta) at epoch $epoch"
no_improve_epochs = 0
end
end
println("Epoch $epoch/$num_epochs | train $(round(train_loss, digits=4)) | val $(round(val_loss, digits=4))")
println("Current LR = $(opt[2].eta)")
push!(train_losses, train_loss)
push!(valid_losses, val_loss)
push!(MAEList, metrics.MAE)
push!(RMSEList, metrics.RMSE)
push!(R2List, metrics.R2)
push!(AccList, metrics[Symbol("Acc@0.05")])
if val_loss == best_val_loss
best_model = deepcopy(model)
@info "Валидационная потеря упала" val_loss=best_val_loss
end
end
catch e
if e isa InterruptException
println("Training interrupted by user.")
else
rethrow(e)
end
finally
best_model = cpu(best_model)
output_dir_models = joinpath(train_dir, "models_bson")
isdir(output_dir_models) || mkdir(output_dir_models)
best_path = joinpath(output_dir_models, "model_best.bson")
@info "Saving best model (val_loss=$(best_val_loss)) to $best_path"
@save best_path best_model
save_plots!(MAEList, RMSEList, R2List, AccList,
train_losses, valid_losses;
model_dir = "model",
run_name = name_new_dir,
plot_name = model_name)
params_path = joinpath(train_dir, "params.json")
open(params_path, "w") do io
JSON.print(io, params)
end
println("All done. Artifacts saved in $train_dir")
end
基于测试数据的推理模型
让我们检查模型如何在测试数据上工作。
In [ ]:
radar, ecg = test_loader[2]
# radar, ecg = batch
radar_ = reshape(radar, size(radar,1), 1, size(radar,2)) |> device
ecg_gpu = ecg |> device
pred_gpu = model(radar_)
radar_cpu = Array(radar)
pred_cpu = Array(dropdims(pred_gpu; dims=2))
ecg_cpu = Array(ecg_gpu)
B, L = size(radar_cpu)
rows, cols = 4, 2
nplots = min(rows*cols, L) # вдруг каналов < 8
plt = plot(layout = (rows, cols), # 4×2
size = (1500, 1500),
margin = 15mm,
legend = :outertopright,
legendfontsize = 8)
for i in 1:nplots
# подписи в легенде один-раз, чтобы не плодить дубликаты
lbls = i == 1 ? ("Радар", "Восстановленная ECG", "Эталонная ECG") : ("", "", "")
plot!(plt, radar_cpu[:, i + 2], lw = 2, alpha = 0.7,
label = lbls[1], subplot = i)
plot!(plt, pred_cpu[:, i + 2], lw = 2,
label = lbls[2], subplot = i)
plot!(plt, ecg_cpu[:, i + 2], lw = 2, linestyle = :dash,
label = lbls[3], subplot = i)
title!(plt, "Пример $i", subplot = i)
xlabel!(plt, "Время [s]", subplot = i)
ylabel!(plt, "Амплитуда, мВ", subplot = i)
end
display(plt)
训练模型的评估
让我们在列车和验证上显示错误图。
In [ ]:
plot(train_losses, label="train loss", lw=3)
plot!(valid_losses, label="valid loss", lw = 3)
xlabel!("Эпоха")
ylabel!("Ошибка")
Out[0]:
我们可视化训练过程中积累的度量-R^2，RMSE，MAE，Accuracy
In [ ]:
plt = plot(layout = (1,4), size = (900,300),
legend = :outertopright,
legendfontsize = 8,
margin = 5Plots.mm)
metricsList = [RMSEList, MAEList, R2List, AccList]
metric_names = ["RMSE", "MAE", "R²", "Accuracy"]
plt = plot(
layout = (1, 4),
size = (1000, 300),
margin = 5Plots.mm,
legend = false,
)
for (i, name) in enumerate(metric_names)
plot!(
plt, metricsList[i], subplot = i,lw = 2,alpha = 0.8,label = "")
title!(plt, name, subplot = i)
xlabel!(plt, "Epoch", subplot = i)
ylabel!(plt, name, subplot = i)
end
display(plt)
结论
在工作中，在预处理雷达数据上训练了一个自定义的U-net类神经网络，用于心电重建。 该模型在计算指标上显示良好的数字方面做得很好。