0.   引言

功能性近红外光谱（functional near-infrared spectroscopy，fNIRS）技术一般指利用人体血液中氧合血红蛋白（oxygenated hemoglobin，HbO₂）对特定范围波长的近红外光相比脑部其他组织有特异吸收的特性来测量并记录近红外光强的变化的脑成像技术，因其具有高时间分辨率、安全无损、便携性高、应用场景广泛等优点^[1]，近年来在脑神经科学、脑机接口、精神及神经类疾病等领域得到广泛应用^[2]，深受研究者的关注。然而，在利用fNIRS记录被试者的脑功能活动时，会产生不可避免的噪声干扰，若不对干扰信号进行合理校正，很容易导致其与大脑活动引起的真实血流动力学响应相混淆，对后期数据分析的结果产生影响^[3-4]。因此，在fNIRS建模分析过程中，滤除冗余信息干扰、提高信号质量成为了一项必要环节。实际场景采集到的fNIRS信号噪声成分较为复杂，含有生理噪声、运动伪迹、系统噪声等多重干扰信号^[5]，生理噪声一般具备一定的规律性，可通过一定频段的带通滤波器滤除，系统噪声可用短距离参考通道记录浅层噪声，然后用回归的方法从神经记录中校正^[6]。而运动伪迹的校正则是较为棘手的问题，实验范式的多样性可能导致运动伪迹含有与脑活动信号类似的特征，需要根据实际情况对伪迹进行识别，不同类型的运动伪迹校正的方法也不同^[7]，因此运动伪迹的校正算法也是近几年研究的热点^[8]。

早期的研究者会采用舍弃运动伪迹较多的数据段和通道，但这样会使数据信息丢失^[9]。还有基于硬件的方法，使用加速度计或短距离通道记录运动伪迹，然后利用回归等算法对伪迹进行校正，但是该方法增加了设备的成本，且基于硬件的方法很难完全抑制运动伪迹的产生^[10]。因此，近年来，多种基于后处理技术的运动伪迹校正算法被陆续提出，较为常用的有运动伪迹减少算法（movement artifact reduction algorithm，MARA）、时间导数分步修复（temporal derivative distribution repair，TDDR）、小波滤波（wavelet filtering，WF）、基于相关性的信号质量提升算法（correlation-based signal improvement，CBSI）等。MARA算法^[11]可以很好地校正较大幅度的运动伪迹，但对信号中的小幅度运动伪迹不能很好地滤除。TDDR算法^[12]依靠对时间导数方差的值来计算鲁棒权值，可以很好地校正尖峰和基线漂移，但高频生理噪声和随机噪声的存在都会使方差计算产生误差，从而降低去伪迹效果。WF算法^[13-14]校正尖峰的效果良好，但其去噪性能会随阈值参数的选择较为复杂，且有时会加剧信号中的基线漂移；CBSI算法^[15]的去伪迹效果受数据实验范式的影响较大。考虑到以上经典算法的局限性，近年来有研究者提出将两种或多种算法相结合的方式来弥补单一算法的不足，如Spine-SG算法将样条插值算法与最小二乘法相结合，先用MARA算法校正基线突变，再用SG滤波校正尖峰，实验结果表明，Spine-SG算法的去伪迹性能相比单一算法有显著提升^[16]。

