脑电5目录1. EEG的逆问题2. 偶极子源定位3. 分布源成像4. 为什么需要单试次分析5. 单试次分析的方法6. 用trial-by-trial数据解码认知状态1. EEG的逆问题1.1 正问题 vs 逆问题importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignal# 直观理解正问题和逆问题的区别print( 正问题和逆问题 )print(正问题容易唯一解已知脑内源的位置和强度 → 计算头皮电位分布)print(逆问题困难无穷多解已知头皮电位分布 → 反推脑内源的位置和强度)print()print(打个比方)print( 正问题 知道石头扔进池塘的位置 → 预测水波传到岸边的样子)print( 逆问题 看到岸边水波的样子 → 反推石头扔在哪个位置)print( 逆问题有无穷多组可能的答案——你需要额外的约束条件)为什么逆问题难同样的头皮电位分布可以由数量、位置、方向都不同的脑内源组合产生。这是数学上固有的不适定性——无穷多组脑内源能产生完全相同的头皮EEG。1.2 解决逆问题的额外约束约束方法思路等效偶极子假设源的数量很少1-3个用等效电流偶极子近似分布源模型假设源分布在整个皮层表面但位置固定在灰质LORETA/sLORETA分布源 空间平滑约束——相邻源强度相似波束成形(Beamformer)对每个位置的源单独估计抑制其他位置的干扰核心权衡偶极子模型稀疏但精确 / 分布源模型全面但模糊。2. 偶极子源定位2.1 等效电流偶极子简化的物理模型把一大片同步激活的锥体神经元看作一个等效电流偶极子。一个偶极子由6个参数描述位置(x,y,z)、方向(2个角度)、强度。defdipole_forward(dipole_pos,dipole_moment,electrode_pos): 计算单个偶极子在头皮电极上产生的电位简化球形头模型 参数: dipole_pos: 偶极子位置 (3,) dipole_moment: 偶极子矩向量 (3,)方向和强度 electrode_pos: 电极位置 (n_elec, 3) 返回: potentials: (n_elec,) 各电极的理论电位 n_elecelectrode_pos.shape[0]potentialsnp.zeros(n_elec)foriinrange(n_elec):# 偶极子到电极的向量relectrode_pos[i]-dipole_pos distancenp.linalg.norm(r)1e-6# 简化公式电位 ∝ (偶极矩·r) / r³potentials[i]np.dot(dipole_moment,r)/(distance**3)returnpotentials# 演示同一偶极子在不同位置的电位分布n_elec32anglesnp.linspace(0,2*np.pi,n_elec,endpointFalse)elec_posnp.column_stack([np.cos(angles)*0.8,np.sin(angles)*0.8,np.zeros(n_elec)])# 两个不同位置的偶极子dipole1np.array([0.3,0.0,0.3])# 偏右前dipole2np.array([-0.15,-0.35,0.3])# 偏左后momentnp.array([0,-1,0])# 向下朝向皮层深处pot1dipole_forward(dipole1,moment,elec_pos)pot2dipole_forward(dipole2,moment,elec_pos)fig,axesplt.subplots(1,2,figsize(12,4.5),subplot_kw{projection:polar})forax,pot,title,posin[(axes[0],pot1,右侧偶极子,dipole1),(axes[1],pot2,左侧偶极子,dipole2)]:ax.fill(angles,pot-pot.min()0.5,alpha0.6)ax.set_title(f{title}\n位置({pos[0]:.1f},{pos[1]:.1f}));ax.set_yticks([])plt.suptitle(偶极子源定位原理不同位置的源→不同的头皮电位分布,fontsize13)plt.tight_layout();plt.show()print(观察偶极子位置变了头皮电位分布明显不同——这就是源定位的基础。)2.2 偶极子拟合已知头皮电位 → 不断调整偶极子的6个参数 → 使理论电位和实际电位的误差最小 → 找到最佳偶极子。本质上是一个非线性优化问题。常用算法Levenberg-Marquardt阻尼最小二乘、单纯形法。deffit_dipole(eeg_topo,electrode_pos,n_dipoles1): 简化版偶极子拟合演示原理非实际可用 真实应用需要真实头模型BEM/FEM、多层组织电导率 fromscipy.optimizeimportminimizedefcost(params):# params: [x, y, z, mx, my, mz] for each dipoletotal_potnp.zeros(electrode_pos.shape[0])fordinrange(n_dipoles):idxd*6posparams[idx:idx3]momparams[idx3:idx6]total_potdipole_forward(pos,mom,electrode_pos)returnnp.sum((total_pot-eeg_topo)**2)# 随机初始化initnp.random.randn(n_dipoles*6)*0.2resultminimize(cost,init,methodNelder-Mead)returnresult.x.reshape(n_dipoles,6)print(拟合本质minimize ||V_forward(dipole_params) - V_measured||²)print(关键挑战初始值敏感 容易陷入局部最优 源数未知)3. 分布源成像3.1 从几个源到几千个源偶极子模型假设只有1-3个源适合癫痫棘波等局灶性活动但认知过程通常涉及广泛分布的脑网络。分布源成像把大脑皮层表面划分为几千个小网格每个网格是一个迷你偶极子反推所有网格的激活强度。3.2 sLORETA——标准的分布源成像方法sLORETA标准化低分辨率电磁断层成像是应用最广泛的分布源成像工具。核心公式J ^ T ⋅ Φ \hat{J} T \cdot \PhiJ^T⋅Φ其中Φ是头皮电位T是逆算子矩阵。# sLORETA的简化演示核心思想defdemo_sloreta_concept():演示分布源成像vs偶极子模型的区别fig,axesplt.subplots(1,2,figsize(12,4.5))# 偶极子模型稀疏1-3个精确的点axes[0].scatter([0.3,-0.2],[0.1,-0.3],s[300,200],c[red,blue],alpha0.7,edgecolorsblack)axes[0].set_xlim(-0.5,0.5);axes[0].set_ylim(-0.5,0.5)axes[0].set_title(偶极子模型1-3个等效源);axes[0].set_aspect(equal)# 分布源模型整个皮层表面都是源np.random.seed(0)grid_xnp.random.uniform(-0.5,0.5,500)grid_ynp.random.uniform(-0.5,0.5,500)activationsnp.exp(-((grid_x-0.2)**2(grid_y-0.1)**2)/0.03)*0.8activationsnp.exp(-((grid_x0.15)**2(grid_y0.2)**2)/0.04)*0.6axes[1].scatter(grid_x,grid_y,cactivations,s15,cmaphot,alpha0.8)axes[1].set_xlim(-0.5,0.5);axes[1].set_ylim(-0.5,0.5)axes[1].set_title(分布源模型数千个皮层网格点);axes[1].set_aspect(equal)plt.suptitle(偶极子 vs 分布源稀疏精确 vs 全面模糊,fontsize13)plt.tight_layout();plt.show()demo_sloreta_concept()3.3 源分析的实际应用场景应用方法选择典型场景癫痫灶定位偶极子(1-2个)发作间期棘波的精准定位认知ERP的源分布源(sLORETA)P3/N400等成分的脑内源分布静息态网络分布源 功能连接默认模式网络、突显网络实时神经反馈波束成形实时估计ROI的源活动冥想状态的脑网络分布源(sLORETA/eLORETA)冥想时α源从枕叶→前额叶转移注源分析不是找到唯一的真实源而是在合理约束下给出最可能的估计。4. 为什么需要单试次分析4.1 叠加平均丢掉了什么传统ERP分析 把同一条件的所有trial叠加平均。这假设了每个trial的脑电响应是相同的。但反应时变异同一被试对同一刺激有时快有时慢脑电响应的时间也相应变化注意波动被试可能走神了几个trial这些trial的ERP实际上是不同的学习效应实验中段和后段的ERP可能因练习而改变个体策略差异不同被试完成同一任务的认知策略可能不同叠加平均把所有这些变异都抹平了。单试次分析就是不去平均直接分析每个trial的特征。4.2 单试次分析能回答的问题问题类型举例脑-行为相关“P3幅度越大的trial反应时越快吗”时间动态“α功率从实验开始到结束如何变化”脑状态解码“仅凭这一个trial的EEG能判断被试看到了什么图片吗”预测建模“刺激出现前的α相位能预测被试能否检测到微弱刺激吗”5. 单试次分析的方法5.1 三大方法路线单试次分析 ├── 时域逐trial提取ERP峰值/均值 ├── 时频域逐trial计算ERSP → 试次间的频谱变异 └── 解码/分类用ML模型从单trial EEG预测行为/刺激类别5.2 时域单试次特征提取defextract_single_trial_features(epochs,times,fs256): 从每个trial提取时域和频域特征——替代叠加平均 参数: epochs: (n_trials, n_channels, n_times) 返回: features: (n_trials, n_features) 每行一个trial的特征向量 fromscipy.signalimportwelch n_trials,n_chepochs.shape[0],epochs.shape[1]features_list[]foriinrange(n_trials):trialepochs[i]# (n_ch, n_times)f{}# 时域P3窗口(250-500ms)的均值p3_mask(times0.25)(times0.5)forchinrange(min(n_ch,8)):# 每个通道f[fch{ch}_P3_mean]trial[ch,p3_mask].mean()f[fch{ch}_P3_peak]trial[ch,p3_mask].max()# 频域α功率forchinrange(min(n_ch,8)):f_psd,psdwelch(trial[ch],fs,npersegfs)alpha_mask(f_psd8)(f_psd13)f[fch{ch}_alpha]np.trapz(psd[alpha_mask],f_psd[alpha_mask])features_list.append(f)returnfeatures_list5.3 单试次解码——用EEG预测被试状态fromsklearn.svmimportSVCfromsklearn.model_selectionimportcross_val_score,StratifiedKFoldfromsklearn.preprocessingimportStandardScalerdefdecode_from_single_trials(epochs,labels,times,fs256): 从单试次EEG解码刺激类别或行为反应 这就是脑机接口(BCI)的基本范式 # 提取特征feat_dictsextract_single_trial_features(epochs,times,fs)Xnp.array([[d[k]forkinsorted(d.keys())]fordinfeat_dicts])# 标准化scalerStandardScaler()X_scaledscaler.fit_transform(X)# 分类clfSVC(kernelrbf,C1.0,random_state42)scorescross_val_score(clf,X_scaled,labels,cvStratifiedKFold(5,shuffleTrue))print(f单试次解码准确率:{scores.mean():.3f}±{scores.std():.3f})print(f随机基线:{1/len(np.unique(labels)):.3f})returnscores.mean(),clf,scaler6.用trial-by-trial数据解码认知状态# 完整实验模拟解码 np.random.seed(42)fs256timesnp.arange(-0.2,0.8,1000/fs)n_ch,n_trials8,200# 生成两类刺激的epochA类强P3, B类弱P3epochsnp.zeros((n_trials,n_ch,len(times)))labelsnp.zeros(n_trials,dtypeint)foriinrange(n_trials):is_class_Ain_trials//2labels[i]int(is_class_A)alpha12*np.sin(2*np.pi*10*timesnp.random.randn()*0.5)noisenp.random.randn(n_ch,len(times))*3# P3A类强(8μV)B类弱(3μV)p3_amp8ifis_class_Aelse3np.random.randn()*1.5p3p3_amp*np.exp(-((times-0.32)/0.06)**2)forchinrange(n_ch):epochs[i,ch]alpha0.6*noise[ch]p3*(1.5-0.15*ch)# 单试次解码acc,clf,scalerdecode_from_single_trials(epochs,labels,times,fs)# 对比叠加平均后的ERPerp_Aepochs[labels1].mean(axis0).mean(axis0)# A类平均erp_Bepochs[labels0].mean(axis0).mean(axis0)# B类平均fig,axesplt.subplots(1,3,figsize(16,4.5))# 单试次前20个foriinrange(10):axes[0].plot(times,epochs[i,0]i*15,blue,lw0.4,alpha0.6)axes[0].plot(times,epochs[i100,0](i10)*15,red,lw0.4,alpha0.6)axes[0].set_title(单试次EEG前10个A类10个B类\n蓝色A(强P3)红色B(弱P3))axes[0].set_xlabel(s);axes[0].set_yticks([]);axes[0].grid(True,alpha0.3)# 叠加平均ERPaxes[1].plot(times,erp_A,blue,lw2,labelA类平均)axes[1].plot(times,erp_B,red,lw2,labelB类平均)axes[1].axvline(0,colork,ls--,alpha0.3)axes[1].set_title(f叠加平均ERP\nP3差异明显但对trial间变异视而不见)axes[1].set_xlabel(s);axes[1].legend();axes[1].grid(True,alpha0.3)# 解码准确率axes[2].bar([随机,解码],[0.5,acc],color[gray,green],edgecolorblack)axes[2].set_ylim(0,1);axes[2].set_title(f单试次解码:{acc:.2%})axes[2].set_ylabel(准确率)fori,vinenumerate([0.5,acc]):axes[2].text(i,v0.03,f{v:.1%},hacenter,fontsize12,fontweightbold)plt.suptitle(单试次分析不叠加平均直接分析每个trial——解码认知状态,fontsize13)plt.tight_layout();plt.show()print(f\n结论解码准确率{acc:.1%} 随机50%)print(说明单试次EEG确实包含足够信息来区分两类刺激——叠加平均掩盖了这些信息。)单试次分析的核心价值发现叠加平均看不出的trial-by-trial变异和脑-行为相关性是脑机接口(BCI)和实时神经反馈系统的技术基础配合解码模型可实现仅凭单个trial的EEG预测被试看到了什么/在想什么