教程 | 近红外数据的预处理和平均（上）

前言

近红外光谱(NIRS)是一种测量流经传感器所在组织的血液中氧合水平的方法。它基于这样一个事实，即含氧血红蛋白和脱氧血红蛋白具有不同的吸收光谱，因此你会看到它有不同的颜色。大多数近红外系统在每个光源光电二极管发射2个波长的光，通常约700和850nm，有些系统则使用3个波长。

光电探测器被放置在离源光极几厘米远的地方。光源和探测器的光电二极管共同组成一个通道，但系统对于一个通道总是有两个空间位置(大多数人认为记录的活动主要来自光源和探测器的中间位置)和两个波长。光散射穿过组织，一小部分光将通过光电探测器在组织的表面被探测到。在头皮上检测到的光强随着位于源和探测器光电二极管之间的组织中含氧和脱氧血红蛋白浓度的变化而有系统地波动。由于(或多或少的随机)散射，大多数到达探测器的光只穿过了组织表面，只有一小部分穿过了表面下更深一点的区域。通常情况下，在头皮上检测到的光穿过了一个组织体积，从侧面看是一个香蕉形状。如果光源和探测器被放置得更远，那么被探测到的光可能在组织中传播得更深。但实际上，这是有限的，因为到达探测器的那部分光会随着光源-探测器距离的增加而急剧下降。当测量大脑中血红蛋白的变化时，颅骨作为一个重要的屏障，吸收了很多光。因此，对大脑深处进行扫描只有在头骨较薄的婴儿中才可行。

氧和脱氧血红蛋白的浓度不能以绝对条件评估，这意味着所有的测量都将反映的是相对变化。因此，像EEG和MEG一样，需要一个基线期。在血流动力学响应函数(HRF)上，记录和后期处理的fNIRS信号与fMRI类似，但具有更高的时间分辨率。因此，至少在理论上，可以测量曲线的初始倾角。与fMRI相比，空间分辨率要差很多，通常我们主要从大脑表面进行记录。因此，绘图通常类似于EEG。EEG的时间分辨率更高，但fNIRS的一个优势是，测量的信号只来自光源和探测器之间的区域。fNIRS是一种非常有前景的方法，可以测量无法通过fMRI扫描仪测试的人群的大脑活动，并且允许被试做出无法在扫描仪中完成的更多动作。但是血压会随着你所做的任务而变化，从而也会影响测量到的信号。

单通道近红外数据的预处理和平均

本文将演示如何分析单个通道的功能性近红外光谱(fNIRS)数据。主要使用FieldTrip对Artinis近红外光谱数据进行基本分析。由于本文分析的是单通道数据，因而没有展示如何处理坏通道，关于如何删除坏通道和分析多通道的信息详见近红外数据的预处理和平均(下)，关注茗创科技，不迷路哦~

本教程所使用的数据集

本文将使用到Artinis Medical Systems的Oxymon MK III系统采集运动皮层信息的近红外数据集。光电极放置在运动皮层上，随后要求被试进行手指敲击。除了来自大脑的记录，还记录了一个带有标注的事件(event)通道，在事件通道中，参与者开始敲手指时的记录为“A”，停止敲手指的记录为“B”。这个运动任务一共重复了12次，所有的数据，包括事件通道，都保存在一个.oxy3文件中。本文中使用的数据可以从https://download.fieldtriptoolbox.org/tutorial/nirs_multichannel/上获得。

分析步骤

根据数据和实验设计的不同，分析可以有许多不同的方式和顺序。这里将介绍分析步骤的标准顺序，随后小伙伴们可以逐步测试这些步骤。以下步骤为分析fNIRS数据提供了一个很好的标准方法，如图1所示:

1、读取连续数据；

2、(优先地)裁剪记录开始和结尾的数据，特别是当这些需要被剪切的数据是噪声或不感兴趣的片段时；

3、删除通道中的坏段(例如，运动伪影)；

4、将光密度转化为含氧血红蛋白(oxyHb)和脱氧血红蛋白(deoxyHb)浓度的变化；

5、将功能响应与系统反应分开，可以使用以下方法其中之一，或者组合使用；

①滤波

②减去参考通道

③每个通道都有反相关的oxyHb/deoxyHb

6、在实验任务中，定义与试次相对应的段；

7、将试次数据平均并可视化。

图1.标准fNIRS分析流程概览。

本文分析的数据集名为“motor_cortex.oxy3”。首先，需要将FieldTrip指向这个文件。注意，FieldTrip没有图形用户界面(GUI)，而是通过脚本来实现最终目标。需要将FieldTrip文件夹添加到MATLAB路径中，然后执行ft_defaults函数，输入：

ft_defaults

通常，FieldTrip会将函数的所有输入参数存储在一个名为“configuration”或“cfg”的本地变量中。虽然MATLAB变量可以取任何名称，但为了方便起见，这里将继续使用名称“cfg”。分析流程中的一个(关键)步骤是读入数据并对数据进行预处理。为此，这里将使用FieldTrip函数ft_preprocessing。因为FieldTrip是一个开源工具箱，可通过输入以下代码来查看代码细节：