近年来，随着计算机科学与人工智能等领域的发展，深度学习技术已成为一种强大的信号处理工具^[17]。相比于传统算法，深度学习方法不依赖于参数主观设定，可以自动提取信号中的伪迹特征并抑制噪声，有着更强的伪迹识别能力和去噪性能^[18]。而由滤波权重和非线性单元组成的卷积自编码框架可以使用池化层和上采样算子对信号进行多尺度分解来提取多分辨率支持下的特征，非常适合处理时间序列信号，可以更加精准地识别fNIRS信号的运动伪迹^[19]。陈建等将卷积神经网络和自编码器相结合，通过卷积的方法构建深度神经网络来学习从含噪心电信号到干净心电信号的端对端映射，结果表明，该方法可在消除心电信号中复杂噪声的同时完整保留心电信号的形态^[20]。Qiu等结合初始模块提出了一种用于不同噪声环境下海事图像去噪的轻量级自动编码器，并利用经过调整的预训练U-Net模型完成了无人水面车辆拍摄的海洋图像的语义分割任务^[21]。滑世辉等在卷积自编码器的基础之上，添加更多的层数并加入跳跃连接构造，从提升网络的深度，能够提取并结合数据中深层与浅层特征，增强对地震数据随机噪声的处理能力，更好地恢复有效信号的细节，经过实验证明该方法对于不同水平的随机噪声的压制均有优异的表现^[22]。本文提出了一种基于卷积自编码机的fNIRS信号运动伪迹校正算法，首先通过生成大量仿真模拟数据集训练网络模型，提升模型识别信号伪迹特征性能，然后与几种经典算法在模拟数据和实验数据进行对比分析，验证了所提出算法——多层卷积自编码机（multilayer convolutional autoencoder，MCAN）算法在校正fNIRS信号中运动伪迹时的有效性。

1.   算法原理

卷积自编码机（convolutional autoencoder，CA）是一种无监督学习的CNN结构，由Vincent等提出^[23]，该模型在图像和信号去噪领域有着广泛的应用。CA通过对网络输入有噪声污染的样本，对其进行编码与解码后，最终重构出不含有噪声的信号。在CA网络中，对于输入数据引入随机噪声，即：

式中：E(ϑ)表示噪声类型。CA网络中，输入是含有噪声的特征向量，损失的输出仍然为预测输出，x_n为不含噪声的输入，即期望输出。CA网络结构层数较少会节省训练时长，但会影响最终的去噪性能^[24]，因此，本文增加网络层数并重新搭建网络架构，提出MCAN算法，设定新的仿真模拟数据的参数来训练网络，将仿真模拟数据集输入编码机进行重构，利用池化层和上采样算子的交替组合，对输入信号进行二尺度分解，提取多分辨率特征，更加精准地识别并校正fNIRS信号的运动伪迹。本文所提出的MCAN算法框图如图 1所示。

图  1  MCAN算法fNIRS信号运动校正算法框图

Figure  1.  Block diagram of the MCAN algorithm fNIRS signal motion correction algorithm

下载: 全尺寸图片 幻灯片

1.1   编码部分

MCAN算法的编码部分由6个卷积层和6个最大池化层及激活函数交替组合而成，该部分可将fNIRS信号进行特征降维，滤除冗余信息。卷积层相当于“特征提取机”，对fNIRS信号特征进行编码，并校正运动伪迹，卷积层的输出可用下式表示：

式中：x_j^l表示第l个卷积层第j个卷积核对应的特征向量；M_j表示当前神经元的接受域；w_ij^l表示x_j^l的第i个加权系数；b_j^l表示对应的该卷积核的偏置系数；f代表非线性函数，其计算公式为：

该部分最大池化层在完成特征降维的同时保留有用的信息，下采样层还具有二次特征提取的作用，其计算公式为：

式中：down(⋅)表示下采样；β_j^l、b_j^l分别表示下采样的加权系数和偏置系数。

1.2   解码部分

解码部分由6个卷积层和6个上采样层交替组合而成，可对fNIRS信号特征的深度表示进行重构，利用卷积操作和上采样操作，输出一个和输入信号同样尺寸的fNIRS信号。上采样层输出可用下式表示：

式中：up(⋅)表示上采样；⊗表示上采样算子；β′_j^l、b′_j^l分别表示上采样的加权系数和偏置系数。第l层的转置卷积层输出可表示为：

