摘要：同时获取脑电图和功能磁共振成像（EEG-fMRI）允许以高时间和空间分辨率对大脑的电生理和血流动力学进行互补研究。其中一个具有巨大潜力的应用是基于实时分析脑电图和功能磁共振成像信号进行目标脑活动的神经反馈训练。这依赖于实时减少严重伪迹对脑电图信号的影响，主要是梯度和脉冲伪迹。已经提出了一些方法来实现这个目的，但它们要么速度慢、依赖特定硬件、未公开或是专有软件。在这里，我们介绍了一种完全开源且公开可用的工具，用于同时进行脑电图和功能磁共振成像记录中的实时脑电图伪迹去除，它速度快且适用于任何硬件。我们的工具集成在Python工具包NeuXus中。我们在三个不同数据集上对NeuXus进行了基准测试，评估了伪迹功率减少和静息状态下背景信号保留、闭眼时α波带功率反应以及运动想象事件相关去同步化的能力。我们通过报告执行时间低于250毫秒证明了NeuXus的实时能力。总之，我们提供并验证了第一个完全开源且与硬件无关的解决方案，用于实时去除同时进行的脑电图和功能磁共振成像研究中的伪迹。
1. 引言
自1999年以来，同时获取脑电图和功能磁共振成像(EEG-fMRI)被用于研究脑功能。这种多模态技术通过结合电活动的精细时间分辨率测量(EEG)和大脑血流动力学的精细空间分辨率测量(fMRI)，为神经元过程提供了两种不同的观点。其中最大的挑战之一是在MR环境下对EEG数据产生的伪影。振幅最大的是梯度伪影(GA)，它是由与电极和受试者头部形成的电路上的图像采集相关的时变磁场引起的。反过来，脉冲伪影(PA)是由心脏周期发生的变化引起的，包括体积头由于心脏将血液泵向头部而导致的运动，由于动脉血管脉动而导致的头皮扩张，以及在磁场存在下在血液离子之间产生电压的霍尔效应。这个伪影已经存在于静态磁场中，即使没有获得MRI。与GA相比，PA表现出更不规则的模式，跟随心脏周期。除了GA和PA外，MRI扫描仪内记录的脑电图还受到自发头部运动以及环境的影响，包括电力线、灯光、通风和氦气泵。
为了在后处理(离线)中减少这些脑电信号伪影，已经开发了几种方法。其中大多数是基于人工平均模板减法(AAS)技术，其工作原理是对具有周期性伪影的信号的一组片段进行平均，以获得伪影模板并从信号中减去它。其他方法包括盲源分离技术，如主成分分析(PCA)或独立成分分析(ICA)，滤波或字典学习方法。另外，也有人提出了基于使用附加传感器测量伪信号的方法。这可以通过使用一个完整的参考电极层来直接实现，该参考电极层与脑电图电极重复，但与大脑绝缘。在一种相对更实用的方法中，可以将有限数量的传感器放置在特定位置(或简单地通过将几个脑电图电极与头皮绝缘来获得)，然后可以将每个电极位置的运动伪像估计为传感器信号时间过程的线性组合。尽管使用基于传感器的技术可以更有效地减少PA和其他运动伪影，尽管碳丝传感器最近已经商业化用于此目的，但大多数EEG-fMRI研究仍然没有使用额外的传感器。尽管提出了各种各样的方法来减少EEG-fMRI中的伪影，但大多数方法都不能实时(在线)使用，因为它们依赖于对完整信号的分析，并且通常需要大量的计算。
在线减少脑电图伪影对于诸如在实验期间监测大脑活动(例如，癫痫活动)或提供神经反馈(例如，训练目标大脑活动)等应用至关重要。已经提出了一些方法来实时减少脑电信号伪影。唯一的商用工具(BrainVision RecView)使用AAS来减少GA和PA，但它是专有的。Shaw等人提出的过滤方法过滤了遗传算法，但引入了时间扭曲，在公开可用之前需要进一步验证。Wu等人提出的PCA方法降低了PA，也不可用。Mayeli 等人中提出的ICA方法仍然需要事先进行GA和PA缩减，并引入了2到4秒的延迟。最后，提出了一些基于传感器的方法来实时减少遗传算法和PA。然而，这些需要额外的(主要是定制的)硬件，而大多数EEG-fMRI研究都不具备这些硬件。在这里，我们提出了一个完全开源的工具，用于在同时进行的EEG-fMRI研究中实时减少梯度和脉冲伪影，该工具速度快，适用于任何EEG-fMRI设置，并且作为NeuXus的一部分公开提供，NeuXus是用Python开发的实时EEG处理工具箱。该方法以原子吸收光谱(AAS)技术为基础，利用长短期记忆(LSTM)人工神经网络检测心电R峰。我们首先简要介绍了NeuXus工具箱，然后描述了实时伪影减少方法，然后对其性能进行了回顾性评估(基于先前收集的数据)，并与成熟的离线方法(EEGLAB的AAS)和商业上可用的在线方法进行了比较。
2. 材料和方法
2.1 NeuXus脑电处理工具包
NeuXus是一个用Python开发的开源工具箱，允许用户通过创建脚本，实例化一组类。每个类实例最多有一个输入和一个输出变量(称为ports)和一个方法(称为update)，该方法以特定的方式将数据从输入转换到输出端口(例如，删除每个通道的梯度工件)。端口内的数据采用表(称为块)的形式，行表示时间，列表示通道。数据块可以包含信号(每个通道一个)或标记(在这种情况下只有一个通道和时间，数据值是标记的名称)。通过创建节点并将一些节点的输出端口设置为其他节点的输入端口，用户可以设计一个管道并执行它，使每个节点在循环中依次更新，并实时处理数据。
2.2 实时伪影减少
脑电图由采集软件采集并传输到NeuXus。在特定情况下，采集软件是BrainVision Recorder，数据由RDA以具有固定数量通道但可变数量时间点的块流式传输，并由NeuXus作为块接收。然后执行工件减少，如图1中的流程图所示。NeuXus伪影算法需要定义的主要参数如表1所示，其中一些也显示在图1中。在一台配备AMD Ryzen 7 3700 U处理器、Radeon Vega Mobile Gfx 2.30 Hz集成显卡和8GB内存的惠普笔记本电脑上运行了这次降低。由于LSTM估计步骤需要密集的计算，因此使用Numba编译器来加速PA缩减。
图1同时记录EEG-fMRI的实时伪影减少算法流程图。当前块(蓝色)与以前的块(灰色)平均到一个遗传算法模板中，该模板用于从当前块中减去遗传算法。结果被下采样，心电图过滤，并连接到一个检测窗口。检测R峰并将其用作触发器，将整个心动周期的EEG平均到PA模板中。最后，这是从检测窗口的最后一部分减去(直到边际)。GA将不断发送块，这些块将被推送到检测窗口，并且每推送一定数量(步长)，它们将触发一个新的检测和随后的PA模板平均和减法。
表1再伪迹去除算法中使用到的参数设置
2.2.1 梯度伪迹
对这些片段进行平均，以取消与TR(生理信号和除GA之外的伪影)无关的信号分量，并创建平均GA模板。在对最小数量的段进行平均之后(表1)，选择当前块所贡献的模板部分并从块中减去，删除其GA，并抛出一个标记来标记减法的开始。在对最大数量的段求平均值之后(表1)，当一个新段进入平均值时，最早的段将被删除。通过这种方式，模板是基于最新的TRs构建的，这使得它能够适应遗传算法的变化(例如由于头部运动)。
2.2.2 脉冲伪迹
经过GA减少后，块被下采样到250 Hz。这一步包括一个切比雪夫I型和8阶的抗混叠低通滤波器，这样可以防止新奈奎斯特频率(125 Hz)以上的频率混叠和污染感兴趣的频带。使用4阶巴特沃斯滤波器对心电信号进行0.5至30 Hz的带通滤波，以去除噪声并改善随后的R峰检测。然后，当它从GA减法开始接收标记时，通过将GA减法，下采样和ECG过滤的块连接到检测窗口中，开始PA还原。随着窗口的填充，每个连接的块被选中并输出，而不改变(图1:最后一次选择用灰色括号表示)，直到检测窗口中的特定位置(表1:距离结束的一个步幅和一个空白)，之后块被保留，直到窗口完成。
2.2.2.1 R峰值检测
当窗口完成时，通过检测相应心电信号上的R峰来发现心电周期。检测算法基于Laitala，并适应实时执行。LSTM网络由两个双向LSTM层和一个密集层组成。每一层将输入与可训练参数结合起来，并将其传递给激活函数。在LSTM层中，输入向前和向后扫描，激活函数是双曲正切。在密集层中，激活函数是一个sigmoid，输出是一个与网络输入大小相同的数组。检测首先对-1和1之间的心电信号进行归一化，并将其传递给LSTM网络，该网络估计每个时间点出现R峰值的概率。然后将该概率与之前所属的所有窗口在该点估计的概率平均。太靠近窗口末端的点的概率，因为假设网络没有足够的点来做出明智的决定。这些概率被设定了阈值，以产生过多的正R峰值标签。然后标签经过两个过滤步骤:(1)只选择被最小数量的正标签包围的正标签; (2)如果在一个心动周期内发现多个阳性标记，则只保留概率最高的一个。最后，利用滤波后的正R峰标签，将检测窗口的EEG信号分割为心动周期段，选择之前持有的窗口段进行相减(图1:蓝色括号，直到边缘)。
2.2.2.2 平均工件减法
每次检测完成后，取心动周期段平均值，生成PA模板。为了适应心率的变化，为模板设置了最大心循环持续时间(表1)，这样它可能包括周期的每个部分的不同数量的片段。当选择检测窗口段进行减法时，选择模板中对应的部分，该部分必须使用最小心动周期数 (图1:暗橙色部分)进行平均，并从该段中减去，去除其PA。然后输出PA减法段，并抛出一个标记来标记PA减法的开始。如果模板中相应的部分有一个部分尚未使用最小循环数进行平均，则该部分将不会被减去，并且输出部分将被部分pa减少。事实上，当模板不够精确时，最好不要减去PA。当遇到异常长的周期时，可能会发生这种情况，可能是自然的，也可能是因为检测器错过了R峰值(在这种情况下，将值不正确地平均到模板的较远部分，从而阻止了减法)。尽管如此，这些情况并不常见，并且可以通过减少最小平均值数量或增加最大平均值数量来最小化。进入检测窗口的新块将移动那里的块并移除最早的块。每次发生一定的位移(表1)，就将LSTM检测器应用于心电，更新心动周期段，并选择下一个窗口段。因此，只有在检测时才会输出片段，这确保了它已正确地进入模板并被它减去。图2显示了实时GA和PA减少的示例。NeuXus还允许使用实验室流层(LSL)发送干净的数据，在查看器中进行监控(例如在StreamViewer中)或保存在LabRecorder中。
图2 同时记录的EEG- fMRI的实时伪影减少的例子。
2.3 数据采集
用于训练NeuXus LSTM网络和评估NeuXus假影减少算法的数据是在3T Vida MRI扫描仪(德国西门子公司)上获取的，使用核磁共振兼容的脑电图系统，该系统具有32通道BrainCap、SyncBox和Triggerbox。根据《赫尔辛基宣言》，数据采集得到委员会的批准，所有受试者都提供了书面知情同意。如有合理要求，本处将提供所有资料。验证数据可从https://doi.org/10.6084/m9.figshare.22561555.v2获取。按10-20制将帽套于头部，背部放置1个ECG电极。使用BrainVision Recorder以5000hz的采样率获取数据。fMRI数据采集采用64通道头部射频线圈，采用2D-EPI序列采集60张全脑轴向切片，各向同性分辨率2.2 mm。
2.3.1 用于训练LSTM网络进行R峰值检测的数据
NeuXus LSTM网络使用从6名女性志愿者收集的心电图数据进行训练，每个受试者在执行一系列认知任务时同时进行EEG-fMRI扫描，共32分钟。该网络在10个epoch中进行训练，从该数据集的ECG中随机采集40批256个窗口，跨度为2s(检测窗口也设置为2s)。
2.3.2 用于评估GA和PA减少的数据
对三组健康受试者的数据进行回顾性评价，每组受试者在不同的条件下:静息状态(RS)、睁眼/闭眼(EO/EC)任务和运动意象(MI)任务。受试者数量、运行次数和每次运行的持续时间见表2。在RS条件下，受试者被要求盯着黑屏上的白色十字2分钟。对于每个受试者，在MRI扫描仪内进行第二次RS EEG采集，但不同时进行fMRI，仅捕获PA而不捕获GA，并评估PA校正。在EO/EC任务中，受试者被要求在黑色屏幕上注视一个白色十字架2分钟，然后再闭上眼睛2分钟。在MI任务中，受试者观看一个屏幕，显示一个人在两排船上的第一人称视角(这种场景称为NeuRow)。范式如图3所示:在每次试验中(每次18次试验)，其中一只手臂开始划动，受试者被指示想象他们是移动手臂的人。
表2用于评估GA和PA算法的数据集
图3运动意象范式的示意图与相应的NeuRow显示
2.4 评估
2.4.1 R峰检测
NeuXus的LTSM R峰值检测算法在EO/ EC数据集上进行评估，使用F1分数与地面真实值进行比较(手动峰值检测)。将其性能与原始LTSM离线算法和常用的Pan-Tompkins算法进行比较。采用Kruskal-Wallis检验确定这些方法之间的统计学显著差异。当零假设被拒绝时(在p<0.05的显著性水平下)，使用Dunn检验进行方法之间的两两比较。
2.4.2 伪迹降低
将提出的实时遗传和PA算法与唯一的实时工具(RecView)以及基于相同算法的成熟且常用的离线方法(EEGLAB AAS校正)进行了比较。在RecView中，使用MRI伪影和脉冲伪影滤波器来减少伪影。在MRI滤波器中，TR设置为1260 ms，在PA滤波器中，平均脉冲数设置为30。在EEGLAB中，使用fMRIb插件。
在GA减少中，GA平均值的TR和窗口数分别设置为1260 ms和20，并且禁用了执行OBS和ANC的选项。在PA约简中，PA平均的窗口数设置为30，约简方法设置为“mean”，只执行AAS。利用RS数据分别对GA和PA约简阶段进行评估，方法是比较每个约简阶段前后伪信号对应频带(伪信号频带)与无伪信号频带(背景频带)的EEG频谱功率。利用MATLAB实现的快速傅里叶变换(FFT)计算频谱功率。对于遗传算法，伪带被定义为0.04 Hz宽的间隔，以1/TR (1/1.26 Hz)为中心，其谐波高达30 Hz，而背景带被定义为连续伪带之间间隔的50%(图4-顶部)。对于PA，在PA峰周围定义伪带，通过在平均心率及其谐波附近拟合洛伦兹量来识别，而背景带的定义方式与GA相同(图4-底部)。
使用EO/EC和MI数据集评估总体伪影减少(包括GA和PA)，通过分析其对典型脑电图闭眼反应性(枕叶皮层α功率增加)和运动皮层激活(甚至在大脑皮层)的影响。
对于MI数据集，对于每个受试者，每次运行的数据都被单独校正，然后通过从1到40 Hz的带通滤波，插值坏通道(使用EEGLA工具)并参考通道的平均值进行串联和预处理。对左臂试验进行分析，选择通道C4(因为它与左初级运动皮层相对应，因此通常用于测量与身体右侧运动相关的脑电图活动)，并计算α - ban中的事件相关去同步(ERD)(因为这通常用作运动激活的脑电图测量)。为此，事件相关谱摄动(ERSP)的计算方法如下，利用短时傅立叶变换(STFT)将每次试验中的信号在时频域转换为功率信号，对各试验中的功率信号进行平均，并将结果按基线周期归一化。
图4 在代表性示例的频谱上，GA(顶部)和PA(底部)的伪影(绿色)和背景(红色)带的说明。在伪影还原前和还原后(NeuXus, RecView, EEGLAB)的光谱上显示了叠加的条带。下面显示了最上面一行的放大部分，以澄清不同方法之间频谱的差异。
2.4.3 实时性
为了评估该方法的实时执行能力，对NeuXus中每个数据点在实时伪影减少的每个阶段所花费的时间进行了一次说明性RS运行。在每个阶段(GA降采样、下采样、心电滤波和PA降采样)之前和之后，每个块的数据点都用当前时间(以时间度量)进行时间戳。在实践中，只有块行(时间实例)有时间戳，因为对于每个时间实例，列(通道)的数量是恒定的，并且假定这些值是并发的。为管道中的每个时刻保存时间实例和时间戳。在NeuXus执行之后，匹配每个阶段前后的时间实例，并减去相应的时间戳，以获得该阶段每个数据点所花费的时间。然后用这些时间来计算中位数、25%和75%的百分位数。通过跟踪数据点(通过它们的时间实例)而不是完整的块，每个块的可变长度不会影响计时，并且在PA缩减中保留的点(直到检测)只在输出时接收它们的最终时间戳，从而确保正确测量它们在该阶段的时间。
3. 结果
3.1 R峰检测
R峰检测结果如图5所示。NeuXus的LSTM网络分类F1评分高于PanTompkins。
3.2 伪迹去除
RS数据集中GA约简的评估结果如图6所示。在伪影带中，不同方法间无显著差异。在背景波段，NeuXus和RecView产生的功率降低比离线校正更小(分别为3.1%和8%)。总的来说，除背景波段的RecView外，不同频道(受试者平均)的减少遵循类似的模式;在这种情况下，在CP2附近发现了更大的下降，可能是由于一个受试者大大偏离了中位数(通道平均值的箱线图也很明显)。
RS数据集中PA约简的评估结果如图7所示。在伪影带中，NeuXus和EEGLAB均比RecView产生更强的功率降低(分别为64和58比39%，跨通道和受试者的平均值，p< 0.05)。在背景波段上，两种方法无显著差异。所有的减少在不同的频道和对象中都遵循相似的模式。
图9显示了EO/EC数据集的结果，显示NeuXus和其他方法对闭眼的alpha功率反应性没有显著差异，只有EEGLAB和RecView之间(136比61%，受试者平均，p< 0.05)。
MI数据集的结果如图10所示。两种方法的ERD差异无统计学意义。计算时序结果如表3所示。使用Numba时，NeuXus的中位执行时间为0.2304s，不使用Numba时为0.6407 s (Numba使执行时间减少64%)。在Numba中，GA约简、下采样、心电滤波和PA约简的中位执行时间分别为0.0013、0.0045、0.0020和0.0020分别为0.2226秒。
图5 R峰检测评价:左) 5 s窗口R峰检测示例，NeuXus漏了一个峰，Laitala和Pan-Tompkins漏了两个峰。右) NeuXus, Laitala和Pan-Tompkins方法在所有科目中的F1分数(每种颜色代表一个科目)。
图6GA还原的评价。每种方法(NeuXus, RecView, EEGLAB)的伪影带(上)和背景带(下)的功率变化:左)平均受试者的地形(显示跨通道的分布); 右)箱形图(显示受试者之间的分布)显示了通道之间的平均值。不同的颜色代表不同的被试。
图7PA还原的评价。每种方法(NeuXus, RecView, EEGLAB)的功率变化随伪影带(上)和背景带(下)的PA减小而变化:左)右)箱形图(显示受试者之间的分布)显示了通道之间的平均值。不同的颜色代表不同的被试。
图8伪影减少对静息状态数据(RS数据集)的说明性时间段的影响。
图9相对于睁眼条件，人工干扰减少对闭眼α功率降低的影响。左)一个说明性受试者在EO(蓝色)和EC(红色)条件下的功率频谱;单个alpha波段用灰色表示。右)所有受试者的alpha功率降低分布(跨通道平均)。
图10伪影减少对运动想象过程中α波段事件相关去同步(ERD)的影响(MI数据集)。左)典型受试者左MI实验C4通道ERSP图，使用三种伪影约简方法。右)ERD在所有受试者中的分布，即ERSP在alpha频带(8-12 Hz)和任务时间间隔(0.5 -4 s)上的平均值。不同的颜色代表不同的受试者。
4. 讨论
这项工作提出并验证了NeuXus开源工具，用于使用传统硬件设置实时减少同时获得EEG-fMRI的梯度和脉冲伪影。我们发现，NeuXus在静息状态下具有背景信号保存的伪影功率降低以及运动图像期间闭眼和去同步对脑电图α功率反应性的影响方面的表现至少与离线还原(EEGLAB)和唯一的商业可用在线工具(RecView)一样好。
其次，NeuXus在构建模板之前不会对数据应用高通滤波器。EEGLAB过滤基线振荡(<1 Hz)以改善伪迹估计(因为它们不是遗传算法所期望的)，然后将模板减去未过滤的数据，保留它们，这取决于实验，可能对分析有用。在实时GA约简中没有实现过滤器，因为它在数据上引入了延迟，这会使模板与未过滤的数据不匹配，从而影响减法。在脱机场景中，具有一半订单的过滤器可以应用两次，一次向前，另一次向后，以取消延迟。然而，在实时情况下，完整的记录是不可用的，并且在单个块上应用这种策略是不可行的，因为它会在边缘扭曲每个块。一种可能的解决方案是使用线性相位的高通滤波器，它对所有频率分量施加恒定的延迟(群延迟)(避免相位失真)，然后将滤波后的信号移回该延迟。线性相位可以用FIR滤波器实现，也可以用IIR滤波器近似实现(计算效率更高)。由于这将增加算法复杂性，结果接近最佳性能方法，并且认为低于1 Hz的基线不具有广泛相关性，因此未实现。值得注意的是，可以在遗传算法约简本身之前应用过滤器，从成为模板一部分的数据中删除基线，但也会被它减去基线。基线不会被保留，并且每个点都有一个小延迟，但是减法不会不匹配，因为过滤后的数据将同时用于模板和减法。
第三，NeuXus不会将模板的振幅缩放到块上。EEGLAB将模板乘以一个常数，以最小化模板和数据之间的最小二乘。虽然进行了测试，但没有发现它会影响结果。RecView也应用TR对齐步骤，通过构建用户定义的模板数量，只从紧密对齐的块中选择最接近当前块的模板进行减法。经遗传算法降阶后，对数据进行下采样，并对心电信号进行带通滤波。滤波器引入了一个延迟，使心电图相对于脑电图不一致。对于8-50 Hz的QRS频率范围，估计平均延迟为0.0012 s。这意味着PA可能在稍微移位的R峰处被分割。与遗传算法一样，对每个块向前和向后应用带通滤波器会扭曲每个块的边缘。作为一种替代方法，它可以应用于R峰检测中串联的心电信号(表1:filter_ecg)。在这个阶段，ECG不仅更长(因此比组成它的块具有更少的边缘)，而且每个边缘对伪影减少的贡献也很低(第一个与当前块不一致，第二个位于边缘内)。然而，PA大多包含较低的频率，并且无论如何都不是由R峰精确计时的。因此，我们认为轻微的延迟以及由此产生的与EEG的不匹配不会影响伪影的减除，保留PA减除前的滤波器。
5. 总结
这项工作提出了一种实时GA和PA约简方法，用于fMRI期间记录的EEG数据，该方法是完全开源和公开可用的，快速且与硬件无关，提供了第一个结合所有这些特性的工具。该方法的性能至少与基于相同算法和唯一可用的与硬件无关的实时缩减(专有)软件的离线校正一样好。它提供了一个可以配置的必要的短延迟。干净的数据可以用来监测或进行神经反馈训练。该方法是开源Python平台NeuXus的一部分，可以很容易地使用和扩展。
参考文献：NeuXus open-source tool for real-time artifact reduction in simultaneous EEG-fMRI.
