数字血压计的模拟和人体血压的分析

目前，有两种最常见的测量血压的方法-Korotkov方法（听诊或声学）和oscillometric。这些方法中的每一种都提供了必要的准确性和方便的再现性，然而，它们中的每一种都有其优点和缺点。

目前，Korotkov方法是世界卫生组织批准的非侵入性血压测量方法。使用眼压计进行测量，并且使用听诊器或麦克风收听来自袖带夹住的动脉的音调。

示波测量方法与Korotkov方法一样，基于记录层流血流变为湍流时的动脉搏动，但不听音调。振荡测量方法使用直接在袖带中记录压力波动（振荡）。因此，在第一个Korotkov音调期间，观察到脉动幅度最显着的增加–记录收缩压。脉动幅度的急剧减小表明湍流血流到层流-舒张压的变化。

使用此示例，您可以使用测量压力的示波法研究数字眼压计的基本原理，并在研究期间考虑人体运动的影响。

安装必要的库以进行进一步建模。

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["Random","Plots", "Statistics", "DSP"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

using Plots, Random, DSP, Statistics

关于心血管系统结构的一般信息

首先，我们将描述一个人（病人）的生理参数，这些参数是建模所必需的

*脉搏率-心脏的收缩（收缩）和松弛（舒张）的周期数，通过每单位时间外周动脉中的脉搏波记录

***收缩压（上）动脉压（SD）**是心脏最大收缩时刻的血压水平，表征左心室心肌的状态。

***舒张（较低）血压（DB）**是心脏最大放松时刻的血压水平，表征动脉壁的音调程度。

# 人体生理参数（paceinta）
age = 20            # 年龄、年份
weight = 80         # 重量，公斤
height = 180        # 身高，厘米

pulse_rate = 75;             # 脉搏率，节拍/分钟。
systolic_bp = 120.0;         # 收缩压，mmHg
diastolic_bp = 80.0;         # 舒张压，mmHg

血压以毫米汞柱为单位测量，缩写为mmHg。血压值为120/80意味着收缩压为120mmHg，舒张压为80mmHg。

收缩压和舒张压值之间的差异称为脉压（PD）。它显示了收缩压超过舒张压的程度，这是在收缩压期间打开半月板主动脉瓣所必需的。通常，脉压为35-55毫米汞柱。

动脉压（BP）以毫米汞柱（mmHg）测量。
写入120/80mmHg的血压值意味着:

*收缩压（DM）=120mmHg。

*舒张压（DD）=80mmHg。

脉压（PD）是收缩压和舒张压值之间的差异。:

pulse_bp = systolic_bp - diastolic_bp;      # 脉冲压力，mmHg

为了确定平均动脉压（MPP），它表示持续血流的能量，可以使用以下Hickam公式:

mean_bp = (systolic_bp + 2*diastolic_bp)/3;   # 平均血压，mmHg

为了评估心血管系统的功能状态，计算心脏的分钟体积（MO）并与适当值（DMO）进行比较。

其中2.2是心脏指数，以升为单位，PT是体表，使用nomogram或各种公式计算。在这个例子中，将使用Dubois公式。:

并且心脏的中风体积（UO）由表达式计算:

其中B是年龄，g。

让我们来介绍一下英文的符号

*PT(体表)→BSA(体表面积)

*UO(冲击体积)→SV(冲程体积)

*MO（分钟体积）→CO（心输出量）

*DMO(应有分钟体积)→RCO(所需心输出量)

BSA =  0.007184 * (weight^0.425) * (height^0.725);                  # 体表面，m^2
SV = 101 + 0.5 * systolic_bp - 1.08 * diastolic_bp - 0.6 * age;     # 冲击体积，ml
CO = SV * pulse_rate;                                               # 分钟体积，ml/min
RCO = 2.2 * BSA * 1000;                                             # 所需的分钟体积，ml/min

我们将显示计算心血管系统（CVS）参数的所得结果。

println("计算心血管系统（CVS）参数的结果:")
println("----------------------------------------------------------------")
println("1. 脉搏压力(PP): ", round(pulse_bp; digits=1), " mmHg")
println("2. 平均动脉压(MPP): ", round(mean_bp; digits=1), " mmHg")
println("3. 体表面积(PT): ", round(BSA; digits=4), " m2")
println("4. 冲击体积(UO): ", round(SV; digits=1), " 毫升")
println("5. 分钟心脏体积(MO): ", round(CO; digits=1), " 毫升/分钟")
println("6. 适当的分钟体积(DMO): ", round(RCO; digits=1), " 毫升/分钟")

Результаты расчета параметров сердечно сосудистой системы (ССС):
----------------------------------------------------------------
1. Пульсовое давление (ПП): 40.0 мм.рт.ст
2. Среднее артериально давление (СрАД): 93.3 мм.рт.ст
3. Площадь поверхности тела (ПТ): 1.9964 м²
4. Ударный объем (УО): 62.6 мл
5. Минутный объем сердца (МО): 4695.0 мл/мин
6. Должный минутный объем (ДМО): 4392.1 мл/мин

用振荡测量方法模拟数字眼压计的操作

具有示波测量方法的数字眼压计通过分析脉搏波从动脉传输时袖带中发生的气压波动（振荡）来确定血压。让我们按顺序分析他的工作阶段。:

增加袖带中的压力
袖带中的压力缓慢释放
波动登记
测量结果的处理
显示结果

建模的准备阶段包括设置后续生成生理数据所必需的信号参数。这个阶段决定了所有信号将建立在其上的时间轴的关键特性。

# 信号参数
fs = 200            # 采样率(200赫兹)
total = 60.0        # 总测量持续时间（60秒）

# 时间表
t = 0:1/fs:total - 1/fs;

第一阶段是模拟眼压计袖带中的压力。该模型使用三个关键参数生成逼真的压力曲线。:

*抽水时间:10.0秒
（舒适的袖带填充的最佳时间，而不会对动脉造成过度负荷）

*高原持续时间:1.0秒
（动脉完全压缩的压力稳定期）

*峰值压力:160.0mmHg
(标准值超过典型收缩压30-40毫米汞柱)

# 眼压计袖带压力测量的阶段
time_pump = 10.0                # 袖带充气的持续时间，秒
time_plato = 1.0                # 泵送后的高原持续时间，秒
peak_pressure = 160.0           # 袖带内最大压力，mmHg


# 眼压计袖带压力模拟
function generate_pressure(time)
    pressure = zeros(length(time))       # 创建一个与时间向量长度相同的零数组
    
    # 将压力泵入袖带
    pump_phase = time .<= time_pump
    pressure[pump_phase] = peak_pressure * (time[pump_phase] / time_pump).^0.85

    # 高原
    plato_start = time_pump
    plato_end = time_pump + time_plato
    plato_phase = (time .> plato_start) .& (time .<= plato_end)
    pressure[plato_phase] .= peak_pressure

    # 袖带减压
    deflate_phase = (time.> plato_end)
    deflate_time = time[deflate_phase] .- plato_end
    tau = 7.0 
    pressure[deflate_phase] = peak_pressure .* exp.(-deflate_time ./ tau)
    
    return pressure, deflate_phase, plato_end
end;

第二阶段包括对眼压计袖带中的振荡进行建模。 Generart_oscillation函数考虑到以下几个方面来实现波动:

*收缩期对应于直接脉冲节拍。

*Decrotic衰退模拟来自闭合主动脉瓣的反射波。

*舒张波动反映血管中的残余搏动。

*高斯调制模拟真实现象：当袖带中的压力等于平均动脉压时，观察到振荡的最大振幅。

function generate_oscillations(pressure, time, pulse_rate, systolic_bp, diastolic_bp, mean_bp, inflation_duration; movement=false)
    pulse_freq = pulse_rate / 60.0
    T = 1 / pulse_freq
    oscillations = zeros(length(time))
    movement_detected = false
    movement_start_time = 0
    movement_end_time = 0
    movement_start = 140.0          # 运动伪影的初始压力，mmHg
    movement_end = 110.0            # 运动伪影的最终压力，mmHg
    
    for i in eachindex(time)
        t_val = time[i]
        p = pressure[i]
        t_mod = mod(t_val, T)
        
        # 心动周期的阶段
        systolic_phase = 0.15 * T       # 收缩压,秒
        dicrotic_phase = 0.18 * T       # Decrotic缺口，sec
        diastolic_phase = 0.45 * T      # 舒张,sec
        
        wave = 0.0
        if t_mod < systolic_phase
            phase_ratio = t_mod / systolic_phase
            wave = 1.6 * phase_ratio * exp(1.8 * (1 - phase_ratio))
        elseif t_mod < dicrotic_phase
            phase_ratio = (t_mod - systolic_phase) / (dicrotic_phase - systolic_phase)
            wave = 1.0 - 2.2 * phase_ratio
        elseif t_mod < systolic_phase + diastolic_phase
            phase_ratio = (t_mod - dicrotic_phase) / (diastolic_phase - (dicrotic_phase - systolic_phase))
            wave = 0.3 * exp(-3.5 * phase_ratio) * sin(2*π * 1.2 * phase_ratio)
        end
        
        # 振幅调制
        spread = (systolic_bp - diastolic_bp) / 2.5
        amp_mod = exp(-0.5 * ((p - mean_bp) / spread)^2)
        
        if t_val < inflation_duration
            amp_mod = 0.0
        elseif p < 60
            amp_mod = 0.0
        elseif p > systolic_bp + 20 || p < diastolic_bp - 15
            amp_mod *= 0.02
        elseif p < 80
            amp_mod *= (p - 60) / 20
        end
        
        # 添加基本噪声
        noise = 0.03 * randn()
        oscillations[i] = 5.2 * amp_mod * (wave + noise)
        
        # 产生运动伪影
        if movement
            if t_val > inflation_duration && p <= movement_start && p >= movement_end
                if !movement_detected
                    movement_detected = true
                    movement_start_time = t_val
                end
                oscillations[i] += 3.5 * sin(2π * 0.8 * t_val)
                if rand() < 0.15
                    oscillations[i] += 8.0 * randn()
                end
            else
                if movement_detected
                    movement_detected = false
                    movement_end_time = t_val
                end
            end
        end
    end
    
    # 如果运动还没有结束
    if movement && movement_detected
        movement_end_time = time[end]
    end
    
    return oscillations, movement_start_time, movement_end_time
end;

第三阶段包括分析获得的结果并使用普遍接受的方法检测DM，DD和SrAD。

function analysis(raw_osc, deflate_phase, pressure, time, fs)
    # 信号滤波
    ttime = 0.15                        # 过滤反应时间
    n = Int(round(ttime * fs))          # 过滤窗口大小
    b = ones(n)/n                       # FIR滤波器系数
    filt_osc = filt(b, [1], raw_osc)    # 应用过滤器
    
    # 检测袖带中的压力放气阶段
    deflate_osc = filt_osc[deflate_phase]
    deflate_pres = pressure[deflate_phase]
    deflate_times = time[deflate_phase]
    
    # 建造信封
    abs_osc = abs.(deflate_osc)
    window_size = Int(round(0.45 * fs))
    window = ones(window_size)/window_size
    envelope = filt(window, [1], abs_osc)
    
    # 我们找到最大振幅
    max_amp, max_idx = findmax(envelope)
    pres_max = deflate_pres[max_idx]
    
    # 收缩压
    systolic_lim = 0.15 * max_amp
    systolic_region = findfirst(x -> x >= systolic_lim, envelope[1:max_idx])
    systolic_idx = something(systolic_region, 1)
    systolic = deflate_pres[systolic_idx]
    
    # 舒张压
    diastolic_lim = 0.55 * max_amp
    diastolic_region = findlast(x -> x >= diastolic_lim, envelope[max_idx:end])
    diastolic_idx = diastolic_region === nothing ? length(envelope) : max_idx - 1 + diastolic_region
    diastolic = deflate_pres[diastolic_idx]
    
    return deflate_osc, deflate_times, envelope, systolic, diastolic, pres_max
end;

第一种情况。

在这种情况下，我们会考虑患者遵守测量血压规则的情况。

# 第一个场景：没有运动伪影
pressure, deflate_phase, plato_end = generate_pressure(t)
time_period = time_pump + time_plato
osc, _, _ = generate_oscillations(pressure, t, pulse_rate, systolic_bp, diastolic_bp, mean_bp, time_period, movement=false)
deflate_osc, deflate_times, envelope, systolic_est, diastolic_est, pres_max = analysis(osc, deflate_phase, pressure, t, fs)
pres_osc = pressure .+ osc

# 创建具有两个子图的图形
plt_1 = plot(layout=(2,1), size=(900, 700), legend=:topright, margin=40*Plots.px)

# 上图：眼压计袖带中的压力
plot!(plt_1[1], t, pres_osc, 
    title = "眼压计袖带中的压力（无伪影）",
    xlabel = "时间，从",
    ylabel = "压力，mmHg",
    label = "袖带压力",
    linewidth = 2,
    color = :blue,
    ylims = (0, 200),
    xlims = (0, total),
    grid = true)

# 关键压力点
hline!(plt_1[1], [systolic_est], linestyle=:dash, color=:red, 
       label="收缩压(ism:$(round(systolic_est,digits=1))")
hline!(plt_1[1], [diastolic_est], linestyle=:dash, color=:purple, 
       label="舒张(ism:$(round(diastolic_est,digits=1))")
hline!(plt_1[1], [pres_max], linestyle=:dash, color=:brown, 
       label="平均值(ism:$(round(pres_max,digits=1))")

# 真正的价值
hline!(plt_1[1], [systolic_bp], linestyle=:dot, color=:red, label="收缩压(real:$systolic_bp)")
hline!(plt_1[1], [diastolic_bp], linestyle=:dot, color=:purple, label="舒张（真实：$diastolic_bp）")

# 下图：压力波动
plot!(plt_1[2], t, osc, 
    title = "压力波动",
    xlabel = "时间，从",
    ylabel = "振幅，mmHg",
    label = "振荡",
    linewidth = 1.5,
    color = :blue,
    xlims = (0, total),
    grid = true)

# 包络仅处于放气阶段
plot!(plt_1[2], deflate_times, envelope, 
    label = "振荡包络",
    linewidth = 2.0,
    color = :red)

# 定义阈值
systolic_lim = 0.15 * maximum(envelope)
diastolic_lim = 0.55 * maximum(envelope)
hline!(plt_1[2], [systolic_lim], linestyle=:dash, color=:green, label="收缩压阈值(15%)")
hline!(plt_1[2], [diastolic_lim], linestyle=:dash, color=:orange, label="舒张阈值(55%)")
hline!(plt_1[2], [maximum(envelope)], linestyle=:dash, color=:purple, label="最大振幅")

# 第一场景的结果的输出
println("第一种情况的分析结果（缺乏运动）:")
println("收缩压:$(round(systolic_est,digits=1))mmHg(real:$systolic_bp)")
println("舒张压:$(round(diastolic_est,digits=1))mmHg(real:$diastolic_bp)")
println("平均血压:$(round(pres_max,digits=1))mmHg")
println("脉搏率：△pulse_rate beats/min")

# 图形显示
display(plt_1)
savefig(plt_1, "tonometer_no_movement.png")

Результаты анализа первого сценария (отсутствие движения):
Систолическое давление: 120.4 мм рт.ст. (реальное: 120.0)
Диастолическое давление: 73.2 мм рт.ст. (реальное: 80.0)
Среднее артериальное давление: 92.2 мм рт.ст.
Пульс: 75 уд/мин

"/user/Demo_public/biomedical/modeling_blood_pressure_monitor/tonometer_no_movement.png"

从第一种情景的仿真结果可以看出，所得模型检测人体血压波动的准确度至少为5%。

第二种情况。

在这种情况下，我们会考虑患者忽略测量血压规则的情况。

# 第二个场景：带有运动伪影
osc_movement, movement_start_time, movement_end_time = generate_oscillations(pressure, t, pulse_rate, systolic_bp, diastolic_bp, mean_bp, time_period, movement=true)
deflate_osc_movement, deflate_times_movement, envelope_movement, systolic_est_movement, diastolic_est_movement, max_pres_movement = analysis(osc_movement, deflate_phase, pressure, t, fs)
pressure_movement = pressure .+ osc_movement

# 创建图形（两个子图形）
plt_2 = plot(layout=(2,1), size=(900, 700), legend=:topright, margin=40*Plots.px)

# 上图：带伪影的眼压计袖带压力
plot!(plt_2[1], t, pressure_movement, 
    title = "眼压计袖带中的压力（带有运动伪影）",
    xlabel = "时间，从",
    ylabel = "压力，mmHg",
    label = "袖带压力",
    linewidth = 2,
    color = :blue,
    ylims = (0, 200),
    xlims = (0, total),
    grid = true)

# 突出显示运动伪像区域
if !isnan(movement_start_time) && !isnan(movement_end_time)
    vspan!(plt_2[1], [movement_start_time, movement_end_time], alpha=0.25, color=:magenta, label="运动伪影")
end

# 关键压力点
hline!(plt_2[1], [systolic_est_movement], linestyle=:dash, color=:red, label="收缩压(ism:$(round(systolic_est_movement,digits=1))")
hline!(plt_2[1], [diastolic_est_movement], linestyle=:dash, color=:purple,label="舒张(ism:$(round(diastolic_est_movement,digits=1))")
hline!(plt_2[1], [max_pres_movement], linestyle=:dash, color=:brown, label="平均值(ism:$(round(max_pres_movement,digits=1))")

# 实际压力值
hline!(plt_2[1], [systolic_bp], linestyle=:dot, color=:red, label="收缩压(real:$systolic_bp)")
hline!(plt_2[1], [diastolic_bp], linestyle=:dot, color=:purple, label="舒张（真实：$diastolic_bp）")

# 底部图表：带有运动伪影的振荡
plot!(plt_2[2], t, osc_movement, 
    title = "带有伪影的压力波动",
    xlabel = "时间，从",
    ylabel = "振幅，mmHg",
    label = "振荡",
    linewidth = 1.5,
    color = :blue,
    xlims = (0, total),
    grid = true)

# 信封
plot!(plt_2[2], deflate_times_movement, envelope_movement, 
    label = "信封",
    linewidth = 2.0,
    color = :red)

# 定义阈值
systolic_lim_mov = 0.15 * maximum(envelope_movement)
diastolic_lim_mov = 0.55 * maximum(envelope_movement)
hline!(plt_2[2], [systolic_lim_mov], linestyle=:dash, color=:green, label="收缩压阈值(15%)")
hline!(plt_2[2], [diastolic_lim_mov], linestyle=:dash, color=:orange, label="舒张阈值(55%)")
hline!(plt_2[2], [maximum(envelope_movement)], linestyle=:dash, color=:purple, label="最大振幅")

# 运动伪影的区域
if !isnan(movement_start_time) && !isnan(movement_end_time)
    vspan!(plt_2[2], [movement_start_time, movement_end_time], 
           alpha=0.15, color=:magenta, label="运动伪影")
end

# 第二种情况的结果输出
println("\第二种情况的分析结果（运动的存在）:")
println("收缩压:$(round(systolic_est_movement,digits=1))mmHg(real:$systolic_bp)")
println("舒张压:$(round(diastolic_est_movement,digits=1))mmHg(real:$diastolic_bp)")
println("平均血压:$(round(max_pres_movement,digits=1))mmHg")
println("脉搏率：△pulse_rate beats/min")

# 图形显示
display(plt_2)
savefig(plt_2, "tonometer_with_movement.png")

Результаты анализа второго сценария (наличие движения):
Систолическое давление: 137.0 мм рт.ст. (реальное: 120.0)
Диастолическое давление: 81.3 мм рт.ст. (реальное: 80.0)
Среднее артериальное давление: 118.4 мм рт.ст.
Пульс: 75 уд/мин

"/user/Demo_public/biomedical/modeling_blood_pressure_monitor/tonometer_with_movement.png"

根据第二种情况的结果，很明显收缩压的确定是不正确的。由于在测量期间由人体运动引起的伪影而发生错误检测。

在不存在和存在由人类运动引起的伪影的情况下测量人类压力的情景的可视化。

# 第一和第二场景建模的结果的可视化
# 第一场景的结果的输出
println("第一种情况的分析结果（缺乏运动）:")
println("收缩压:$(round(systolic_est,digits=1))mmHg(real:$systolic_bp)")
println("舒张压:$(round(diastolic_est,digits=1))mmHg(real:$diastolic_bp)")
println("平均血压:$(round(pres_max,digits=1))mmHg")
println("脉搏率：△pulse_rate beats/min")

# 图形显示
display(plt_1)

# 第二种情况的结果输出
println("\第二种情况的分析结果（运动的存在）:")
println("收缩压:$(round(systolic_est_movement,digits=1))mmHg(real:$systolic_bp)")
println("舒张压:$(round(diastolic_est_movement,digits=1))mmHg(real:$diastolic_bp)")
println("平均血压:$(round(max_pres_movement,digits=1))mmHg")
println("脉搏率：△pulse_rate beats/min")

# 图形显示
display(plt_2)

Результаты анализа первого сценария (отсутствие движения):
Систолическое давление: 120.4 мм рт.ст. (реальное: 120.0)
Диастолическое давление: 73.2 мм рт.ст. (реальное: 80.0)
Среднее артериальное давление: 92.2 мм рт.ст.
Пульс: 75 уд/мин

Результаты анализа второго сценария (наличие движения):
Систолическое давление: 136.3 мм рт.ст. (реальное: 120.0)
Диастолическое давление: 81.1 мм рт.ст. (реальное: 80.0)
Среднее артериальное давление: 129.0 мм рт.ст.
Пульс: 75 уд/мин

结论

在这个例子中，考虑了具有振荡压力测量方法的数字眼压计模型。还考虑了两种情况:

*无人工（病人）运动的袖带压力;
*带人体（患者）运动伪影的袖带压力。

分析所获得的结果，可以说，为了正确确定一个人的收缩压和舒张压，必须遵循一些规则。:

测量应在舒适，平静的环境中进行，房间应在室温下;
只有在至少五分钟的休息期后才能测量血压。;
应该记住，肩膀不应该被衣服挤压，特别是因为通过衣服测量血压是不正确的。;
在测量过程中不要移动或说话。

在初始测量期间，应在双手上确定血压，然后在压力较高的手上测量血压。（手上的血压差异高达10-15mmHg是正常的。)