式中：x′_j^l表示第l层反卷积层第j个卷积核对应的特征向量；M_j表示当前神经元的接受域；⊗表示反卷积操作；w′_ij^l表示第l个反卷积层第j个卷积核的第i个加权系数，b′_j^l表示对应的偏置系数；而f为非线性函数。

1.3   网络训练

MCAN算法的训练流程示意如图 2所示。MCAN中使用均方误差损失函数计算Loss值，计算公式为：

图  2  MCAN算法训练流程示意图

Figure  2.  Schematic diagram of the MCAN algorithm training process

下载: 全尺寸图片 幻灯片

式中：y_i代表真实值；$ {\hat y_i} $代表预测值。最后，损失函数最小时，输出数据中超出标准差阈值和幅值阈值的残余运动伪迹数量。

本文以8:1:1的比例将仿真模拟数据集随机分成训练、验证和测试数据集。MCAN在训练数据集上训练，并在验证数据集上完成验证，训练完成后，使用测试数据集评估其性能，在生成的模拟数据集上完成训练和测试后，再使用真实实验数据对MCAN算法的性能进行进一步验证。本文的采用Adam优化算法代替传统的随机梯度下降算法来最小化损失函数，Adam优化算法为不同的参数设计独立的自适应性学习率，通过梯度的一阶矩估计和二阶矩估计的计算来实现，且Adam算法在保持较低内存需求的同时具有很高的计算效率。使用反向传播方法来训练网络，学习率设置为0.0001，训练100个epoch，当Loss最小时保存相应参数，MCAN训练的可视化过程如图 3所示。

图  3  MCAN训练可视化过程

Figure  3.  Visualization of MCAN training

下载: 全尺寸图片 幻灯片

2.   算法验证

2.1   仿真数据验证

2.1.1   数据模拟

模拟生成的仿真实验数据由真实的血流动力学响应、运动伪迹信号、静息态fNIRS信号分量3部分叠加而成：

式中：Φ_s^(λ)(t)代表当波长为λ时的近红外光强度信号；Φ_HRF^(λ)(t)代表被试者真实的血流动力学响应；Φ_rs^(λ)(t)、Φ_MA^(λ)(t)分别表示静息态fNIRS信号和运动伪迹信号。Φ_HRF^(λ)(t)由Gamma函数模拟获得^[25]，计算公式为：

式中：

参数α用于调整HRF(t)的幅值；β设置反向脉冲的幅度；τ₁和τ₂用于调整HRF(t)的形状，其中t为时间；ρ₁和ρ₂用于调整HRF(t)的阈值。本文中α＝1282，β＝0.17，τ₁＝1，τ₂＝1，ρ₁＝－0.5，ρ₂＝0.35，HRF(t)的阈值从34~37 μm⋅mm之间的均匀分布中随机生成，这样可以使MCAN更好地从样本中学习；Φ_rs^(λ)(t)为静息态fNIRS信号，由一个五阶的自回归模型进行模拟；Φ_MA^(λ)(t)为运动伪迹信号，包括尖峰、基线漂移、短低频振荡等，由Laplace函数模拟生成^[26]，仿真数据模拟过程示意如图 4所示。

图  4  仿真数据模拟示意图

Figure  4.  Schematic diagram of the simulation data

下载: 全尺寸图片 幻灯片

2.1.2   数据处理

用所提出的MCAN算法和MARA、CBSI、WF、TDDR、Spline-SG算法分别对模拟数据进行伪迹校正，对比去噪效果。CBSI和MCAN处理的是血红蛋白浓度值，需将光强信号转化为光密度（optical density，OD）信号，之后再检测剩余运动伪迹数量。TDDR和Spline-SG分别用Homer2工具包中hmr Motion Correct TDDR和hmr Motion Correct Spline SG函数实现，设置参数分别为：p＝0.99，iqr＝1.5，CBSI算法使用NIRS_KIT工具包实现。MARA算法的去伪迹性能受p值影响，而WF算法的性能也受iqr参数影响^[9]。因此，为了公平比较各个算法的性能，应确定MARA和WF的最佳参数p和iqr，对这两种算法进行灵敏度分析^[6]，图 5绘制了两种算法的概率阈值参数和p值关于剩余伪迹数量n和关于均方误差（mean square error，MSE）的变化曲线图，最终确定MARA算法中p＝0.99，WF算法中iqr＝1.5。

