脑电信号爆发抑制比计算方式学习笔记

提示：本贴只做个人学习记录用

---

前言

脑电信号的爆发抑制比(Burst Suppression Ratio, BSR)是麻醉深度监测中的重要指标，用于量化脑电图中爆发活动与抑制活动的比例。

本文内容主要来源于对论文《Real-time segmentation of burst suppression patterns in critical care EEG monitoring》的阅读和理解，相当于本论文的阅读笔记。

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、爆发抑制比是什么？

爆发抑制特征：

爆发(Burst)：高振幅、相对高频的脑电活动

抑制(Suppression)：低振幅、近乎平坦的脑电活动段

BSR = 抑制时间 / 总时间 × 100%

爆发抑制模式在新生儿发育过程中以及各种新生儿脑损伤以及所有年龄段的弥漫性缺氧性脑损伤患者中都会自发出现。爆发抑制模式也可能由麻醉药物的使用诱导，例如在手术麻醉期间或用于治疗难治性癫痫持续状态，或者通过控制性低温疗法，这在心脏手术中的完全循环停止手术中很常见。

二、计算公式及代码

1.论文研究方式 & 计算公式：

实际采集临床脑电信号，并医生手动标记出抑制段和爆发段，最后用代码处理后的结果和医生手动标记的结果Y1 Y2做对比。

文章中使用的是递归方差分析（阈值+局部信号方差），具体的公式如下图：

上述公式中，μ是当前信号段的均值，σ是当前信号段的方差，x是信号数据，其中t、t-1表示不同的时间点的数据点索引，当前的计算结果与上一刻的数据有关。公式（1）和（2）就表示了自适应均值与自适应方差的计算过程。

其中δ表示某函数关系，即方差小于阈值θ的时候，就幅值为1，表示抑制段；否则（≥的时候），就赋值为0，表示爆发段。

上述公式中最关键的两个参数，遗忘因子 β 和分类阈值 θ。下面单独介绍两个参数的计算思路和流程。

2.遗忘因子 β

遗忘因子有两个特点：
（1）介于0~1之间；
（2）与遗忘时间 T 成反比：
其中Fs是信号采样率。

根据论文思路，进行了如下步骤操作：
（1）临床依据表明合适的遗忘时间在 0.01s~3s 之间；
（2）在（1）中的范围内进行逐点计算（精细化网格搜索），计算出各个 T 值对应的 β值，并参与 1 中的递归方差分析方程；
（3）求出各个 β值对应的最优阈值 θ；
（4）看最终结果，分类误差最小的 T（104.7ms），就是最优遗忘因子β（0.9534）。

3.最优阈值 θ & 分类误差 e

根据论文思路，进行了如下步骤操作：
（1）求出临床实际数据的每段的方差，按照方差降序排列；
（2）每个方差都视作一个阈值，进行 δ 函数关系的分类；
（3）每个 θ值都对应计算如下参数：
a. 识别为 Burst 的点数 Tb；
b. 识别为 suppression 的点数 Ts；
c. 实际的 Burst 的点数 Bb；
d. 实际的 suppression 的点数 Bs。
（4）当前遗忘因子 β值 和阈值 θ值对应的分类误差为：

找到误差 e值最小的阈值 θ （一般临床常用固定的5uV做阈值）和遗忘因子 β，即是最优结果。

论文中还有更进一步的优化结果，但是本贴中不进行细说，性能优化不大，下面上matlab脚本代码实例。

4.matlab脚本

%% EEG爆发-抑制模式检测脚本 (Matlab版本)
% 本脚本演示了如何在ICU脑电图(EEG)通道上检测爆发-抑制模式(burst-suppression patterns)
% 使用了一种基于局部方差阈值的简单方法
% 该方法在以下论文中被描述：
% 'Real-time segmentation of burst suppression patterns in critical care EEG monitoring',
% by Westover et al.
clear;
% close all;
fclose('all');
clc;


%% 参数设置
% 使用模拟数据进行演示（这样不需要实际的EDF文件）
use_simulation = 0;

if use_simulation
    % 生成模拟的EEG信号（包含爆发-抑制模式）
    fprintf('使用模拟数据进行演示...\n');
    samplingrate = 250;  % 采样率（Hz）
    duration = 60;       % 模拟信号时长（秒）
    
    % 生成模拟信号
    sigbuf = generate_simulation_eeg(samplingrate, duration);
    
    % 要绘制的时间间隔（单位：秒）
    plot_interval = [5, 25];  % 绘制5-25秒的信号
else
    % 读取真实脑电信号
    load('...\Annotations_1.mat');
    load('...\p1_data.mat')
    samplingrate = Fs;  % 采样率（Hz）
    sigbuf = data';
    plot_interval = [5, 25];  % 绘制5-25秒的信号

end

% 抑制阈值参数theta
theta = 50;  % 按照5uV设置即可，考虑ADC增益，视情况而定
min_suppr_len = samplingrate/2;     % 最小抑制段长度-0.5s

%% 使用递归方差估计对信号进行分段处理
[res, sig_segmented, sig_recvar] = recursive_variance_estimation(sigbuf(:,1), 0.9633, theta, 50);

sigbuf2 = (sigbuf - mean(sigbuf)) / (5 * rms(sigbuf));      % 信号归一化，画图方便

%%% Y1 Y2都是0是爆发段，1是抑制段，要处理一下
Y1AFT = ~Y1;
Y2AFT = ~Y2;

figure,
plot(sig_segmented,'r'); hold on
plot