2️⃣使用独立成分分析(ICA)清理伪影

type
status
date
slug
summary
tags
category
icon
password

引言

在本教程中,我们将演示如何在ArtifactMEG.zip示例MEG数据集上进行ICA清理。这是一个没有明确实验任务的静息态记录,但在记录过程中,受试者特意制造了某些类型的伪影,比如眨眼、来回看、做头部运动,以及咬牙以在颌肌中产生肌电活动。
该数据集不仅包含MEG数据,还包含来自EOG和ECG的独立双极记录,以及记录颈部和颌部肌肉活动的两个双极EMG通道。原则上,我们可以使用这些额外的通道来查找伪影,或促进甚至自动识别ICA伪影成分。在接下来的内容中,我们将使用EOG、ECG和EMG通道,我们只展示如何在MEG通道上使用ICA。

背景

独立成分分析(ICA)是一种时空分解策略,它假设EEG或MEG的潜在信源具有固定的空间投射到各个通道,并且在时间上最大程度地独立。使用ICA可以估计与通道数量相等的成分数。我们知道神经元信源的数量远多于通道数量,特别是在记录时间很长的情况下。因此,ICA分解只能对最显著的独立成分进行近似。我们也知道,数据中显著的伪影贡献通常相当有限,这意味着我们可能只需要少数几个成分就能解释这些伪影。在ICA分解之后,我们可以识别出伪影成分,并将其他所有成分反投影到通道层面,从而排除伪影。
由于ICA假设信源在空间上是固定的,且只能估计有限的信源集,我们应该尽量预处理数据,以避免ICA成分被简单的数据伪影"占用"。如果你的受试者在实验尚未开始的前几分钟仍在移动,你会希望将该部分数据从ICA分解中排除。
出于效率考虑,我们通常通过首先识别感兴趣的试次并仅将这些数据读入内存来进行预处理。然而,某些典型的伪影可能在试次间隔中更为频繁,例如当受试者不需要保持注视时更可能眨眼或做眼跳。因此,包含这些试次间隔可以有助于识别与眼动相关的成分。

程序

要使用ICA清理MEG数据,我们将遵循以下程序
notion image
在流程示意图中,您可以看到我们将使用ft_rejectvisualft_databrowser来识别非典型伪影。数据浏览器将帮助您深入理解数据的特征,这在您不是亲自记录数据的情况下尤其重要。可视化概要模式可以非常快速且高效,这对于具有大量通道或较长时间的数据集特别重要。
此外,流程示意图显示,如果我们在完成ICA之后发现非典型伪影或坏通道,我们可能需要返回到ft_rejectvisualft_databrowser。在这种情况下,我们需要重复初始清理并重新进行ICA。

预处理

MEG数据是使用151通道CTF系统记录的。虽然数据是连续的,但CTF数据集被组织成10秒的试次。由于试次之间没有不连续性,我们可以将其视为连续记录。
要将此数据集加载到MATLAB中并使用FieldTrip进行预处理,请使用:
ICA分解可能需要很长时间,特别是当您需要多次执行时。如果耗时太长,您可以考虑将数据降采样到较低的采样率,这里从1200Hz降至300Hz。使用降采样后的数据集,您可以估计ICA成分,然后从原始数据中移除这些成分。但是,您应该意识到降采样可能会影响某些数据特征和伪影。例如,高频肌肉伪影的捕获效果会降低。这导致这些伪影的去除可能也不那么理想。
出于本教程的目的,并且由于我们想让您尝试各种设置,我们将继续使用降采样后的数据。如果您在处理自己的数据时,请尽量使用原始数据而不是降采样数据;这将提高分解的质量。

去除非典型伪影

ICA假设有固定成分的混合,且无法估计超过通道数量的成分。如果您有一些不频繁的非典型伪影,这些将在成分中表现出来。这会导致损失一些用于有趣内容的成分,并可能导致次优的分解。因此,我们首先要去除稀疏的非典型伪影

使用ft_databrowser

我们可以使用ft_databrowser来查看数据。由于我们在寻找可能出现在数据任何位置的不频繁伪影,我们应该查看所有通道和完整的时间过程。在时间上"缩小"视图会有帮助,这样我们可以一次看到较大的数据时间窗口。
注意,我们对ft_databrowser的输出cfg感兴趣:它将包含我们标记的伪影。
notion image
在记录的最后部分出现了一些异常现象:CTF采集软件以10秒为块写入数据。如果在最后一个数据块完成之前结束采集(这种情况总是会发生),该块的剩余部分将被写入磁盘为全零值。因此,您应该将数据的最后部分标记为伪影。
在MLT024通道(以及同时在其他一些通道)中,您可以观察到所谓的SQUID跳变。这些是由SQUID硬件不稳定导致的,造成信号中出现一个磁通量子的巨大变化。这些表现为大幅度跳变。请识别MLT024中的每个跳变并将其标记为伪影。
作为参考,以下是我在降采样数据中识别的伪影:
以下是这些伪影在连续记录中出现的时间(以秒为单位):
请注意,如果在初始预处理时没有应用高通滤波器,并且我们没有对数据进行降采样,这些跳变会更容易识别。特别是高通滤波器在这里比较麻烦:它使跳变扩散到更长的时间范围。

