江苏省住房和城乡建设厅 官方网站企业网络搭建

---

前言

        传统的典型相关分析（CCA）方法在SSVEP信号识别中取得了一定成效，但容易受到自发脑电活动的干扰，且未充分利用相位信息。为此，研究者提出了TRCA方法，通过最大化任务期间神经影像数据的可重复性，提取任务相关成分，增强信号的信噪比（SNR），抑制背景脑电噪声。

---

一、TRCA是什么？

        在SSVEP-BCI系统中，任务相关成分分析（Task-Related Component Analysis，TRCA）被用于提高SSVEP信号的检测性能。通过增强任务相关成分，TRCA能够提高信号的信噪比，抑制背景脑电活动的干扰，从而提升系统的准确性和稳定性。

二、TRCA方法

        TRCA（任务相关成分分析）是一种通过最大化任务期间的可复现性来高效提取任务相关成分的方法。假设观测到的多通道 EEG 信号X(t)由以下线性生成模型表示：
$$x_j(t)=a_{1,j}s(t)+a_{2,j}n(t),j=1,2,\dots ,N_c$$其中，$j$ 是通道索引， $a_{1,j}$ 和 $a_{2,j}$ 是混合系数，用于将源信号投影到 EEG 观测信号中。

        问题在于如何从观测信号 x(t) 的线性组合中恢复任务相关成分 s(t)，即：
$$y(t)=\sum _{j=1}^{N_c}{w}_j{x}_j(t)=\sum _{j=1}^{N_c}\left(w_j{a}_{1,j}s(t)+w_j{a}_{2,j}n(t)\right)$$

        理想情况下，该问题的解满足：${\sum }_{j=1}^{N_c}{w}_j{a}_{1,j}=1$ 以及 ${\sum }_{j=1}^{N_c}{w}_j{a}_{2,j}=0$
从而得到最终解：
$$y(t)=s(t)$$

        该问题可以通过跨试次协方差最大化（inter-trial covariance maximization）来求解。
$$C_{h_1{h}_2}=\text{Cov}\left(y^{(h_1)}(t),y^{(h_2)}(t)\right)$$
       第 $h$ 次试验的 EEG 信号和估计的任务相关成分分别表示为
$$x^{(h)}(t)和y^{(h)}(t),其中h=1,\mathrm{2.....}{N}_t$$

       这里的$y^{(h)}(t)$的时间区间固定为 $t\in [t_h,t_h+T]$，此处 𝑇 是每个任务试次的持续时间。$y(t)$ 在试次$h_1$和 $h_2$ 之间的协方差表示：
$$C_{h_1{h}_2}=\text{Cov}\left(y^{(h_1)}(t),y^{(h_2)}(t)\right)=\sum _{j_1,j_2=1}^{N_c}{w}_{j_1}{w}_{j_2}\text{Cov}\left(x_{j_1}^{(h_1)}(t),x_{j_2}^{(h_2)}(t)\right)$$

所有可能的试次组合被总结如下：
$$\sum _{\begin{array}{c}{h}_1,h_2=1\\ h_1\ne h_2\end{array}}^{N_t}{C}_{h_1{h}_2}=\sum _{\begin{array}{c}{h}_1,h_2=1\\ h_1\ne h_2\end{array}}^{N_t}\sum _{j_1,j_2=1}^{N_c}{w}_{j_1}{w}_{j_2}\text{Cov}\left(x_{j_1}^{(h_1)}(t),x_{j_2}^{(h_2)}(t)\right)=w^{T}Sw$$

这里，矩阵$S=(S_{j_1{j}_2}{)}_{1\le j_1,j_2\le N_c}$会被定义为$S_{j_1{j}_2}={\sum }_{\begin{array}{c}{h}_1,h_2=1\\ h_1\ne h_2\end{array}}^{N_t}\text{Cov}\left(x_{j_1}^{(h_1)}(t),x_{j_2}^{(h_2)}(t)\right)$

为了获得有限解，$y(t)$的方差被约束为：

$$\text{Var}(y(t))=\sum _{j_1,j_2=1}^{N_c}{w}_{j_1}{w}_{j_2}\text{Cov}(x_{j_1}(t),x_{j_2}(t))=w^{T}Qw=1$$

该约束优化问题可以求解为：

$$\hat{w}=\arg \underset{w}{\max }\frac{w^{T}Sw}{w^{T}Qw}$$

