使用连续辐射雷达和深度学习的心电监测（第二部分）

心电图用于监测心脏：电极放置在皮肤上并记录心脏的活动。然而，由于记录的持续时间短、贴纸的不适感和在移动过程中出现伪影的风险，这种方法并不总是方便的。

雷达系统提供非接触式监测：它们从心跳和呼吸中检测胸部的微运动，而不限制患者的运动，也不需要皮肤接触。使用神经网络，可以从雷达信号的机械振动重建心电图，将胸部的运动转变为通常的P-QRS-T波。

作为演示示例的一部分，我们将重点介绍用于心脏监测的生物雷达以及使用AI算法将雷达信号转换为ECG的等效信号。

要了解该项目的第一部分，专门用于ECG信号特征和反射回波雷达的分析和比较，请转到ссылке.

工作准备

导入必要的库

我们从Julia和Python（skimage，scipy，neurokit2）导入必要的库，并导入utils。jl文件到脚本，其中包含所有必要的功能。

include("$(@__DIR__)/install_packages.jl");

为了使脚本工作，您需要安装几个Python包，其功能通过PyCall调用。要安装，请转到命令提示符并单击";"符号。接下来，shell打开，您需要在其中注册:

pip安装-r要求.txt脧脗脭脴

重要的一个位于工作目录中。

using Glob
using DSP
using CUDA
using Flux 
using Statistics
using cuDNN
using BSON: @save, @load

using Measures

由于使用的函数具有相当多的代码行数，因此它们已被移动到单独的文件中。下面导入的jls

的模型。jl文件实现所用模型的结构。的数据集。jl文件用于实现数据集模块的逻辑，就绪数据。jl包含为训练准备数据的函数和utils。jl文件包含项目工作所需的所有其他功能。

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小时同步记录。从正在调查的数据集中提取了一个变异

创建一个文件夹，其中信号处理的结果将被保存。

prepareRadarData = joinpath(Dirpath, "prepareRadarData")
isdir(prepareRadarData) || (mkdir(prepareRadarData); println("Папка prepareRadarData создана"));

接下来，我们将对演示示例中执行的原始ECG和雷达数据执行所有转换，以分析这些数据，即:

椭圆的重建
重要信号的提取-呼吸，脉搏
计算胸部加速度和随后的过滤
心电图过滤

该功能的实现位于就绪。jl文件。

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的负载，并且使用重叠允许您考虑时间依赖性。

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);

让我们看看心电图和雷达信号是什么样子，切成重叠的窗口。

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的访问，增加将加速几个数量级。

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文件。初始化模型并将其传输到正在使用的设备

模型、优化器、损失函数和其他参数的初始化

model = Radar2ECG_UNet_SE(32);
device == gpu ? model = gpu(model) : nothing;

让我们定义一些参数：学习率，epochs的数量和模型的名称需要在工作区中自动创建一个目录，其中模型的权重和各种度量的图表将在学习过程中存储。

我们还将初始化优化器，loss函数，它在utils中实现。jl文件，并使用辅助函数get_training_dir，它将创建必要的目录来保存学习结果。

fs_ecg = 2000;

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)

/user/Диплом/AI/model/trainedModel/train1 train1

让我们定义计算指标的数据将存储的列表。在这个例子中，我们根据tau水平计算MAE，RMSE，R2，精度。

我们还初始化了params变量，它存储了训练周期的一些参数。

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文件。完成培训后，结果将保存到创建的文件夹中

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

基于测试数据的推理模型

让我们检查模型如何在测试数据上工作。

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)

训练模型的评估

让我们在列车和验证上显示错误图。

plot(train_losses, label="train loss", lw=3)
plot!(valid_losses, label="valid loss", lw = 3)
xlabel!("Эпоха")
ylabel!("Ошибка")

我们可视化训练过程中积累的度量-R^2，RMSE，MAE，Accuracy

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类神经网络，用于心电重建。该模型在计算指标上显示良好的数字方面做得很好。