脑电数据预处理十六：自动化阈值检测原理与实践
在脑电图（EEG）数据分析中，伪迹是影响信号质量和分析结果准确性的主要因素。传统的伪迹处理方法通常依赖于人工检查和标记，耗时且主观性强。为了解决这一问题，**自动化阈值检测（Automated Thresholding）**应运而生，它是一种简单、高效且可重复的伪迹识别技术。
自动化阈值检测的核心思想是，基于信号幅度的统计学异常值来识别和标记伪迹。这些伪迹通常具有比正常脑电信号高得多的振幅，因此可以通过设定一个合理的阈值来将其筛选出来。一旦伪迹被标记，它们可以被直接剔除（Epochs Rejection）或被其他方法（如插值）所取代。
本文将深入探讨自动化阈值检测的原理、常用方法，并结合 MNE-Python 库提供完整的代码实现，旨在帮助你掌握这一基础而实用的预处理技术。
自动化阈值检测的核心原理
自动化阈值检测基于一个简单的统计学假设：伪迹是异常值。
信号振幅：正常脑电信号的振幅通常在微伏（µV）级别，而伪迹（如眨眼、肌肉活动）的振幅可以达到数百甚至数千微伏。
阈值设定：通过对信号振幅或其变化率设定一个上限值，任何超过这个值的样本或时间窗都会被标记为伪迹。
统计学方法：阈值的设定可以是固定的经验值（例如，pm100muV），也可以是基于数据自身统计特征的自适应值（例如，信号标准差的3倍）。
常用方法与应用场景
1. 幅度阈值（Amplitude Thresholding）
原理：直接检查每个时间点的信号幅度，如果超过预设的最大/最小值，则标记为伪迹。
应用：最简单直接的方法，适用于处理振幅剧烈的伪迹，如电极接触不良、剧烈眨眼。
缺点：对于幅度相对较小但频率异常的伪迹（如高频肌肉活动），效果不佳。
2. 峰度阈值（Kurtosis Thresholding）
原理：**峰度（Kurtosis）**是描述数据分布“尖峰”程度的统计量。高振幅、短暂的伪迹会导致数据分布的峰度显著增加。通过检测峰度异常的时间窗，可以识别出伪迹。
应用：比幅度阈值更鲁棒，可以更有效地识别出突然出现的伪迹，尤其是在数据整体振幅变化不大的情况下。
3. 变化率阈值（Rate-of-Change Thresholding）
原理：伪迹通常表现为信号的急剧跳变，这意味着其**变化率（一阶导数）**会非常大。通过对信号的变化率设定阈值，可以识别出这些尖锐的伪迹。
应用：适用于识别尖峰伪迹，如电极突然移动或信号饱和。
Python实现：基于MNE-Python
MNE-Python 提供了多种阈值检测和数据剔除的方法，通常集成在 mne.Epochs 对象的创建过程中，或用于识别坏通道。
1. 环境准备与数据加载
首先，确保你已经安装了所需的库。
Bash
pip install mne
Python
import mne
import numpy as np
import matplotlib.pyplot as plt
from mne.datasets import sample# 1. 加载MNE示例数据
data_path = sample.data_path()
raw_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_raw.fif'
event_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_raw-eve.fif'raw = mne.io.read_raw_fif(raw_fname, preload=True)
events = mne.read_events(event_fname)# 仅选择EEG通道，并进行初步滤波
raw.pick_types(meg=False, eeg=True)
raw.set_eeg_reference('average', projection=True)
raw.filter(1., 40., fir_design='firwin')
2. 应用阈值检测进行分段剔除
我们将定义一个分段（Epochs），并在分段过程中使用 reject 参数进行自动化阈值剔除。
Python
# 2. 定义事件和分段参数
event_id = {'auditory/left': 1, 'visual/left': 3}
tmin, tmax = -0.2, 0.5 # 分段时间窗# 3. 设置自动化阈值剔除
# 任何一个EEG通道的幅度超过 ±100 µV，则整个分段被剔除
reject_criteria = dict(eeg=100e-6) # 单位为伏特 (V)# 4. 创建epochs并应用阈值
# MNE将自动识别并剔除不符合标准的分段
epochs_with_rejection = mne.Epochs(raw, events, event_id, tmin, tmax,preload=True,reject=reject_criteria
)print(f"原始分段数: {len(epochs_with_rejection.drop_log)}")
print(f"被剔除的分段数: {len(epochs_with_rejection.drop_log) - len(epochs_with_rejection)}")
3. 可视化剔除效果
通过绘制分段数据的头皮图，我们可以看到那些因为伪迹而被剔除的试验。
Python
# 可视化一个被剔除的分段，以了解原因
# 查找第一个被剔除的分段索引
dropped_indices = np.where([len(d) > 0 for d in epochs_with_rejection.drop_log])[0]if len(dropped_indices) > 0:first_dropped_idx = dropped_indices[0]# 重新创建原始分段（不剔除）以进行对比epochs_raw = mne.Epochs(raw, events, event_id, tmin, tmax, preload=True)# 绘制被剔除分段的通道波形图epochs_raw[first_dropped_idx].plot(picks='eeg')plt.title(f"被剔除的分段 (索引: {first_dropped_idx}) - 伪迹幅度超标")plt.show()# 绘制平均ERP波形图对比
# 如果有足够多的分段被保留，可以绘制平均波形
if len(epochs_with_rejection) > 10:evoked_raw = epochs_raw.average(picks='eeg')evoked_cleaned = epochs_with_rejection.average(picks='eeg')fig, ax = plt.subplots(1, 1, figsize=(8, 6))evoked_raw.plot(spatial_colors=True, axes=ax, show=False)evoked_cleaned.plot(spatial_colors=True, axes=ax, show=False)ax.set_title("平均ERP波形对比 (红色为剔除后)")ax.legend(['原始', '剔除后'])plt.show()
总结与应用场景
自动化阈值检测是一种简单、快速且高效的脑电伪迹处理方法。尽管它不如ICA或ASR等高级算法那样灵活，但它在处理振幅剧烈、突发的伪迹时非常有效。
优点：
快速：计算量小，适用于大规模数据处理。
可重复：阈值一旦设定，处理结果是完全确定的。
简单直观：易于理解和实现。
应用场景：
初步数据清理：在任何复杂处理之前，使用阈值剔除可以快速移除最明显的伪迹。
分段剔除：在ERP研究中，这是一种标准方法，用于确保平均结果不受个别伪迹分段的影响。
BCI：在某些脑机接口应用中，快速剔除高振幅伪迹可以提高信号质量。
将自动化阈值检测与其他预处理方法（如滤波、ICA）结合使用，可以构建一个鲁棒的脑电数据处理流程。