最优的系数向量可以通过矩阵 $Q^{-1}S$的特征向量获得。

       一般来说，从一个 $N_c\times N_c$矩阵中可以得到 $N_c$个特征值和特征向量。矩阵$Q^{-1}S$的特征值 $\lambda$ 可以按降序排列，它表示相应特征向量 $\hat{w}$ 的成本函数值。因此，这些特征值反映了多个试次之间的任务一致性。

       图 (a) 显示了基于 TRCA 方法的流程图。在基于 SSVEP 的脑机接口（BCI）中，TRCA 可用于设计空间滤波器，以去除头皮记录中的背景 EEG 活动。
       针对第 $n$ 个刺激的空间滤波器， $w_n^{(m)}\in {\mathbb{R}}^{N_c}$ 通过 TRCA 从个体校准数据${\chi }_n^{(m)}$中获取，并在应用滤波器组分析后计算得到。在这里，公式$\text{Var}(y(t))=w^{T}Qw=1$中的 $Q$ 矩阵是使用所有训练试次的拼接矩阵${\chi }_n^{(m)}$计算得到的。

       单试次测试数据 $X^{(m)}\in {\mathbb{R}}^{N_c\times N_s}$与跨试次平均的第$n$个视觉刺激训练数据${\bar{\chi }}_n^{(m)}\in {\mathbb{R}}^{N_c\times N_s}$之间的相关系数计算如下：
$$T{r}_n^{(m)}=\rho \left({\left(X^{(m)}\right)}^T{w}_n^{(m)},{\bar{\chi }}_n^{(m)}{w}_n^{(m)}\right)$$
其中，$p(a,b)$ 表示信号 $a$ 和 $b$ 之间的皮尔逊相关分析。

---

三、matlab实战

1.数据集

使用清华大学脑机接口团队的公开数据集，benchmark dataset
下载地址： https://bci.med.tsinghua.edu.cn/
关于该数据集的介绍和使用，可以看这篇文章：
清华大学脑机接口研究组Benchmark dataset数据集介绍

2. test_trca.m

trca的测试代码如下，主要是读取数据文件，然后对数据滤波处理后，送入trca算法。

% 加载所需的数据文件
load ('Freq_Phase.mat');  % 载入频率和相位数据
data = load('S1.mat');    % 载入EEG数据
eeg_data = data.data;    % 获取EEG数据矩阵% 处理数据
% 原始数据，维度是 [64, 1500, 40, 6]（64 个电极，1500 个时间点, 40目标刺激, 6个试次）
% 提取顶叶和枕叶区域（Pz、PO5、PO3、POz、PO4、PO6、O1、Oz和O2）的SSVEP信号
stim_data = eeg_data([48, 54, 55, 56, 57, 58, 61, 62, 63], 251:1250, :, :);
% 对stim_data进行5-90Hz带通滤波
% stim_data的维度是 [9, 1500, 40, 6]
% 设置滤波器参数
low_freq = 5;   % 下边界频率
high_freq = 90; % 上边界频率
fs = 250;      % 采样频率为250Hz% 获取数据的维度
[N_channels, N_timepoints, N_targets, N_blocks] = size(stim_data);% 对每个数据块和每个目标进行滤波
for block_idx = 1:N_blocksfor target_idx = 1:N_targetsfor channel_idx = 1:N_channels% 提取每个通道、目标和数据块的信号signal = squeeze(stim_data(channel_idx, :, target_idx, block_idx));% 使用bandpass滤波器对信号进行5-90Hz滤波stim_data(channel_idx, :, target_idx, block_idx) = bandpass(signal, [low_freq, high_freq], fs);endend
end
eeg = stim_data;
% 获取EEG数据的维度：通道数、时间点数、目标刺激数、数据块数
[N_channel,~, N_target, N_block] = size(eeg);%% ------------分类部分-------------
tic
% 留一交叉验证（LOO Cross-Validation）
for loocv_i = 1:N_block% 提取测试数据（留一交叉验证中，当前块作为测试数据）Testdata = squeeze(eeg(:, :, :, loocv_i));% 提取训练数据（去除当前块数据作为训练集）Traindata = eeg;Traindata(:, :, :, loocv_i) = [];% 对每个目标频率，计算训练数据的平均值for targ_i = 1:N_targetaver_Traindata(:, :, targ_i) = squeeze(mean(squeeze(Traindata(:,:,targ_i,:)),3));  % 对每个目标频率进行平均end % end targ_i% 标签赋值，真实标签是频率truelabels = freqs;% 获取测试数据的试次数量N_testTrial = size(Testdata, 3);% 对每个测试试次进行处理for trial_i = 1:N_testTrialcoefficience = zeros(1, length(truelabels));  % 初始化系数矩阵% 对每个频率进行处理for targ_j = 1:length(freqs)             % 使用训练数据计算空间滤波器 wnwn = TRCA(squeeze(Traindata(:,:,targ_j,:)));  % 调用TRCA函数计算空间滤波器% 计算加权训练数据与测试数据的相关性weighted_train = wn' * aver_Traindata(:,:,targ_j);  % 加权训练数据weighted_test = wn' * Testdata(:,:,trial_i);  % 加权测试数据% 计算测试数据和训练数据之间的相关系数coefficienceMatrix = corrcoef(weighted_test, weighted_train);coefficience(targ_j) = coefficienceMatrix(1, 2);  % 获取相关系数end % end targ_i% 目标检测：选择具有最大相关系数的频率作为输出标签[~, index] = max(coefficience);  % 获取最大相关系数对应的频率outputlabels(trial_i) = freqs(index);  % 输出预测的标签end % end trial_i% 计算该次测试的正确分类数量trueNum = sum((outputlabels - truelabels) == 0);  % 统计预测标签与真实标签相同的个数acc(loocv_i) = trueNum / length(truelabels);  % 计算该次交叉验证的准确率% 打印当前交叉验证的准确率fprintf('第%d试次的交叉验证准确率是: %.4f, 样本数: %d/%d\n', loocv_i, acc(loocv_i), trueNum, N_testTrial)
end % end looCv_i% 计算总的运行时间
t = toc;
% 打印总的运行时间和平均准确率
fprintf('\n-----------------------------------------\n')
disp(['总运行时间: ', num2str(t), ' s']);
fprintf('6试次数据交叉验证的平均准确率是: %.4f\n', mean(acc));