edit ft_preprocessing

读取和裁剪数据

在变量cfg.dataset中输入数据的文件名。

cfg = [];cfg.dataset = 'motor_cortex.oxy3';

通常需要指定完整路径，包括驱动器号。这可以确保代码将独立于当前的MATLAB位置执行。将原始数据、处理过的数据和脚本保存在三个不同的位置是一个很棒的做法(推荐)。在本例中，将使用当前工作目录中的数据。

现在，尝试执行ft_preprocessing：

[data] = ft_preprocessing(cfg);

在处理.oxy3文件时，需要加载包含光极布局的光极模板。如果MATLAB无法在你的路径上找到模板文件，它将会弹出一个图形用户界面对话框，要求你定位XML文件。默认情况下，该文件位于C:\Program Files(x86)\Artinis Medical Systems BV\Oxysoft 3.0.103。需要选择的文件是optodetemplates.xml。鉴于FieldTrip经常会需要搜索此函数，最好将此函数复制到存放MATLAB分析脚本或fNIRS数据的文件夹中。关于如何选择最佳模板的信息，请参阅Artinis网页上的这篇博客文章(网址：https://www.artinis.com/blogpost-all/2017/6/27/how-do-i-choose-the-correct-fibers-and-template-for-my-oxymon)。

当你在屏幕上看到类似这样的内容时，FieldTrip就完成了。

>> the call to "ft_preprocessing" took 9 seconds

仔细观察workspace，就会注意到添加了一个名为“data”的新变量。这个名称与ft_preprocessing前面输入的名称相同。如果写入[something]，那么变量'something'就会被添加到workspace。

如果按照以上描述进行预处理，数据将具有以下字段：

data =        hdr: [1x1 struct]      label: {48x1 cell}       time: {[1x4462 double]}      trial: {[48x4462 double]}    fsample: 5 sampleinfo: [1 4462]       opto: [1x1 struct]        cfg: [1x1 struct]

“hdr”字段包含关于数据的所有最高级信息，例如原始采样率、通道数量等。字段“label”列出了你决定读入的所有通道的名称。字段“time”，表示数据集的时间轴。“trial”字段包含所有通道的数据，它被称为“试次”，因为通常情况下，数据被分割成与实验中的不同试次相对应的片段。字段“fsample”描述了“trial”字段中数据的采样率。“sampleinfo”字段描述了相对于磁盘上原始测量值的每个试次的样本编号。“opto”包含了有关通道和光电极组成的信息，例如使用了什么波长，光点极位置等等。最后，“cfg”字段与刚才使用的cfg相同，只是通过一些默认值进行了扩展。

为了快速查看数据，可以使用ft_databrowser函数。databrowser不仅仅是一个简单的“数据浏览器”，不过现在将利用这个功能来实现查看数据的目的。在配置中添加'ylim ='maxmin'来调整y轴，使图中的最小值决定y轴上的最低点，图中的最大值决定y轴上的最高点。

cfg = [];cfg.ylim = 'maxabs';ft_databrowser(cfg, data);

图2.在databrowser中显示原始数据。

使用ft_databrowser，还可以删除不需要的数据片段。例如，如果你在放置光电极的同时开始记录，可能会在记录的开始有一大段不需要分析的数据，这些数据包含与大脑无关的值，并且快速波动。把这些片段剪掉更有利于后续的分析。本示例数据不需要裁剪。

此外，这里先只选择一对通道，以减少分析过程的复杂性。

cfg = [];cfg.ylim = 'maxmin';cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'}; % you can also use wildcards like 'Rx4b-Tx5*'ft_databrowser(cfg, data);

图3.在databrowser中显示一对通道(两个波长)。

去除伪影

近红外光谱数据不仅包含了氧和脱氧血红蛋白浓度的变化，而且还包含了其他生理信号、随机噪声和来自测量环境的变化。因此，可能存在较为明显的运动伪影，它是由光电二极管和头皮之间接触的临时变化产生的，通常是由头部运动引起。这些运动伪影可能存在于所有通道中，但也可能只有一个或几个通道受到影响。比如，数据中短的、意想不到的峰值被认为由运动导致，可以通过ft_artifact_zvalue检测和删除这些伪影。

指定一些预处理选项来检测伪影，如下所示：

cfg = [];cfg.artfctdef.zvalue.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'};cfg.artfctdef.zvalue.cutoff = 5;cfg.artfctdef.zvalue.hpfilter = 'yes';cfg.artfctdef.zvalue.hpfreq = 0.1;cfg.artfctdef.zvalue.rectify = 'yes';cfg.artfctdef.zvalue.artpadding = 2;% cfg.artfctdef.zvalue.interactive = 'yes'; % the interactive display makes more sense after segmentating data in trials[cfg, artifact] = ft_artifact_zvalue(cfg, data);

运行上述代码可以发现，FieldTrip在此过程中识别了8个伪影。请注意，只是标识出来了这些伪影，还没有被删除。将其缓存，并在数据滤波分段之后调用ft_rejectartifact。

转化为oxyHB/deoxyHB浓度的变化

此时的数据是光密度(OD)值，而不是含氧和脱氧血红蛋白浓度。可以使用ft_nirs_transform_ODs将数据转换为浓度变化。在使用ft_nirs_transform_ODs时，要做出的一项选择是差分路径长度因子(DPF)，它取决于参与者的年龄和组织类型(例如，当分析的是肌肉组织而不是大脑的血氧变化时)。

cfg = [];cfg.dpf = 5.9;cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'};data_conc = ft_nirs_transform_ODs(cfg, data);

将功能响应与系统反应分开

除了出现在数据中的运动伪影，fNIRS测量还包含系统反应，这是研究中不感兴趣的，并希望从数据中能够过滤掉这些反应。有多种方法可以提取功能响应(大脑的血流动力学响应)，这里将采用一种简单但相对有效的方法，即应用带通滤波器来消除频带以外的信号值，该频带捕获了大脑血流动力学响应中的动态。由于实验任务和血流动力学响应函数，这里预计活动低于0.1Hz(对应~10s的波动)。此外，将较低的范围设置为0.01Hz(波动~100s)，以消除缓慢漂移。

cfg = [];cfg.bpfilter = 'yes';cfg.bpfreq = [0.01 0.1];data_filtered = ft_preprocessing(cfg, data_conc);

定义兴趣段

如果想知道大脑是否会对实验中发生的事件做出特定的反应，所以需要把分析的重点放在事件发生的时间窗上(这些时间窗通常被称为epoch或trials)。

此前，通过调用ft_preprocessing时不指定cfg.trl，读取了所有epoch。在帮助文档中，cfg.trl指向ft_definetrial，现在将使用它来定义试次。接下来先看一下函数：

help ft_definetrial

假设我们对数据中的triggers没有任何线索。因此，现在将利用隐藏的功能来检索这些信息。

cfg = [];cfg.dataset = 'motor_cortex.oxy3';cfg.trialdef.eventtype = '?';

问号表示不确定事件triggers，因此函数ft_trialfun_general将输出在数据集中找到的所有事件。接下来，调用ft_definetrial

ft_definetrial(cfg);

在命令窗口中，将看到发现了24个事件，并且使用了两种类型的事件，即事件'A'和事件'B'。此外，需要注意的是，现在没有定义输出变量'cfg'，因此变量'cfg'将保持不变。如果我们想看看从刺激前10s(事件'A'之前)开始到'A'之后35s(刺激后)结束的epoch/trials。可以定义试次，然后使用它来调用ft_preprocessing：

cfg.trialdef.eventtype = 'event';cfg.trialdef.eventvalue = 'A';cfg.trialdef.prestim = 10;cfg.trialdef.poststim = 35;cfg = ft_definetrial(cfg);cfg.channel = {'Rx4b-Tx5 [860nm]', 'Rx4b-Tx5 [764nm]'};data_epoch = ft_redefinetrial(cfg, data_filtered);

现在已经选择了一对通道，并在12个试次中剪切了数据。使用databrowser进行检查，使用以下这些设置，可以使绘制图形看起来更整洁，同时也可视化之前被识别出的伪影:

cfg = [];cfg.ylim = [-1 1];cfg.viewmode = 'vertical';cfg.artfctdef.zvalue.artifact = artifact;ft_databrowser(cfg, data_epoch);

图4.Databrowser显示了(12个试次之一的)滤波后的数据。

现在可以删除确定为伪影的试次。

cfg = [];cfg.artfctdef.zvalue.artifact = artifact;cfg.artfctdef.reject = 'complete';data_epoch = ft_rejectartifact(cfg, data_epoch);

锁时分析

使用ft_timelockanalysis计算试次的平均值。随后，可以使用FieldTrip绘图函数来查看数据，也可以简单地使用MATLAB绘图函数进行可视化。

cfg = [];data_timelock = ft_timelockanalysis(cfg, data_epoch);

输出是带有以下字段的data_timelock数据结构

data_timelock =      avg: [2x225 double]      var: [2x225 double]     time: [1x225 double]      dof: [2x225 double]    label: {2x1 cell}   dimord: 'chan_time'     opto: [1x1 struct]      cfg: [1x1 struct]

接下来，绘制A-10s到A+35s的平均O2Hb和HHb轨迹。根据fNIRS惯例，O2Hb为红色，HHb为蓝色。

time = data_timelock.time;O2Hb = data_timelock.avg(1,:);HHb = data_timelock.avg(2,:);figure;plot(time,O2Hb,'r'); hold on;plot(time,HHb,'b');legend('O2Hb','HHb'); ylabel('\DeltaHb (\muM)'); xlabel('time (s)');