Beamformer 溯源实战：空域滤波器

一、Beamformer 直觉：给每个体素装一台"定向麦克风"

还是用"听声"类比：会议室里有几十个人同时说话，你手里有几十个话筒（电极），每个话筒录到的都是所有人的混音。 定向麦克风的思路是：把所有话筒的信号按不同权重加起来，让某个方向的声音原样通过、其他方向被抵消。

Beamformer（空域滤波器）对大脑做同样的事：对每个候选体素 v，从数据里学出一组权重 w_v， 使得 w_v 加权后的传感器信号"只反映体素 v 的活动"。这组权重与传感器信号的乘积， 就是该体素的虚拟电极（virtual electrode）——仿佛你把一根电极直接插进了那个位置。

关键点在于：权重完全从数据协方差矩阵里学来——哪些电极组合是噪声、哪些组合像信号， 不是人为规定，而是数据自己教给算法的。这也是 beamformer 家族名称的由来： 用波束（beam）的指向性"扫描"全脑。

二、LCMV：协方差驱动的空域滤波器

LCMV（Linearly Constrained Minimum Variance，线性约束最小方差）是最常用的 beamformer 变体。它的滤波器要同时满足两条：

1. 约束：该位置的单位源必须"原样通过"——滤波器不能扭曲目标位置的信号幅度；
2. 目标：输出总功率最小——在满足约束的前提下，把"其他位置 + 噪声"的贡献压到最低。

直觉版公式（不必推导）：w = C⁻¹k / (kᵀC⁻¹k)。 其中 C 是数据协方差矩阵（"各电极如何一起波动"的统计，来自你的试次）， k 是该体素的前置场列向量（"这个位置的源长什么样"）。 一句话：用噪声结构的逆，去乘上这个位置的特征，再做归一。 分母的归一化正是为了满足上面第 1 条约束。

两个必懂的延伸：

• 正则化：C 求逆在病态时会爆炸（见第六节），所以实践里总是加 lambda 正则化，FieldTrip 写作 cfg.lcmv.lambda = '5%'（按协方差迹的百分比自动定）；
• 深度偏置：beamformer 也有（浅部"更容易通过"），可用 cfg.lcmv.weightnorm = 'unitnoisegain' 等权重归一化缓解。

三、MNE vs Beamformer：什么时候用哪个

四、FieldTrip 实战：协方差 → LCMV → 虚拟电极

头模型与前置场直接复用 7.2 存好的那套。先算全脑激活图，再在感兴趣位置提取虚拟电极：

% ---- 第 1 步：估计数据协方差矩阵（beamformer 的燃料）----
cfg = [];
cfg.trials = 'all';
cfg.channel = 'EEG';
cfg.covariance = 'yes';                    % 让 ft_timelockanalysis 顺手算协方差
cfg.covariancewindow = 'all';              % 也可 [-inf 0]：只用刺激前基线
cfg.demean = 'yes';                        % 协方差对慢漂移极其敏感
timelock = ft_timelockanalysis(cfg, data); % timelock.cov 即协方差矩阵

% ---- 第 2 步：全脑 LCMV 激活图 ----
cfg = [];
cfg.method = 'lcmv';
cfg.sourcemodel = leadfield;               % 复用 7.2 的头模型 + 前置场
cfg.headmodel = headmodel;
cfg.elec = timelock.elec;
cfg.lcmv.lambda = '5%';                    % 协方差正则化：防病态（见第六节）
source = ft_sourceanalysis(cfg, timelock);
% source.avg.pow 是全脑激活图，用 7.2 学的方法定位感兴趣的位置

% ---- 第 3 步：在两个种子点提取虚拟电极 ----
% peakA / peakB：第 2 步激活图上选定的体素编号（或解剖先验位置）
seedpos = source.pos([peakA, peakB], :);   % 2x3 坐标

cfg = [];
cfg.method = 'lcmv';
cfg.sourcemodel.pos = seedpos;             % 只含 2 个种子点的"迷你源模型"
cfg.headmodel = headmodel;
cfg.elec = timelock.elec;
cfg.lcmv.keepfilter = 'yes';               % 保留滤波器，下一步要用
cfg.lcmv.fixedori = 'yes';                % 每个种子点只留一条主朝向时序
cfg.lcmv.lambda = '5%';
seed = ft_sourceanalysis(cfg, timelock);    % 在平均数据上算好滤波器

% ---- 第 4 步：把滤波器应用到每个试次 → 分试次虚拟电极 ----
cfg = [];
cfg.method = 'lcmv';
cfg.sourcemodel.pos = seedpos;
cfg.sourcemodel.filter = seed.avg.filter;  % 复用刚算好的滤波器（官方推荐写法）
cfg.headmodel = headmodel;
cfg.elec = timelock.elec;
cfg.keeptrials = 'yes';                    % 保留每试次结果
vsource = ft_sourceanalysis(cfg, data);    % 输入分试次数据
% vsource.trial{n}.mom 为 [2 x 时间]：两条虚拟电极在每个试次的时序

五、虚拟电极与连接分析：呼应模块 6.3

模块 6.3 结尾我们留了一个伏笔：头皮电极之间的"连接"有很大一部分是容积传导的假象—— 同一个源传到两个电极，它们当然高度相关。正确的思路是先源后连： 先把源定位做好，再在源与源之间算连接。虚拟电极正是这条路的入场券。

% ---- 第 5 步：虚拟电极进入模块 6 的连接分析 ----
% 先把虚拟电极包装回 FieldTrip 的 raw 数据结构
vdata = [];
vdata.label = {'seedA', 'seedB'};
vdata.fsample = data.fsample;
for n = 1:numel(vsource.trial)
    vdata.trial{n} = vsource.trial(n).mom;   % [2 x 时间]
    vdata.time{n}  = data.time{n};          % 时间轴与原试次一致
end

% 然后走模块 6.3 的老路：时频 → 连接
cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'fourier';
cfg.keeptrials = 'yes';
cfg.foilim = [8 13];                        % 例：alpha 频段
vfreq = ft_freqanalysis(cfg, vdata);

cfg = [];
cfg.method = 'wpli';
vconn = ft_connectivityanalysis(cfg, vfreq);  % 两个源之间的 wPLI

六、常见坑：病态协方差、参考依赖、信号抵消

1. 协方差矩阵病态：通道多、数据少时，协方差矩阵接近奇异，求逆数值爆炸，表现为激活图上出现离谱的尖峰或全 0。
   • 判据：条件数（cond(timelock.cov)）过大（经验上超过 1e12 很危险）；
   • 对策：增加试次/时长；提高正则化（cfg.lcmv.lambda 从 '5%' 升到 '10%'）；检查是否有坏通道混进协方差。
2. 参考依赖：beamformer 的滤波器建立在协方差上，而 EEG 的协方差依赖参考电极选择——换参考等于换数据。 传感器水平相对宽容的参考选择问题，在 beamformer 这里会被放大。用平均参考（模块 2.2）是最稳妥的做法，且不同被试、不同条件要用同一参考。
3. 相关源的"信号抵消"：beamformer 的目标函数假设"别处的活动都是要压制的干扰"。 当两个源高度相关（如双侧镜像放电）时，它们的线性组合会被滤波器当成"公共干扰"一起压制——真实存在的双侧激活在图上消失。
   • 自查办法：与 MNE 结果交叉对比（MNE 不做独立性假设）；分别用单侧数据重算协方差看激活是否"回来"；
   • 综述与癫痫文献里这是 beamformer 最著名的失败模式，报告阴性结论前务必排除。