图  5  MARA和WF算法灵敏度分析结果

Figure  5.  Sensitivity analysis results of MARA and WF algorithms

下载: 全尺寸图片 幻灯片

2.2   真实数据验证

本文采用NITRC网站（https://www.nitrc.org/frs/downloadlink.php/11882）中的真实公开数据集的HbO₂信号作为实验数据，该数据使用TechEn CW6设备采集5名被试者任务态的fNIRS信号，实验任务要求被试者进行多种不同类型的动作（上下点头、侧向点头、大声朗读、右扭、左扭、快速左右摇头和扬起眉毛），每个运动试验进行3 s，每种运动类型试验重复5次，试验间随机间隔为5~10 s。设备采样频率为50 Hz，包含2个光源和4个探测器，长距离探测器和短距离探测器到光源的距离分别为30 mm和8 mm，探头位于每个被试者头部的运动区域。

分别用本文提出的MCAN算法和MARA、CBSI、WF、TDDR、Spline-SG五种算法对含有运动伪迹的信号进行校正。其中MARA、WF、TDDR、Spline-SG处理的是OD信号，CBSI和MCAN处理的是血氧浓度信号，使用Homer2中的hmrConc2OD函数将其转换为OD值，再用hmr Motion Artifact函数检测剩余运动伪迹。6种算法的处理流程如图 6所示。

图  6  六种算法的数据处理流程

Figure  6.  Data processing flow diagram for six algorithms

下载: 全尺寸图片 幻灯片

2.3   评价指标

采用块平均HbO₂浓度信号的均方误差（mean square error，EMS）、信噪比（signal-to-noise ratio，RSN）、皮尔逊相关系数的平方R²和峰峰误差E_p来比较上述算法的校正运动伪迹的性能。4种指标的计算公式为：

式中：$ {s_x} = \sqrt {\frac{1}{{N - 1}}{{\sum\limits_{i = 1}^N {\left[ {x\left( {{t_i}} \right) - \left\langle {x\left( t \right)} \right\rangle } \right]} }^2}} $；

x(t)和$ \tilde x\left( t \right) $分别为不含运动伪迹的信号和已经过校正的信号；x_p、$ {\tilde x_p} $分别表示信号x(t)和$ \tilde x(t) $的最大；N为信号x(t)和$ \tilde x(t) $的总采样点数。EMS和RSN用来评价信号x(t)和$ \tilde x(t) $的一致性，R²评价两者的相似性。EMS和E_p越小，算法去运动伪迹的性能越好；RSN越大，去运动伪迹的性能越好；R²越接近1，去运动伪迹的性能越好。

3.   实验结果与讨论

3.1   仿真数据验证结果

将上述6种算法预处理后的仿真块平均HbO₂信号分别计算其EMS，RSN，R²和E_p四个评价指标来综合评价各个算法的去伪迹性能。使用Homer2工具箱内的hmrMotionArtifact函数检测每种算法剩余运动伪迹，结果如图 7(d)所示。对每组数据的EMS，RSN，R²和E_p进行平均，结果如图 8所示。

图  7  MCAN算法处理后的仿真数据结果

Figure  7.  The simulation data results processed by the MCAN algorithm

下载: 全尺寸图片 幻灯片

图  8  仿真数据参数对比

Figure  8.  Comparison of simulation data parameters

下载: 全尺寸图片 幻灯片