练习

预处理和滤波用于减少伪影(如漂移),因此会使伪影不太明显。未经预处理的数据能最好地展示伪影的特征。
在原始数据上重复使用ft_databrowser进行检查,不要使用高通滤波器和重采样。这需要您重新调用ft_preprocessing。查看记录开始后约178秒处以MTL*开头的通道,并将跳变与您在经过滤波和降采样的数据中观察到的跳变进行比较。
请注意,伪影(和试次)是以采样点表示的,所以您可以查看记录中的相同时间点,但不能使用相同的采样点数来引用原始数据和降采样数据中的伪影。
在识别非典型伪影之后,您可以将它们从后续分析中移除。出于本教程的目的并且是在处理连续数据,我们不会将其剔除,而是用NaN值(非数字值)填充。在连续数据上,也可以使用ft_rejectartifact,配合选项cfg.artfctdef.reject = 'partial'。这里我们将使用选项nan
注意,这里我们以不同方式重复使用datadata_clean变量,因为我们在来回演示不同的数据处理和清理方法。不要对清理后的数据代表什么感到困惑。
要回到原始数据,您随时可以执行data = data_orig

使用ft_rejectvisual

移除不频繁和非典型伪影的另一种策略是使用ft_rejectvisual。但是,这需要将数据分段为试次。在这种情况下,我们可以将数据分段为一秒长的连续片段。
注意,这里我们依赖于预处理期间的cfg.hpfilter='yes'选项。信号中的低频漂移并不重要,而且会影响ICA分解;因此我们无论如何都要去除它。由于我们处理的是连续数据,我们使用高通滤波器。如果我们最初就是从基于试次的数据开始,那么在最初预处理阶段使用cfg.demean='yes'会是一个更简单的选择。
现在数据被切分成每秒一段的片段(即试次),我们可以识别具有非典型伪影的片段。为此,我们使用带有summary方法的ft_rejectvisual。这将直接返回清理后的数据。由于我们不想从后续分析中排除任何通道,我们指定要保留所有通道。由于我们希望数据之后能够再次连续表示,我们指定用NaN替换伪影。
notion image
现在我们可以将分段数据重新"缝合"成连续表示:
如果您在ft_databrowser中检查data_clean,您会看到数据的某些部分不可见,这些部分已被NaN值替代。

ICA分解

我们使用ft_componentanalysis进行ICA分解。它有许多选项,支持不同的数据分解方法,包括PCA和不同的ICA算法。这里我们将使用来自EEGLAB的Extended Infomax算法,即runica方法。您不需要安装EEGLAB,所需的函数都包含在fieldtrip/external目录中。
您可以考虑其他ICA算法。例如,fastica算法速度快,不需要对数据进行完整分解;它也可以只识别少量成分。fastica首先识别的成分是具有最大方差的成分,这些通常是伪影。amica算法是识别生物物理学上合理的ICA成分的最佳算法之一。它产生的成分比runica更加独立。
Delorme等人的论文独立EEG源是偶极子的。PLoS One (2012)比较了多种ICA算法,并得出结论:更独立的成分也更具偶极子特性,这与许多最大独立EEG成分是单个紧凑皮层区域内部分同步局部皮层场活动的体积传导投影的解释相符。
另一种不完全盲目的成分分析策略是"去噪源分离"或dss,这种方法在ICA分解过程中最大化特定的显式特征。如果您预先知道EOG或ECG伪影的样子,您可以使用dss来利用这一点并高效地清理数据。
要执行ICA分解,您可以使用以下代码。这里我们在分解之前使用PCA降维来加快处理速度,以适应本教程的目的。请注意,PCA降维可能会对分解质量产生负面影响,所以通常最好在计算运行时保持更多的耐心。
您应该采用类似的方式,这些共同代表了底层源的线性混合。如果您有使用相同参考的EEG和EOG通道,EEG和EOG通道可以一起使用。如果您的EEG使用一个参考而EOG通道是双极性的,则不应将它们组合在单个ICA分解中。当您同时有MEG和EEG数据时,MEG和EEG通道都会检测到大脑和伪影源,但对它们的敏感度不同,对头部相对于MEG头盔运动的效应也有不同的敏感度;在这种情况下,我们也建议分别对MEG和EEG通道使用ICA。

识别伪影成分

