ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

YASA自动化多导睡眠图分析:从EDF到睡眠分期与纺锤波检测

YASA自动化多导睡眠图分析:从EDF到睡眠分期与纺锤波检测 简介这是一份面向睡眠研究人员、脑电数据分析者及Python开发者的YASA工具箱完整源码包。YASA专注于多导睡眠图PSG的自动分析涵盖自动睡眠分期、纺锤波/慢波/快速眼动事件检测、伪影剔除、频谱分析及催眠图统计等功能需要结合NumPy、Pandas、MNE与Jupyter Lab使用。资源共254个文件压缩包约115.78MB主体为Python脚本27个py和Jupyter Notebook示例16个ipynb同时包含HTML文档、CSS/JS网页资源、PNG/SVG图示、RST说明及少量测试数据npz、fif、joblib等目录结构完整便于直接参考与二次开发。目前已有803人学习下载适合具备一定Python基础、希望快速搭建睡眠分期或事件检测流程的研究者。通过源码、示例文档与可视化图表读者可掌握YASA的调用方式、参数配置思路及典型分析流程节省从零摸索环境与接口的时间。1. YASA 是来给多导睡眠图做自动分析的不是又一个人工智能框架睡眠实验室最耗人的一件事是夜里导出的七八个小时多导睡眠图PSG要一段一段人工判读。脑电、眼电、肌电十几个通道叠在一起先分睡眠分期再回头找纺锤波、慢波和快速眼动光靠肉眼拖波形一晚记录的处理时间按小时算。YASAYet Another Spindle Algorithm是一个用 Python 写的 PSG 分析包名字听上去像又造了一个纺锤波检测轮子实际它把自动睡眠分期、睡眠统计指标、纺锤波/慢波/REM 事件检测和频谱分析全收进了同一套接口。它能直接吃 EDF 转出来的 DataFrame输出的是带时间戳的事件表和完整睡眠参数而不是一张没法复核的图。适合刚接触脑电数据想快速出指标的工程师也适合已经在用 MNE 做预处理、想把睡眠事件指标接进现有管道的分析人员。有一点要提前说清楚YASA 计算的是可复现的量化指标不等同于临床诊断结论用它做科研和筛选可以做临床报告要有人工复核环节。2. 睡眠记录的数据准备EDF、通道命名与 hypno 时间轴2.1 多导睡眠图里的三类信号与 30 秒分期规则PSG 记录里通常包含脑电EEG常见 C3-M2、C4-M1、Fz-Cz、Pz-Oz、眼电EOG左右眼各一个通道、下颌肌电EMG另外还有心电、口鼻气流、血氧饱和度这些辅助通道。YASA 真正用到的是前三类辅助通道可以留着但不会参与睡眠分期和事件检测。国际通用的睡眠分期把一夜切成一个个 30 秒的 epoch每个 epoch 打一个标签清醒W、N1、N2、N3 或 REM。整套分析的时间轴都挂在这个 30 秒网格上所以数据准备的核心不是训练什么模型而是把原始信号和分期标签在时间上对齐。很多第一次跑 YASA 的人坑都不在算法上而是 hypno 的时间索引起点跟数据起点错了几秒最后算出来的睡眠潜伏期和睡眠效率整体漂移。2.2 用 MNE 读 EDF、统一单位与采样率YASA 的输入约定很明确数据是一个 pandas DataFrameindex 为采样时间戳columns 为通道名hypno 是一个 Seriesindex 为每个 epoch 的起始时间戳值为 0-4 或 W/N1/N2/N3/R。从 EDF 文件到这种结构最常见的路径是先用 MNE 读入再转成 DataFrameimport mne import pandas as pd raw mne.io.read_raw_edf(night1.edf, preloadTrue) # PSG 里 EEG 和 EOG 采样率不一致时先统一 if raw.info[sfreq] ! 500: raw.resample(500) # 过滤掉直流漂移和肌电高频段保留睡眠分析的频带 raw.filter(0.3, 35, picks[eeg, eog]) raw.notch_filter(50, picks[eeg, eog]) # 市电干扰按本地频率选 50 或 60 df raw.to_data_frame() # MNE 默认单位是 VYASA 期望 uV必须换算 for col in df.columns: if any(k in col.lower() for k in (eeg, eog, emg)): df[col] * 1e6这段代码有三个地方不能省。第一是重采样到统一采样率YASA 内部对事件检测采用固定时间窗采样率不一致会让事件起止时间在不同通道间对不齐。第二是 50 Hz 陷波纺锤波频带是 12–15 Hz市电干扰不会直接踩在频带里但它的谐波会污染包络检测出来的事件会偏多。第三是单位换算YASA 内部所有幅度阈值都基于微伏设计EDF 文件里多数采集系统导出的是 V不换算的话幅度阈值完全失去意义慢波和纺锤波事件数量会成倍膨胀或归零。2.3 构建 hypno时间索引与标签语义hypno 的来源一般是采集软件导出的分期文本或者由睡眠技师在配套工具里逐 epoch 标注后导出。拿到的是每行一个标签的文本文件要转成带时间索引的 Seriesimport pandas as pd stages pd.read_csv(stages.txt, headerNone, names[stage]) hypno pd.Series(stages[stage].values) # 索引必须跟数据起点对齐第一个 epoch 从数据第一秒开始 epoch_start pd.date_range(df.index[0], periodslen(hypno), freq30s) hypno.index epoch_start hypno.name Hypnogram标签语义上要注意YASA 接受的 hypno 值是 0W、1N1、2N2、3N3、4REM也接受直接的文本标签。不同设备导出的分期文件可能用 AASM 标准的 R 表示 REM用 W 表示清醒处理文本时先做一次映射表统一成 0-4能省掉后面很多判断逻辑。一个常见的隐蔽错误是 epoch 起始时间戳写成了 epoch 中点或者第一段记录前有几分钟的静息态没算 epoch这两种情况会让事件检测阶段把事件错误挂到相邻分期上。如果记录里完全没有人工标注又想快速摸底可以用 YASA 自带的自动分期生成一版 hypno 先跑通流程但正式的睡眠报告参数不能直接采信这个边界后面会展开。2.4 通道命名的约定与检查YASA 判断通道类别靠的是列名里的关键字EEG 通道名要带 eeg 或直接给出脑电导联名并辅以 eeg_name 参数EOG 要带 eogEMG 要带 emg。如果从 EDF 读进来通道名是 Cz、LOC 这种需要手工补上类别标记。我通常在拿到 DataFrame 后先打印通道列表确认每个通道能被正确归类再进入检测流程否则 spindles_detect 里找不到 EOG 或 EMG 通道时报错信息会很隐晦。3. 睡眠分期与睡眠统计参数的计算逻辑和边界3.1 sleep_stats 返回的核心睡眠指标数据准备好之后第一步通常是算整夜睡眠统计。YASA 提供了 sleep_stats输入数据、hypno 就能得到一串常用临床指标import yasa stats yasa.sleep_stats(dat, hypno) stats输出是一个表格其中比较关键的字段列在下面字段含义计算口径Total time in bed卧床总时长hypno 第一个到最后一个 epoch 的时间跨度Sleep period time睡眠周期时长入睡到最终醒来的时长Total sleep time总睡眠时间所有非清醒 epoch 的累加时长Sleep onset latency入睡潜伏期从记录开始到第一个非清醒 epoch 的时间WASO入睡后清醒时间入睡之后所有清醒 epoch 的累加Number of awakenings觉醒次数入睡后转为清醒的连续片段数Sleep efficiency睡眠效率TST 除以 TBT 的百分比REM latencyREM 潜伏期入睡到第一个 REM epoch 的时间Stage durations各期时长与占比按 W/N1/N2/N3/REM 分别统计这些字段在 sleep_stats 的结果里直接可以取到但口径要按研究需求确认一遍。比如 Sleep onset latency 的定义有些 protocol 要求必须连续出现三个非清醒 epoch 才算真正入睡YASA 默认取第一个非清醒 epoch两边会差几分钟到十几分钟。如果你所在的项目采用的是连续三 epoch 标准需要在拿到结果后自己再算一版不能直接改参数让 YASA 换口径。3.2 bandpower 是分期和事件检测的共同底座睡眠分期本质上是在看不同频带功率随睡眠周期的变化。YASA 里最直接的工具是 bandpower它对每个通道每个 epoch 做 Welch 功率谱估计再按给定频带积分bands [(0.5, 4, Delta), (4, 8, Theta), (8, 12, Alpha), (12, 15, Sigma), (15, 30, Beta), (30, 45, Gamma)] bp yasa.bandpower(dat, hypnohypno, bandsbands) print(bp.head())返回结果里每行是一个通道列包含每个频带的绝对功率、总绝对功率、相对功率和相对总功率。做纺锤波分析时Sigma 频带的相对功率是最该先看的指标如果受试者整夜 Sigma 相对功率明显偏低说明这一夜可能根本没有多少典型纺锤波或记录中 EMG 噪声把信号盖掉了。这种情况下直接跑检测器得到的零星事件很难用于统计。我一般先用 bandpower 画一遍各通道的功率谱变化确认要检测的频带确实有信号再进事件检测这一步能省掉后面大量无意义的调参。3.3 自动分期的模型边界与工程兜底YASA 的自动睡眠分期用的不是深度学习而是人工设计的时域加频域特征加随机森林分类器。特征包括各频带相对功率、频谱边缘频率、频谱熵、EMG 功率等分类器按 30 秒 epoch 逐段预测最后输出整夜 hypno。好处是计算快、依赖少、特征可解释坏处是训练数据以正常成人为主模型对儿童、老年人、睡眠呼吸暂停患者和部分神经系统疾病记录的泛化能力有限。实际项目中这几类记录经常出现 N1 被低估、N3 被高估的情况因为病理脑电的背景频率和慢波形态与训练集差异较大。把自动分期当成预标注用是完全合理的直接拿它的输出做最终睡眠报告风险就高了。我的做法是自动分期结果永远保留概率或置信度字段只把置信度高的段落入自动统计边界段交给人工复核。4. 纺锤波检测的参数、算法与慢波/REM 事件输出4.1 纺锤波检测带通、滑动阈值与事件合并YASA 名称里的 Spindle Algorithm 是整个包最核心的部分。纺锤波是睡眠 N2/N3 期出现的 12–15 Hz 震荡典型时长在 0.3 到 3 秒幅度在脑电背景上并不总是明显。自动检测的经典做法不是神经网络而是信号处理三步带通滤波、包络阈值、事件合并。首先是 12–15 Hz 带通把信号限制在 Sigma 频带然后对包络计算滑动平均设置高、低两档相对阈值低阈值确定事件边界高阈值确认图中峰值最后把被间隙分开但实际连续的事件合并。YASA 内部用 numba 做加速单通道一整夜数据跑下来通常只需要几秒到几十秒比纯 Python 实现快得多。调用入口很简洁sp yasa.spindles_detect( dat, hypnohypno, include(2, 3), # 只在 N2 和 N3 阶段检测 freq_sp(12, 15), # Sigma 频带 min_duration0.3, # 最短事件时长秒 remove_edgesFalse, # 是否剔除记录首尾的事件 verboseTrue ) print(sp)sp 是一个 pandas DataFrameindex 是通道名每一行代表检测到一个纺锤波事件包含 Start、Peak、End 时间戳Duration、Frequency、RMS、AbsPower、RelPower 这些特征列如果传了 hypno还会带上 Stage 列标记这个事件落在哪个睡眠分期。这个输出结构是之后做统计和可视化的基础建议在项目开始阶段就把事件表规范成统一的 parquet 格式别每跑一次就换一种导出。4.2 关键参数表与调整场景spindles_detect 的默认参数对绝大多数成人整夜记录是合理的但项目里真正花时间的地方是理解每个参数在什么场景下要动。参数默认值需要调整的场景freq_sp(12, 15)儿童纺锤波频率偏高可放宽到 (12, 16)include(2, 3)只想统计 N2 时改为 (2,)min_duration0.3信号噪声大时调高到 0.5 过滤假阳remove_edgesFalse分析完整夜时建议保持 Falsethresholds相对功率两档阈值低幅纺锤波为主时下调高阈值档同时接受更多假阳阈值参数是误用高发区。YASA 的阈值是相对功率而不是绝对微伏值理解为相对背景功率的高出倍数更准确因此不同受试者之间不需要因为放大器增益不同而反复调整。真正需要调的场景是纺锤波幅值整体偏低的记录比如额叶导联检出的纺锤波比中央导联低这时把阈值下调会找回更多事件但要同步用可视化确认找回的是真实纺锤波而不是 EMG 碎片。另一个高频错误是在有周期性肢体运动的通道上硬跑检测运动伪迹会在 Sigma 频带产生类似事件跑之前先看原始波形的伪迹情况。4.3 慢波与 REM 检测的入口和输出列慢波检测的入口是 sw_detect原理是对负向峰值幅度做自适应阈值匹配。与纺锤波不同慢波的幅度绝对值更大检测器会先估计整夜负峰幅度的分布再基于分布确定阈值所以它不像纺锤波那样对相对功率阈值敏感sw yasa.sw_detect( dat, hypnohypno, include(3,), # 慢波主要出现在 N3 freq_sw(0.5, 4.5) ) print(sw)sw 输出的特征列包括负峰出现时间和幅度、事件时长、频率、相位斜率等。REM 检测由 rem_detect 负责它依赖 EOG 通道识别快速眼动同时用 EMG 通道确认肌电抑制状态以排除单纯清醒时的眼动。调用时要确保 EOG 通道带 eog 标记或者直接用 eog_name 参数指定rem yasa.rem_detect(dat, hypnohypno, include(4,), eog_name(LOC, ROC))REM 检测的实际产出不只是一次次眼动事件REM 密度每分钟眼动次数是一个常用指标后面做组间对比时可以直接由事件表算出来。4.4 按睡眠分期统计纺锤波密度检测完的事件表带 Stage 列统计密度时就不需要再回头对着 hypno 逐个匹配了# 统计 N2 和 N3 各自的总分钟数 n2_n3_epochs hypno.isin([2, 3]).sum() # epoch 数 n2_n3_min n2_n3_epochs * 0.5 # 30 秒一个 epoch sp_density sp.groupby(Stage).size() / n2_n3_min print(sp_density)这里有一个容易算错的地方hypno 是 30 秒一个 epoch统计时长时要把 epoch 数乘以 0.5 换算成分钟直接用 epoch 数当分钟数会把密度夸大 2 倍。事件表和 hypno 的 Stage 列偶见不一致例如检测时 include 限定了 (2, 3)但 hypno 本身有分期待修正的片段这时事件表里的 Stage 会显示实际归属按 groupby 统计即可不要再用 hypno 重新映射。5. 批量处理多晚睡眠记录缓存、并行与结果合并5.1 单夜处理函数中的缓存设计单晚记录的处理流程可以收敛成独立的函数再缓存中间结果。事件检测是重计算睡眠统计和频谱分析相对轻量按文件粒度做缓存是最合适的选择。我通常的做法是每个受试者一个目录中间结果只存 parquet 事件表和统计表原始 EDF 和人工分期文件不进入产物目录from pathlib import Path import pandas as pd import mne import yasa def analyze_one(edf_path, stages_path, out_dir): raw mne.io.read_raw_edf(edf_path, preloadTrue) if raw.info[sfreq] ! 500: raw.resample(500) raw.filter(0.3, 35, picks[eeg, eog]) dat raw.to_data_frame() * 1e6 stages pd.read_csv(stages_path, headerNone)[0] hypno pd.Series(stages.values) hypno.index pd.date_range(dat.index[0], periodslen(hypno), freq30s) hypno.name Hypnogram stats yasa.sleep_stats(dat, hypno) stats.to_json(out_dir / f{Path(edf_path).stem}_stats.json) sp yasa.spindles_detect(dat, hypnohypno) sp.to_parquet(out_dir / f{Path(edf_path).stem}_spindles.parquet)缓存设计上parquet 比 CSV 更适合事件表因为它保留时间戳和多数值列的类型重新读入后不需要再处理字符串转时间的问题。如果一个受试者跑完又改了参数只重跑事件检测这一层统计表不用动这正是把函数拆成单夜独立流程的好处。5.2 多进程跑批与 numba 的上下文注意事项YASA 的事件检测在底层依赖 numba 的 JIT 编译批量并行时对进程启动方式有要求。写多进程批处理时直接在命令行跑没有保护会看到 numba 对象无法 pickled 的报错from concurrent.futures import ProcessPoolExecutor def run_batch(): jobs [(edf_path, stages_path, out_dir) for edf_path in edf_list] with ProcessPoolExecutor(max_workers4) as ex: list(ex.map(lambda x: analyze_one(*x), jobs)) if __name__ __main__: run_batch()并行进程数不能只看 CPU 核数。YASA 检测单通道事件时会用 numba 多线程一个 analyze_one 进程可能吃掉两到四个逻辑核工作站核再多max_workers 设成 8 也会把内存和算力打满。我一般按机器物理核数的四分之一到二分之一起步。整夜记录在 500 Hz 采样率下16 通道数据量约一两百 MB加载进内存后叠加 numba 中间变量单进程内存占用到 1 GB 很常见内存比核数更容易成为瓶颈。5.3 结果合并与后续分析批量跑完后把所有 parquet 读进来拼成大表做组间分析是标准操作。事件表的 index 是通道名合并时要把受试者 ID 作为一列写进每个事件行否则几十个文件拼在一起后没法区分来源files sorted(out_dir.glob(*_spindles.parquet)) all_sp [] for f in files: tmp pd.read_parquet(f) tmp[subject_id] f.stem.split(_)[0] all_sp.append(tmp) events pd.concat(all_sp, ignore_indexTrue)到这里整夜的统计表、事件表、频谱表已经可以进入统计建模环节。事件表里每个事件的频带功率、幅度、时长作为特征可以按受试者聚合也可以按通道聚合具体口径取决于研究问题。6. 验证检测结果的三个具体操作事件分布、可视化与误用排查检测器跑完不等于工作结束先把三个验证动作固定到流程里能拦截大部分异常结果。第一个动作是检查事件时长分布。纺锤波正常时长分布在 0.3 到 3 秒之间睡眠波检测器给出的慢波一般更长。跑完检测后第一件事是打印描述统计d sp[Duration] print(d.describe()) print((d 0.3).mean())如果事件数量异常多时长小于 0.3 秒的事件占比超过两成基本可以判定阈值过低检测到了 EMG 碎片或滤波振铃。这种情况不要急着调阈值先看原始波形确认噪声类型。第二个动作是可视化抽查。把检测到的事件叠加到原始脑电和滤波包络上肉眼确认事件峰值是否落在 Sigma 包络的强震荡段。YASA 自带 plot_events传入原始数据、hypno 和事件表就能快速抽查抽查的 epoch 不要只选第一段要覆盖前半夜、后半夜和各个分期。第三个动作是排查三类高频误用。单位错误最常见EDF 读出 V 直接进 YASA所有幅度阈值失效事件数量要么爆炸要么归零排查方式是打印数据列的最大值和标准差正常 EEG 的幅值应该在几十到几百微伏量级。hypno 与数据错位是第二个高频问题事件表里的 Stage 列和实际波形对不上时先检查 hypno 的索引起点和 freq 参数别先怀疑检测器。第三是先降采样再检测有人为了省内存把数据降到 100 Hz纺锤波频带上限 15 Hz 在 100 Hz 下虽然名义上满足奈奎斯特条件但包络细节严重失真事件起止时间会整体模糊。保留 250 Hz 以上再跑事件检测是稳妥的工程底线。把这三个检查写进你的分析管线开头比事后对着一堆异常事件排查要省时间。本文还有配套的精品资源点击获取
返回列表