一、什么是源定位1.1 通俗理解到目前为止我们分析的是头皮上的脑电 头皮电极 → 记录头皮表面的电位 ↓ 这就像在地球表面测量地震波 我们想知道的是震源在哪里多深 源定位 从头皮电位反推大脑内部的激活位置 头皮电位已知 → 逆推 → 大脑激活源未知 ↑ ↑ 测量到的 我们想知道的1.2 正问题和逆问题正问题Forward Problem 已知大脑中某处有激活 计算头皮上各电极会测到什么电位 方向大脑 → 头皮 难度简单唯一解 逆问题Inverse Problem 已知头皮上各电极测到的电位 计算大脑中哪些地方有激活 方向头皮 → 大脑 难度困难无唯一解需要额外假设1.3 形象比喻正问题 你知道手电筒的位置和方向 计算墙上哪个位置会被照亮 → 简单有唯一答案 逆问题 你看到墙上有一个光斑 推测手电筒在哪里 → 困难可能有多个位置都能产生同样的光斑二、源定位需要的原料2.1 三个关键模型┌─────────────────────────────────────────┐ │ 源定位 源空间 头模型 传感器位置 │ ├─────────────────────────────────────────┤ │ │ │ 1. 源空间Source Space │ │ 大脑皮层上可能的激活位置 │ │ 几千个候选点 │ │ 像在大脑表面画满小格子 │ │ │ │ 2. 头模型Head Model / BEM │ │ 头部的导电模型 │ │ 描述电流如何穿过大脑、颅骨、头皮 │ │ 颅骨导电差 → 信号被模糊化 │ │ │ │ 3. 传感器位置已通过 Montage 设置 │ │ 头皮上电极的 3D 坐标 │ │ │ │ 正向解 头模型 源空间 传感器 │ │ 逆解 正向解 实际数据 │ └─────────────────────────────────────────┘2.1 三个关键模型标准源定位需要被试的 MRI 扫描 MRI → 提取大脑皮层表面 → 在上面撒候选点 → 源空间 没有 MRI 怎么办 使用标准模板fsaverage 一个平均大脑的模型 MNE 的 sample 数据集自带 MRI三、环境准备与依赖安装3.1 源定位需要的额外库源定位比前面的分析需要更多依赖库一次性安装所有依赖pip install nibabel h5io pyvistaqt -i https://pypi.tuna.tsinghua.edu.cn/simple3.2 导入库和设置# 环境设置 # 关键3D 可视化由 pyvistaqt 处理基于 Qt # matplotlib 用 Agg 后端不显示窗口只用于保存图片 # Agg Anti-Grain Geometry一个非交互式后端 # 为什么用 Agg 而不是 TkAgg # 因为 stc.plot() 使用 Qt 引擎做 3D 渲染 # TkAgg 和 Qt 不能同时运行会冲突 import matplotlib matplotlib.use(Agg) # pyplotmatplotlib 的绘图接口 # 虽然用 Agg 后端但仍可保存图片 import matplotlib.pyplot as plt # mne脑电分析核心库 import mne # numpy科学计算库 import numpy as np # os操作系统路径处理 import os # warnings警告控制 import warnings warnings.filterwarnings(ignore) # 中文字体设置 plt.rcParams[font.sans-serif] [Microsoft YaHei, SimHei] plt.rcParams[axes.unicode_minus] False print(*60) print(MNE-Python 第9天源定位基础) print(*60)四、加载数据并预处理4.1 获取数据路径# ---------- 1. 获取数据路径 ---------- # mne.datasets.sample.data_path() # 返回 sample 数据集的本地路径 # sample 数据集包含被试的 MRI 扫描数据 sample_data_folder mne.datasets.sample.data_path() # subjects_dir存放被试 MRI 数据的目录 # 里面按被试名字分文件夹 subjects_dir os.path.join(sample_data_folder, subjects) print(fMRI 数据路径: {subjects_dir}) # 被试名称 # 这里用 sample实际研究中通常是编号如 sub-01 subject sample # 原始脑电数据路径 raw_fname os.path.join( sample_data_folder, MEG, sample, sample_audvis_raw.fif )4.2 加载数据# ---------- 2. 加载数据并预处理 ---------- print(\n加载数据...) # 加载原始数据 raw mne.io.read_raw_fif(raw_fname, preloadFalse) # 提取通道EEG EOG STIM # eegTrue脑电通道分析目标 # eogTrue眼电通道可用于 ICA # stimTrue刺激通道提取事件需要 raw_eeg raw.copy().pick_types(eegTrue, eogTrue, stimTrue)4.3 重命名通道并设置 Montage# 找出所有以 EEG 开头的通道名 # 列表推导式[表达式 for 变量 in 列表 if 条件] eeg_names [ch for ch in raw_eeg.ch_names if ch.startswith(EEG)] # 创建标准 10-20 蒙太奇 montage mne.channels.make_standard_montage(standard_1020) # 取相同数量的标准电极名 standard_names montage.ch_names[:len(eeg_names)] # dict(zip(A, B)) # zip 将两个列表配对 → [(EEG 001,Fz), ...] # dict 转为字典 → {EEG 001:Fz, ...} raw_eeg.rename_channels(dict(zip(eeg_names, standard_names))) # set_montage()设置电极在头皮上的 3D 坐标 # 源定位需要知道电极的精确位置 raw_eeg.set_montage(montage)4.4 滤波和重参考# load_data()将数据加载到内存 # 滤波等操作必须在加载后进行 raw_eeg.load_data() # 陷波滤波去除 60Hz 工频干扰 raw_eeg.notch_filter(freqs60, pickseeg, verboseFalse) # 带通滤波保留 1-40 Hz # l_freq1高通去除慢速漂移 # h_freq40低通去除高频肌电噪声 raw_eeg.filter(l_freq1, h_freq40, pickseeg, verboseFalse) # 关键用投影方式做平均参考 # projectionTrue使用投影算子而非直接修改数据 # 为什么必须用投影方式 # 源定位需要在计算过程中动态处理参考问题 # 直接修改数据会破坏这种灵活性 raw_eeg.set_eeg_reference(average, projectionTrue) print(✅ 数据预处理完成)projectionTrue的含义普通重参考 (projectionFalse) 直接修改数据值 Ch1_new Ch1 - average Ch2_new Ch2 - average ... → 数据被永久改变了 投影重参考 (projectionTrue) 附加一条规则投影算子 不直接改数据 MNE 在需要时自动应用 → 数据保持原始状态更灵活 源定位必须用投影方式五、创建源空间5.1 什么是源空间源空间 大脑皮层上的一堆候选点 想象你在玩猜位置游戏 把大脑皮层表面画满小格子 每个格子是一个候选激活位置 源定位就是从这几千个格子中 找出哪些格子最可能产生了你测到的头皮电位 spacing 参数控制格子密度 oct6 → 每个半球约 4098 个点推荐 ico4 → 每个半球约 2562 个点较快 ico5 → 每个半球约 10242 个点精确但慢5.2 代码实现# ---------- 3. 创建源空间 ---------- print(\n创建源空间...) # setup_source_space()在皮层表面创建候选源点 src mne.setup_source_space( subjectsubject, # 被试名称 spacingoct6, # 源点密度oct6每半球约4098点 subjects_dirsubjects_dir, # MRI 数据存放目录 add_distFalse # 不计算源点间距离节省时间 ) # src 是一个列表包含左右两个半球 # src[0]左半球源空间 # src[1]右半球源空间 # vertno源点的顶点编号 print(f✅ 源空间: {len(src[0][vertno])} {len(src[1][vertno])} 源点) print(f 左半球: {len(src[0][vertno])} 个源点) print(f 右半球: {len(src[1][vertno])} 个源点)六、头模型BEM6.1 什么是 BEMBEM Boundary Element Method边界元方法 头模型描述电流如何穿过不同组织层 大脑皮层 → 脑脊液 → 颅骨 → 头皮 各层导电率不同 大脑导电好~0.3 S/m 颅骨导电差~0.006 S/m← 关键 头皮导电好~0.3 S/m 颅骨导电差 电流被阻挡和扩散 头皮电位被模糊化 这就是为什么头皮脑电空间分辨率不高约5-10cm6.2 代码实现# ---------- 4. 加载 BEM ---------- # BEM 模型文件路径 # 命名规则{被试名}-{头皮}-{颅骨内}-{颅骨外}-bem-sol.fif # 5120 是每层的三角形数量 bem_fname os.path.join( subjects_dir, subject, bem, f{subject}-5120-5120-5120-bem-sol.fif ) # read_bem_solution()加载预计算的 BEM 模型 bem mne.read_bem_solution(bem_fname) print(✅ BEM 模型加载完成)七、计算正向解7.1 什么是正向解正向解Forward Solution 传导矩阵Lead Field Matrix 对于每个候选源点计算 如果这里激活了头皮上每个电极会测到什么电位 结果是一个矩阵 源点 × 电极 传导关系 源点1激活 → 电极A:1.0μV, 电极B:0.5μV, 电极C:-0.3μV 源点2激活 → 电极A:0.3μV, 电极B:1.2μV, 电极C:0.1μV 源点3激活 → 电极A:-0.5μV, 电极B:0.8μV, 电极C:1.5μV ...7.2 代码实现# ---------- 5. 计算正向解 ---------- # trans 文件MRI 坐标系 ↔ 头部坐标系的转换矩阵 # 告诉 MNE大脑MRI 中的位置和电极头皮上的位置的关系 trans_fname os.path.join( sample_data_folder, MEG, sample, sample_audvis_raw-trans.fif ) trans mne.read_trans(trans_fname) # make_forward_solution()计算正向解 fwd mne.make_forward_solution( raw_eeg.info, # 通道信息电极位置、类型 transtrans, # 坐标转换矩阵 srcsrc, # 源空间候选激活点 bembem, # 头模型导电特性 eegTrue, # 使用 EEG 通道 megFalse, # 不使用 MEG 通道 mindist5.0, # 源点离内表面最小距离mm n_jobs1, # 并行核心数 verboseFalse # 不打印详细信息 ) print(f✅ 正向解: {fwd[nsource]} 源 × {fwd[nchan]} 通道)正向解的维度说明fwd[nsource] 8196 个源点左右半球总和 fwd[nchan] 60 个 EEG 通道 传导矩阵大小 60 × 8196 每一列 一个源点对所有通道的影响模式 每一行 一个通道对所有源点的敏感度八、计算噪声协方差8.1 为什么需要噪声协方差不同通道的噪声水平不同 电极1贴得很紧 → 噪声小 电极2有点松动 → 噪声大 通道间的噪声可能相关 相邻电极可能同时受同一噪声源影响 噪声协方差 描述这些噪声特性的矩阵 在逆解时白化数据 噪声大的通道 → 权重降低 噪声小的通道 → 权重升高 → 让所有通道平等地贡献8.2 代码实现# ---------- 6. 噪声协方差 ---------- # 提取事件 events mne.find_events(raw_eeg, stim_channelSTI 014) # 使用基线期刺激前的数据估计噪声 # tmin-0.5, tmax0只用刺激前 0.5 秒 # 这段时间没有刺激信号 纯粹的噪声 epochs_noise mne.Epochs( raw_eeg, events, event_id{听觉/左耳: 1}, tmin-0.5, tmax0, # 只用基线期 baselineNone, # 不做基线校正 preloadTrue, verboseFalse ) # compute_covariance()计算协方差矩阵 # methodempirical直接用数据计算 noise_cov mne.compute_covariance( epochs_noise, methodempirical, # 经验方法 verboseFalse ) print(✅ 噪声协方差计算完成) print(f 协方差矩阵形状: {noise_cov.data.shape})九、计算逆解9.1 什么是逆解逆解 从头皮电位反推大脑激活 MNE 方法 Minimum Norm Estimate最小范数估计 原理奥卡姆剃刀 在所有能解释头皮电位的源配置中 选择总能量最小的那个 最节俭的解释 就像 看到墙上有光斑 假设是最近的手电筒照的 而不是远处的高功率探照灯9.2 创建逆解算子# ---------- 7. 逆解算子 ---------- # make_inverse_operator()创建逆解算子 # 结合正向解和噪声协方差 inverse_operator mne.minimum_norm.make_inverse_operator( raw_eeg.info, # 通道信息 fwd, # 正向解传导矩阵 noise_cov, # 噪声协方差 loose0.2, # 源朝向约束 depth0.8, # 深度加权 verboseFalse ) print(✅ 逆解算子创建完成)参数详解十、应用逆解到数据# ---------- 8. 应用逆解 ---------- # 创建 Epochs event_id {听觉/左耳: 1} epochs mne.Epochs( raw_eeg, events, event_idevent_id, tmin-0.2, tmax0.5, baseline(-0.2, 0), rejectdict(eeg150e-6), preloadTrue, verboseFalse ) # 叠加平均得到 ERP evoked epochs[听觉/左耳].average() # apply_inverse()将逆解应用到 Evoked 数据 # 从头皮电位 → 大脑皮层激活 stc mne.minimum_norm.apply_inverse( evoked, # ERP 数据 inverse_operator, # 逆解算子 lambda21.0 / 9.0, # 正则化参数 methoddSPM, # dSPM 噪声归一化 verboseFalse ) print(f✅ 源估计完成形状: {stc.data.shape}) print(f 源点数 × 时间点数 {stc.data.shape})method参数选择十一、3D 源激活可视化# ---------- 9. 3D 可视化 ---------- print(\n绘制源激活3D 大脑窗口...) print(提示可以用鼠标旋转、缩放大脑) print( 关闭 3D 窗口后在控制台按回车退出程序) # stc.plot()在 3D 大脑模型上显示源激活 # 使用 pyvistaqt 创建交互式 3D 窗口 stc.plot( hemisplit, # 左右半球分开显示 subjects_dirsubjects_dir, # MRI 数据目录 subjectsubject, # 被试名称 initial_time0.12, # 初始显示 120ms 时间点 time_unitms, # 时间单位毫秒 climdict( # 颜色范围设置 kindvalue, # 基于数值 pos_lims[3, 6, 9] # [最小值, 中间值, 最大值] ) ) # 关键用 input() 保持程序运行 # stc.plot() 打开的是独立 3D 窗口Qt # 程序继续执行就会退出窗口随之关闭 # input() 让程序暂停等待用户输入 input(\n关闭 3D 窗口后按回车键退出程序...) print(\n *60) print(第9天学习完成) print(*60)十二、第9天完整代码# 环境设置 import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt import mne import numpy as np import os import warnings warnings.filterwarnings(ignore) plt.rcParams[font.sans-serif] [Microsoft YaHei, SimHei] plt.rcParams[axes.unicode_minus] False print(*60) print(MNE-Python 第9天源定位基础) print(*60) # ---------- 1. 获取数据路径 ---------- sample_data_folder mne.datasets.sample.data_path() subjects_dir os.path.join(sample_data_folder, subjects) subject sample raw_fname os.path.join(sample_data_folder, MEG, sample, sample_audvis_raw.fif) # ---------- 2. 加载数据并预处理 ---------- print(\n加载数据...) raw mne.io.read_raw_fif(raw_fname, preloadFalse) raw_eeg raw.copy().pick_types(eegTrue, eogTrue, stimTrue) eeg_names [ch for ch in raw_eeg.ch_names if ch.startswith(EEG)] montage mne.channels.make_standard_montage(standard_1020) standard_names montage.ch_names[:len(eeg_names)] raw_eeg.rename_channels(dict(zip(eeg_names, standard_names))) raw_eeg.set_montage(montage) raw_eeg.load_data() raw_eeg.notch_filter(freqs60, pickseeg, verboseFalse) raw_eeg.filter(l_freq1, h_freq40, pickseeg, verboseFalse) raw_eeg.set_eeg_reference(average, projectionTrue) print(✅ 数据预处理完成) # ---------- 3. 创建源空间 ---------- print(\n创建源空间...) src mne.setup_source_space( subjectsubject, spacingoct6, subjects_dirsubjects_dir, add_distFalse) print(f✅ 源空间: {len(src[0][vertno])} {len(src[1][vertno])} 源点) # ---------- 4. 加载 BEM ---------- bem_fname os.path.join(subjects_dir, subject, bem, f{subject}-5120-5120-5120-bem-sol.fif) bem mne.read_bem_solution(bem_fname) print(✅ BEM 模型加载完成) # ---------- 5. 计算正向解 ---------- trans_fname os.path.join(sample_data_folder, MEG, sample, sample_audvis_raw-trans.fif) trans mne.read_trans(trans_fname) fwd mne.make_forward_solution( raw_eeg.info, transtrans, srcsrc, bembem, eegTrue, megFalse, mindist5.0, verboseFalse) print(f✅ 正向解: {fwd[nsource]} 源 × {fwd[nchan]} 通道) # ---------- 6. 噪声协方差 ---------- events mne.find_events(raw_eeg, stim_channelSTI 014) epochs_noise mne.Epochs(raw_eeg, events, event_id{听觉/左耳: 1}, tmin-0.5, tmax0, baselineNone, preloadTrue, verboseFalse) noise_cov mne.compute_covariance(epochs_noise, methodempirical, verboseFalse) print(✅ 噪声协方差计算完成) # ---------- 7. 逆解算子 ---------- inverse_operator mne.minimum_norm.make_inverse_operator( raw_eeg.info, fwd, noise_cov, loose0.2, depth0.8, verboseFalse) print(✅ 逆解算子创建完成) # ---------- 8. 应用逆解 ---------- event_id {听觉/左耳: 1} epochs mne.Epochs(raw_eeg, events, event_idevent_id, tmin-0.2, tmax0.5, baseline(-0.2, 0), rejectdict(eeg150e-6), preloadTrue, verboseFalse) evoked epochs[听觉/左耳].average() stc mne.minimum_norm.apply_inverse( evoked, inverse_operator, lambda21.0/9.0, methoddSPM, verboseFalse) print(f✅ 源估计完成形状: {stc.data.shape}) # ---------- 9. 3D 可视化 ---------- print(\n绘制源激活3D 大脑窗口...) print(提示可以用鼠标旋转、缩放大脑) print( 关闭 3D 窗口后在控制台按回车退出程序) stc.plot( hemisplit, subjects_dirsubjects_dir, subjectsubject, initial_time0.12, time_unitms, climdict(kindvalue, pos_lims[3, 6, 9]) ) input(\n关闭 3D 窗口后按回车键退出程序...) print(\n *60) print(第9天学习完成) print(*60) print(\n 明日预告第10天 - 结果报告与可视化)十三、今日总结 核心概念️ 掌握的技能 预处理铁律源定位版1. 提取通道 pick_types() 2. 重命名通道 rename_channels() 3. 设置蒙太奇 set_montage() 4. 加载到内存 load_data() 5. 陷波滤波 notch_filter() 6. 带通滤波 filter() 7. 投影重参考 set_eeg_reference(average, projectionTrue) ← 必须用投影