分解后的数据结构data_comp表示成分的地形图,即每个源如何投射到通道,以及表示每个源的时间激活。这两者都可用于识别对应于伪影的成分。
notion image
有了一些经验,您就可以相对快速地识别可疑的伪影成分。眼动相关的成分在前额通道上呈现局部空间分布,眨眼和垂直眼跳呈对称分布,而水平眼跳则呈现明显的左右模式。MEG中与心跳相关的成分表现为一个非常深的信源,在头盔的左右两侧呈现双极性投射。对于眼动和心跳成分来说,通常都会观察到多个相关成分。您可能需要记下每个可疑成分的编号。
在上图中,成分10和12看起来是眼动相关的成分。成分5和6可能仍与MLT24及其相邻通道的SQUID跳变有关。在这种情况下,没有明显的心跳相关成分。
接下来您可以查看成分的时间进程。心跳相关的成分会显示规律的心跳,在分解的最初几秒内就应该能够识别出来。眼动相关的成分只有在受试者眨眼或做眼跳时才会在时间序列中出现偏移;要看到这些,您可能需要滚动浏览数据或缩小查看更大段的数据。
notion image
再次记下代表伪影的成分。
请注意,MEG中的眼动成分看起来与EEG中的眼动成分不同,这是由于MEG通道的方向所致:头部左侧的通道指向右侧,右侧的通道指向左侧,这导致它们的"极性"发生翻转。而对于EEG电极,所有额叶电极具有相同的极性。结合成分时间序列,您可以确定哪个成分反映水平运动,哪个反映眨眼。

识别不良通道

EEG和MEG都对生理活动给出相对模糊的表示。在空间上高度局部化的成分,即仅在单个或极少数通道上活跃的成分,不太可能代表来自大脑或心脏的生理源。如果您之前没有移除SQUID跳变,这些会表现为高度局部化的伪影。移动的EEG电极或阻抗突然变化的电极也可能表现为高度局部化。
如果发生这种情况,您可以使用ft_databrowser查看成分时间序列,以确定跳变或其他短暂伪影出现的时间,您可以将该时间窗标记为视觉伪影,并使用ft_databrowser的输出cfg返回并在ICA分解之前从数据中删除该部分(或像我们在这里做的那样用NaN填充该部分)。随后,您需要重新进行ICA分解,并再次检查成分地形图和时间序列。
也可能会发现一个或几个成分在空间上高度局部化,这是由于相应的通道在较长时间段内状态不佳或有噪声。在这种情况下,您需要从数据中剔除这些通道,并在没有这些通道的情况下进行ICA分解和反投影。如果许多受试者都有不良通道,且/或这些不良通道在受试者之间变化很大,您可能需要使用ft_channelrepair对这些通道进行插值;这是在成分反投影之后获得清洁的通道级数据时要做的事情,但要在进行时锁ERP或频率/时频分析之前完成。
使用高密度EEG时,您有时可能会看到非常局部化的肌肉抽搐,特别是在颞区,但也可能出现在头皮的其他部位。这些可能较少或较多发生,这取决于您的受试者和任务。尽管这些在空间上相当紧凑,但它们确实代表了一个生理源,ICA是去除它们的适当技术。

移除伪影成分

在识别出伪影成分后,您可以使用ft_rejectcomponent将所有成分反投影到数据的通道级表示,同时排除伪影。
如果您在数据重采样版本上计算了成分,您也可以使用ft_rejectcomponent将伪影从原始数据中投影出去。这只需要成分地形图,然后将其应用于以原始采样率解混和重混数据。

总结和结论

在本教程中,我们研究了如何准备数据进行ICA分解,如何处理不频繁和非典型伪影。我们简要讨论了不同的ICA方法,并使用特定方法展示了如何通过降采样数据和减少估计的成分数量来加快ICA速度。给定ICA分解后,我们演示了如何可视化它们并识别伪影,并讨论了如果得到非常局部的成分(在空间和/或时间上),您可能需要删除数据的该部分并重新开始。在获得干净的ICA分解后,我们解释了如何将成分反投影到通道级表示以进行进一步分析,或如何使用成分来清理原始数据。
在本教程中我们没有展示的是,您还可以将ICA成分解释为大脑中的源。正如可以对伪影成分使用所有FieldTrip绘图和分析策略一样,也可以分析大脑成分。所有成分都表示为通道时间序列,可以应用所有通道级分析。例如,在将数据分段为刺激锁定试次后,您可以使用ft_timelockanalysisft_freqanalysisft_componentanalysis来研究每个感兴趣的成分,并使用统计方法比较条件之间的差异。
分解不仅给出了成分时间序列,还给出了它们的地形图。除了用这些来识别伪影外,您还可以清楚地识别代表大脑活动的成分。您可以使用ft_dipolefitting用简单的偶极子模型定位这些成分,或使用ft_sourceanalysis用分布式模型定位它们。您不能(轻易)应用的唯一源重建策略是波束成形,因为在成分级别您只有一个地形图,没有数据协方差矩阵。作为空间滤波技术,波束成形与独立成分分析太相似(尽管基于生物物理模型和非相关而不是独立源的假设)。

建议进一步阅读

关于如何在FieldTrip中一般处理伪影的介绍,您应该看看伪影处理简介教程。在本教程使用ICA清理数据之后,您可能想回去看看视觉伪影拒绝自动伪影拒绝教程。
 
上一篇
自动伪影拒绝
下一篇
参考文档
Loading...