一、什么是源定位？
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') # pyplot：matplotlib 的绘图接口 # 虽然用 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(f"MRI 数据路径: {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, preload=False) # 提取通道：EEG + EOG + STIM # eeg=True：脑电通道（分析目标） # eog=True：眼电通道（可用于 ICA） # stim=True：刺激通道（提取事件需要） raw_eeg = raw.copy().pick_types(eeg=True, eog=True, stim=True)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(freqs=60, picks='eeg', verbose=False) # 带通滤波：保留 1-40 Hz # l_freq=1：高通，去除慢速漂移 # h_freq=40：低通，去除高频肌电噪声 raw_eeg.filter(l_freq=1, h_freq=40, picks='eeg', verbose=False) # 🔑 关键：用投影方式做平均参考 # projection=True：使用投影算子而非直接修改数据 # 为什么必须用投影方式？ # 源定位需要在计算过程中动态处理参考问题 # 直接修改数据会破坏这种灵活性 raw_eeg.set_eeg_reference('average', projection=True) print("✅ 数据预处理完成")projection=True的含义：
普通重参考 (projection=False)： 直接修改数据值 Ch1_new = Ch1 - average Ch2_new = Ch2 - average ... → 数据被永久改变了 投影重参考 (projection=True)： 附加一条"规则"（投影算子） 不直接改数据 MNE 在需要时自动应用 → 数据保持原始状态，更灵活 源定位必须用投影方式！五、创建源空间
5.1 什么是源空间？
源空间 = 大脑皮层上的一堆"候选点" 想象你在玩"猜位置"游戏： 把大脑皮层表面画满小格子 每个格子是一个"候选激活位置" 源定位就是从这几千个格子中 找出哪些格子最可能产生了你测到的头皮电位 spacing 参数控制格子密度： 'oct6' → 每个半球约 4098 个点（推荐） 'ico4' → 每个半球约 2562 个点（较快） 'ico5' → 每个半球约 10242 个点（精确但慢）5.2 代码实现
# ---------- 3. 创建源空间 ---------- print("\n创建源空间...") # setup_source_space()：在皮层表面创建候选源点 src = mne.setup_source_space( subject=subject, # 被试名称 spacing='oct6', # 源点密度（oct6=每半球约4098点） subjects_dir=subjects_dir, # MRI 数据存放目录 add_dist=False # 不计算源点间距离（节省时间） ) # 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'])} 个源点")六、头模型（BEM）
6.1 什么是 BEM？
BEM = Boundary Element Method（边界元方法） 头模型描述电流如何穿过不同组织层： 大脑皮层 → 脑脊液 → 颅骨 → 头皮 各层导电率不同： 大脑：导电好（~0.3 S/m） 颅骨：导电差（~0.006 S/m）← 关键！ 头皮：导电好（~0.3 S/m） 颅骨导电差 = 电流被阻挡和扩散 = 头皮电位被"模糊化" = 这就是为什么头皮脑电空间分辨率不高（约5-10cm）6.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, # 通道信息（电极位置、类型） trans=trans, # 坐标转换矩阵 src=src, # 源空间（候选激活点） bem=bem, # 头模型（导电特性） eeg=True, # 使用 EEG 通道 meg=False, # 不使用 MEG 通道 mindist=5.0, # 源点离内表面最小距离（mm） n_jobs=1, # 并行核心数 verbose=False # 不打印详细信息 ) 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_channel='STI 014') # 使用基线期（刺激前）的数据估计噪声 # tmin=-0.5, tmax=0：只用刺激前 0.5 秒 # 这段时间没有刺激，信号 = 纯粹的噪声 epochs_noise = mne.Epochs( raw_eeg, events, event_id={'听觉/左耳': 1}, tmin=-0.5, tmax=0, # 只用基线期 baseline=None, # 不做基线校正 preload=True, verbose=False ) # compute_covariance()：计算协方差矩阵 # method='empirical'：直接用数据计算 noise_cov = mne.compute_covariance( epochs_noise, method='empirical', # 经验方法 verbose=False ) 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, # 噪声协方差 loose=0.2, # 源朝向约束 depth=0.8, # 深度加权 verbose=False ) print("✅ 逆解算子创建完成")参数详解：
十、应用逆解到数据
# ---------- 8. 应用逆解 ---------- # 创建 Epochs event_id = {'听觉/左耳': 1} epochs = mne.Epochs( raw_eeg, events, event_id=event_id, tmin=-0.2, tmax=0.5, baseline=(-0.2, 0), reject=dict(eeg=150e-6), preload=True, verbose=False ) # 叠加平均得到 ERP evoked = epochs['听觉/左耳'].average() # apply_inverse()：将逆解应用到 Evoked 数据 # 从头皮电位 → 大脑皮层激活 stc = mne.minimum_norm.apply_inverse( evoked, # ERP 数据 inverse_operator, # 逆解算子 lambda2=1.0 / 9.0, # 正则化参数 method='dSPM', # dSPM = 噪声归一化 verbose=False ) 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( hemi='split', # 左右半球分开显示 subjects_dir=subjects_dir, # MRI 数据目录 subject=subject, # 被试名称 initial_time=0.12, # 初始显示 120ms 时间点 time_unit='ms', # 时间单位：毫秒 clim=dict( # 颜色范围设置 kind='value', # 基于数值 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, preload=False) raw_eeg = raw.copy().pick_types(eeg=True, eog=True, stim=True) 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(freqs=60, picks='eeg', verbose=False) raw_eeg.filter(l_freq=1, h_freq=40, picks='eeg', verbose=False) raw_eeg.set_eeg_reference('average', projection=True) print("✅ 数据预处理完成") # ---------- 3. 创建源空间 ---------- print("\n创建源空间...") src = mne.setup_source_space( subject=subject, spacing='oct6', subjects_dir=subjects_dir, add_dist=False) 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, trans=trans, src=src, bem=bem, eeg=True, meg=False, mindist=5.0, verbose=False) print(f"✅ 正向解: {fwd['nsource']} 源 × {fwd['nchan']} 通道") # ---------- 6. 噪声协方差 ---------- events = mne.find_events(raw_eeg, stim_channel='STI 014') epochs_noise = mne.Epochs(raw_eeg, events, event_id={'听觉/左耳': 1}, tmin=-0.5, tmax=0, baseline=None, preload=True, verbose=False) noise_cov = mne.compute_covariance(epochs_noise, method='empirical', verbose=False) print("✅ 噪声协方差计算完成") # ---------- 7. 逆解算子 ---------- inverse_operator = mne.minimum_norm.make_inverse_operator( raw_eeg.info, fwd, noise_cov, loose=0.2, depth=0.8, verbose=False) print("✅ 逆解算子创建完成") # ---------- 8. 应用逆解 ---------- event_id = {'听觉/左耳': 1} epochs = mne.Epochs(raw_eeg, events, event_id=event_id, tmin=-0.2, tmax=0.5, baseline=(-0.2, 0), reject=dict(eeg=150e-6), preload=True, verbose=False) evoked = epochs['听觉/左耳'].average() stc = mne.minimum_norm.apply_inverse( evoked, inverse_operator, lambda2=1.0/9.0, method='dSPM', verbose=False) print(f"✅ 源估计完成，形状: {stc.data.shape}") # ---------- 9. 3D 可视化 ---------- print("\n绘制源激活（3D 大脑窗口）...") print("提示：可以用鼠标旋转、缩放大脑") print(" 关闭 3D 窗口后在控制台按回车退出程序") stc.plot( hemi='split', subjects_dir=subjects_dir, subject=subject, initial_time=0.12, time_unit='ms', clim=dict(kind='value', 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', projection=True) ← 必须用投影！