采用MCAN算法对3类不同血流动力学响应波形信号的去噪结果如图 7(a)~(c)所示，从图 7(a)~(c)中可以看出MCAN算法可以校正大部分尖峰、基线漂移和短低频振荡，并且可以很好地还原有用信号。从图 8可以看出，MCAN算法和MARA、CBSI、WF、TDDR、Spline-SG五种算法对运动伪迹均有不同校正效果，与校正前相比，MCAN、Spline-SG、TDDR、WF的EMS，RSN和R²都有很大改善，其中MCAN算法表现最优；TDDR和Spline-SG在E_p指标上相比于校正前也有很大提升；MARA、TDDR、MCAN在R²指标中表现较为优异，WF表现最差；CBSI在各项综合指标中表现较为普通；综合EMS，RSN，R²和E_p四项指标来看，MCAN算法表现最为优异，性能最好。

3.2   真实数据验证结果

图 9(a)为用MCAN算法处理含有运动伪迹的真实信号的示例，由图 9可看出MCAN很好地校正了信号中尖峰、短低频振荡、基线漂移运动伪迹，校正后的信号接近真实任务态信号。用MCAN和MARA、CBSI、WF、TDDR、Spline-SG五种算法对真实数据校正运动伪迹后，计算块其平均HbO₂信号的EMS，RSN，R²和E_p四项评价指标，并将每个被试者数据的所有通道平均，如图 10所示。由图 10可以看出，与校正运动伪迹之前相比，TDDR、Spline-SG、MCAN的EMS、RSN和E_p都有明显的改善，其中MCAN算法最佳；其中Spline-SG在R²上的表现最佳，MCAN和MARA也有较好的表现，WF表现最差；在该数据集上CBSI的表现较为中等。相比于其他算法，MCAN算法在EMS，RSN，R²和E_p四项综合评价指标表现最优，只有在R²指标上比Spline-SG算法表现稍差。

图  9  实验数据数据处理结果

Figure  9.  Diagram of experimental data processing results

下载: 全尺寸图片 幻灯片

图  10  真实数据参数对比

Figure  10.  Comparison of real data parameters

下载: 全尺寸图片 幻灯片

4.   结论

本文提出了一种基于卷积神经网络的fNIRS信号运动伪迹校正算法——MCAN算法。为了使MCAN网络更好地学习伪迹特征，生成大量模拟fNIRS数据来训练网络。当MCAN应用于模拟和真实fNIRS数据时，通过将训练的MCAN网络与其他传统伪迹校正算法进行比较来验证其校正伪迹的性能。实验结果表明，与其他算法相比，MCAN算法能够校正信号中大部分尖峰、基线漂移和短低频振荡运动伪迹，且处理数据后产生的剩余运动伪迹数量最低，在EMS、RSN、E_p和R²四项综合指标上表现突出。

在各个算法处理数据过程中，TDDR的去噪性能比MCAN算法较差，因为TDDR对于数据中时间导数值较大的噪声（如尖峰，基线漂移）有很好的校正效果，而对于短低频振荡以及较小尖峰的校正效果则较为普通，而经过模拟数据训练的MCAN则更容易将这类伪迹校正。MARA和Spline-SG都需要对运动伪迹进行检测，然后用样条曲线拟合运动伪迹，并从信号中减去该曲线达到去运动伪迹的目的，但这类算法容易将部分有效信号去除，导致校正过度^[27]，MCAN算法则更能起到学习和识别伪迹和真实血流动力学响应的作用，不易产生校正过度现象，更容易保留有用的信号细节。CBSI算法在EMS，RSN，R²和E_p四项指标上表现较为普通，这是因为该算法的去噪性能容易受到数据的实验任务的影响。WF算法在指标R²上表现较差，该算法处理基线漂移和短低频振荡时可能导致信号平滑，对后续数据处理和分析产生影响。

经过仿真数据和实验数据的验证，与MARA、CBSI、WF、TDDR、Spline-SG五种算法对比发现，MCAN算法去运动伪迹性能较好，可作为一种新算法应用于fNIRS信号的预处理阶段。