3. trca.m

trca的计算代码

% -------------------------------------------------------------------
% Task-related component analysis[1].extract task-related components 
% by maximizing the reproducibility during the task period.
%
% Input:         
%   eeg : input eeg training data
%         (# channels, # points, # trials)
% Output:  
%   W: w means eigenvectors corresponding to max eigenvalues
%
% Reference:
%   [1] H. Tanaka, T. Katura, H. Sato,
%       "Task-related component analysis for functional neuroimaging and 
%        application to near-infrared spectroscopy data",
%       NeuroImage, vol. 64, pp. 308-327, 2013.
% --------------------------------------------------------------------
function [w] = TRCA(eeg)% 输入：eeg：一个三维矩阵，包含 EEG 数据（通道数 × 样本数 × 试次数）% 输出：w：任务相关成分的空间滤波器（特征向量）% 获取 EEG 数据的维度[num_chans, num_smpls, num_trials] = size(eeg);% 初始化协方差矩阵 S 为零矩阵S = zeros(num_chans);% 计算协方差矩阵 S，双重循环遍历每一对试次for trial_i = 1:1:num_trials-1  % 遍历第一个试次x1 = squeeze(eeg(:,:,trial_i));  % 提取第 trial_i 试次的数据（去掉第三维）x1 = bsxfun(@minus, x1, mean(x1, 2));  % 对每个通道去均值，避免偏差for trial_j = trial_i+1:1:num_trials  % 遍历第二个试次（只考虑 trial_i 和 trial_j 的组合）x2 = squeeze(eeg(:,:,trial_j));  % 提取第 trial_j 试次的数据x2 = bsxfun(@minus, x2, mean(x2, 2));  % 对每个通道去均值% 计算协方差矩阵的增量，并更新 SS = S + x1*x2' + x2*x1';end % end trial_jend % end trial_i% 将 EEG 数据重塑为一个二维矩阵，每一列表示一个试次的所有样本UX = reshape(eeg, num_chans, num_smpls*num_trials);% 对重塑后的数据进行去均值处理UX = bsxfun(@minus, UX, mean(UX, 2));% 计算 UX 的自相关矩阵 QQ = UX*UX';% 通过广义特征值问题求解空间滤波器 w[V, ~] = eigs(S, Q);  % 求解 S 和 Q 的广义特征值问题，V 为特征向量% 选择最大特征值对应的特征向量作为空间滤波器w = V(:, 1);  % 返回第一个特征向量作为空间滤波器（任务相关成分）end %end function

4.运行结果

---

总结

TRCA作为一种有效的任务相关信号提取方法，在SSVEP信号检测中发挥了重要作用。通过最大化任务相关成分的可重复性，TRCA提高了信号的信噪比，增强了SSVEP信号的检测性能。与其他方法相比，TRCA在处理背景脑电噪声方面表现出色，为SSVEP-BCI系统的优化提供了新的思路。

---

参考文献

M. Nakanishi, Y. Wang, X. Chen, Y. -T. Wang, X. Gao and T. -P. Jung, “Enhancing Detection of SSVEPs for a High-Speed Brain Speller Using Task-Related Component Analysis,” in IEEE Transactions on Biomedical Engineering, vol. 65, no. 1, pp. 104-112, Jan. 2018, doi: 10.1109/TBME.2017.2694818.