基于fNIRS信号的运动伪影去除算法综述

A Review of Motion Artifact Removal Algorithms Based on fNIRS Signals

摘要: 近红外脑功能成像技术(functional near-infrared spectroscopy, fNIRS)是一种新型、无创的脑功能检测技术。由于光源探测器与头皮之间的相对运动会引起信号的运动伪影，严重影响信号质量，这也成为了如今该技术的一个重要研究方向。经过多年的研究与发展，人们已经提出小波滤波、样条插值、主成分分析、基于相关的信号改进、深度学习等运动伪影去除技术。本文旨在系统地总结已有的运动伪影去除技术，并阐述了当前普遍使用的定量评价指标：均方误差(MSE)、皮尔逊相关性(R)、平均绝对百分比误差(MAPE)、曲线下面积(AUC)。最后总结了在运动伪影去除领域研究的不足，并提出展望。

Abstract: Functional near-infrared spectroscopy (fNIRS) is a new, non-invasive technique for detecting brain function. Because the relative motion between the light source detector and the scalp causes the motion artifact of the signal, the signal quality is seriously affected, which has become an important research direction of this technology. After years of research and development, people have proposed wavelet filtering, spline interpolation, principal component analysis, Deep learning, correlation based signal improvement and other motion artifact removal techniques. The purpose of this paper is to systematically summarize the existing motion artifact removal techniques, and describe the commonly used quantitative evaluation indexes: mean square error (MSE), Pearson correlation (R), mean absolute percentage error (MAPE), and area under the curve (AUC). Finally, the shortcomings of the research in the field of motion artifact removal are summarized, and the prospect is put forward.

1. 引言

近20年以来，功能性近红外脑成像(functional near-infrared spectroscopy, fNIRS)近作为一种无创、便携的神经成像技术[1]，日益受到大众的关注。相较于脑电图(EEG)和核磁共振成像(fMRI)，具有更高的空间分辨率与更低的成本，并且对伪影的敏感程度较低[2]-[4]。它使用了两段不同的波长(一般使用735 mm与850 mm)的近红外光，从紧贴在头皮的光源发出，并使用修改的比尔–朗伯定律将其转换为氧合血红蛋白(HbO)和脱氧血红蛋白(HbR)的浓度变化。浓度变化体现为表层皮质层中的神经元激活变化，来监测大脑神经活动相关的血流动力学变化[5]-[7]。

运动伪影是指由患者移动引起的生物电位测量误差，这在非接触式测量中是不可避免的。然而，许多因素会导致伪影的产生。一些研究人员进一步发现，面部肌肉运动，包括扬起眉毛，可以导致运动伪影[8]。尤其是在头部有极端运动的人群中，例如婴儿的头部运动、前额抖动以及任何皮肤的移动都有可能引起的运动伪影，导致接收信号的幅度发生突变，从而造成不同类型的伪影，这会显著地影响fNIRS记录的光学信号的质量[9]。人们普遍认为，运动伪影源于两个主要来源：电容变化和摩擦电压[10] [11]。尖峰伪影通常是在很短的时间内幅值出现了较大变化，并且很快又恢复到原始幅值，可以视为脉冲型信号。如图1所示，一般来说尖峰型伪迹的幅度最大可以比正常信号大两个数量级，对于这类伪迹可以采用平滑滤波的方式来校正。基线漂移类似像是一个阶跃信号，这类伪迹又可以细分为基线突变与缓慢漂移两种类型。出现这类伪迹时，信号基线将从一个幅值变化到另一个幅值，且不再恢复到原幅值[12]。

随着时代的发展，产生了各种不同的抑制伪影的方法。在众多研究者不断地开发下，滤波器和去伪技术已经运用在实际生活中。早期的模拟滤波器随着电子元件的发展和电路理论的演进而逐渐形成。早期，当测量的信号被运动伪影显著破坏时，研究人员跳过分析或丢弃数据集。在上世纪60年代，Kálmán首次引入了卡尔曼滤波的概念，该方法运用最优化自回归数据处理算法[13]。随后，Stearns提出的自适应滤波开始兴起，出现了一种动态调整滤波器参数以适应环境变化的滤波技术[14]。Morlet和Grossmann提出了小波变换的概念，解决了傅里叶变换的困难问题，为滤波技术的发展产生了重大的突破[15]。在21世纪初，研究者又陆续提出一些不需要额外输入的矫正方法，例如样条插值法等。随后，深度学习在滤波中的应用开始广泛，通过卷积神经网络(CNN)来实现[16]。

如今，为了进一步提高伪影去除的效果，还原真实的血液动力学反应，更多人将运用多种混合算法来处理应对不同的伪影特点，不同的技术可能在处理特定类型的运动伪影时表现更好，因此更多考虑结合多种技术可以弥补各自的局限性，提高信号处理的效果。

Figure 1. The characteristics of different types of motion artifacts

图1. 不同类型运动伪影的特征

2. 近红外光谱数据

首先，需要准备fNIRS设备，包括近红外光源和探测器。将光源和探测器的阵列放置在被试者头部的特定区域，通常以一定间隔形成阵列。如图2所示，这个阵列的位置通常按照国际10-20系统的位置进行定位[17] [18]，以确保数据的标准化和可比性。这些设备通常以探测头的形式放置在被试者的头皮表面上。头盔或头带等装置通常用于确保探测器的正确定位和固定。

Figure 2. International 10-20 system location positioning

图2. 国际10-20系统位置定位

其次是实验设计，根据研究目的，设计相应的实验任务或刺激。例如，可以通过观察被试者执行特定的认知任务、感觉刺激或情绪激发，来引起脑区域的活动变化。

采集到的fNIRS数据需要进行后续的分析。这可能包括对数据进行预处理，如去噪、去趋势、运动伪影校正等，然后使用数学模型来推断血红蛋白和氧合血红蛋白的浓度变化。统计分析通常用于确定在实验条件下是否存在显著的脑活动变化。

3. 运动伪影矫正技术

3.1. 维纳滤波

维纳滤波是Wiener于1940年提出的技术，已广泛应用于许多不同的应用领域，包括生物医学信号，语音等，可以用于去除图像中的噪声，提高图像的清晰度，并恢复图像中丢失的细节。维纳滤波与自适应滤波一样，也是一种均方意义上的最优滤波技术，但它利用信号的统计特性来估计滤波器系数而不需要附加的传感器信息[19] [20]。

维纳滤波器基于最小均方(MMSE)准则，旨在最小化原始信号与估计信号之间的均方误差。其关键思想是通过估计信号的条件期望值来滤除噪声，从而重建或增强原始信号[21]。首先需要将失真fNIRS信号$y\left(n\right)$ 进行建模：

$y\left(n\right)=w\left(n\right)+x\left(n\right)$ (1)

其中，$x\left(n\right)$ 是真实测量的脑信号，$w\left(n\right)$ 是运动伪影。然后做出假设：$x\left(n\right)$ 与$w\left(n\right)$ 为静止信号，并且不相关[19]。维纳滤波器使得真实近红外信号$x\left(n\right)$ 和估计信号$\widehat{x}\left(n\right)$ 之间的均方误差$e\left(n\right)$ 最小化即：

$E\left[x\left(n\right),w\left(n\right)\right]=0$ (2)

由此，利用正交性原理可以得到使均方误差最小的最佳滤波器[22]。

$E\left[e\left(n\right),y\left(n\right)\right]=E\left[\left(x\left(n\right)-\widehat{x}\left(n\right)\right),y\left(n\right)\right]$ (3)

式(3)把原式转化为傅里叶域上的代数方程[19]。最后，在傅里叶域上的维纳滤波器可以表示为：

$G\left(w\right)=\frac{P_x\left(w\right)}{P_x\left(w\right)+P_w\left(w\right)}$ (4)

$P_x\left(w\right)$ 与$P_w\left(w\right)$ 为功率谱密度。

随着对噪声模型和信号特性理解的提高，未来的研究可能会更加精确地建模噪声和信号的统计特性，以提高维纳滤波器在复杂场景中的适用性和性能。维纳滤波器在特定条件下是一种有效的信号处理工具，但在面对噪声模型未知或信号特性复杂的情况下，可能需要考虑其他滤波方法或者结合多种方法来取得更好的效果。

3.2. 样条插值

样条插值的方法首次由Scholkmann在2010年提出，在去除基线漂移伪影方面有着重要的作用。该技术主要步骤分为以下三步：(1) 运动伪影的识别；(2) 运动伪影的去除；(3) 时间序列的重构。它只作用于检测到的运动伪影，而不改变信号的其余部分[23]。首先计算离散时间序列的移动标准差s(t)，公式如下：

$S\left(t\right)=\frac{1}{2k+1}{\left[{\displaystyle \sum _{j=-k}^k{x}^2\left(t+j\right)}-\frac{1}{2k+1}{\left({\displaystyle \sum _{j=-k}^{k}x\left(t+j\right)}\right)}^2\right]}^{\frac{1}{2}}$ (5)

其中，W为滑动窗口大小，W = 2k + 1。$t=k+1,k+2,\cdots ,N-k$ ，其中N是时间序列x(t)的长度。然后开始检测信号片段中运动伪影的起点与终点。需要人为设定一个阈值T，将MSD计算出的结果与阈值T作比较。S(t)中所有大于t的部分将被识别为运动伪影部分。但这种检测方法需要依赖于人为设置参数T和W，所以会导致算法相对复杂[10]。

在识别完伪影之后，需要对伪影部分分别进行样条插值处理。样条插值依赖于一个参数p，它决定了样条函数滤波的程度效果。如果p = 0，则为直线插值，如果p = 1，则为三次样条插值[24]。建议对应于参数p = 0.99，在去除运动伪影方面可靠有效[25]。运用插值法对伪影部分进行拟合后，做样条插值函数的减法。用原始信号减去拟合之后的信号值，得到一段新的过滤信号。最后，为了保持过滤后的信号与其余信号的连续性，需要对整个信号序列进行重构。根据前一段的均值和目标段的均值，对每一段伪影部分进行平移，变成一段连续的时间序列信号。完整流程如下图3所示。

Figure 3. Flowchart of spline interpolation algorithm

图3. 样条插值算法流程图

3.3. 主成分分析

主成分分析(PCA)是一种常见的数据分析方式，常用于高维数据的降维，可用于提取数据的主要特征分量。PCA通常用于降低大型数据集的维数，方法是数据集中的指标数量变少，并且保留原数据集中指标的大部分信息[26]。该方法对由N个测量值组成的原始数据集进行正交变换，产生N个不相关的分量。由于运动伪影的振幅通常比正常的生理fNIRS信号大得多，它们应该构成数据方差的很大一部分[24]。假设前M个分量表示运动伪影引起的方差。因此，从信号中去除前M分量可以校正运动伪影。如下图4、图5所示，随机生成一组时序信号，信号数据的投影到主成分空间的过程。方差按PC1 > PC2 > PC3的顺序递减，假设PC1-PC2平面已包含大部分运动伪影信息，总方差保持不变。

Figure 4. The projection process of signal data

图4. 信号数据的投影过程

Figure 5. Data distribution of the three principal component planes

图5. 三个主成分平面的数据分布

Zhang等人将PCA应用于fNIRS信号的基线，以识别系统干扰的主要空间特征向量。他们假设基线信号主要包含全局干扰的空间模式，并且相应的空间分量具有更高的能量，更具有全局性。然后通过将fNIRS数据投影到所识别的空间特征向量的正交子空间上来进行空间滤波[27]。主成分分析算法可以分为以下几个步骤：

通过将来自所选择的通道的光强度时间过程放入N_s = N_t矩阵中来形成原始数据矩阵，其中N_s是通道的数量，N_t是在不同时刻的数据点的数量。将估计的基线和刺激浓度的矩阵表示为H_base和H_stim，空间相关矩阵为C_base，可由下式表示：

$C_{base}=\left(1/N_t\right)H_{base}{H}_{base}^{\text{T}}$ (6)

C_base的特征分解式：

$C_{base}=U_{base}{\Sigma }_{base}{U}_{base}^{\text{T}}$ (7)

其中U_base被称为基线空间特征向量的列u_i的正交矩阵。这里的矩阵${\Sigma }_{base}$ 是一个包含特征值的对角矩阵。

全局空间干扰子空间是基线数据中的主导，假设对于较小的r值，是由C_base的前r个空间特征向量所决定。然后，通过将刺激数据投影到$U_{base,r}$ 的正交子空间表示为$U_{base,r}^{\perp }$ 上来实现空间滤波，以获得滤波后的矩阵${\tilde{H}}_{stim}$ 。

$\begin{array}{c}\tilde{H}=\left(1-U_{base,r}{U}_{base,r}^{\text{T}}\right)H_{stim}\\ =U_{base}^{\perp }{U}_{base,r}^{\perp }{}^{\text{T}}{H}_{stim}\end{array}$ (8)

最后，滤波的结果进行相干块平均，以提高信噪比相对于一个单一的块响应。

3.4. 小波滤波

小波滤波是Molavi提出的运动伪影去除技术[28]，小波滤波技术作为一种强大的信号处理方法，在多个领域都发挥着重要作用。随着技术的不断进步，小波滤波算法的改进和应用领域的拓展将进一步推动这一领域的发展。

小波滤波技术是试图去除所有噪声，保留所有信号，并不考虑它们的频率范围。其次，是在小波域对小波系数进行处理。其基本步骤包括以下三步[29]：

1) 小波分解：如图6所示，将其信号分解成不同尺度和频率的小波系数，通常通过多层级的分解来获取更详细的信息。

2) 阈值处理：对小波系数进行阈值处理，保留具有较大幅值的系数，而将较小的系数置零。这一步可以实现信号的稀疏表示，从而实现降噪的效果。

3) 小波逆变换：利用经过阈值处理后的小波系数重构原始的信号，得到滤波后的结果。

Figure 6. Schematic diagram of the wavelet decomposition process

图6. 小波分解过程示意图

小波变换基本上分为连续小波变换和离散小波变换。连续小波变换(CWT)将连续信号分解为小波。离散小波变换(DWT)与连续小波变换的不同之处在于小波是离散采样的。与连续小波变换相比，离散小波变换速度快，易于实现。离散小波变换技术包括通过一系列具有不同截止频率的不同尺度的滤波器来传递信号[30]。

阈值的选取也十分重要，它决定了该方法对信号滤波的效果[31]，以去除与噪声相关的低系数[32]。目前，已经提出了众多阈值确定的方法。BayesShrink阈值法、通用阈值法、极小化风险阈值、多假设检验法等[33] [34]。自适应阈值根据小波系数的局部特征来进行阈值化处理，最简单的是可以适应于不同尺度上系数特征的阈值[29]。

小波逆变换是滤波过程中的重要步骤，它用于将经过小波变换得到的小波系数重构成原始信号。小波变换将信号分解成不同频率的小波系数，而逆变换则是将这些系数重新组合以恢复原始信号。逆变换从最底层(最细节的尺度)开始，逐步进行逆向的重构。这个过程是一个逐层上采样和滤波的过程，逆序地应用滤波器，将小波系数重新组合[35] [36]。

小波滤波技术在多尺度分析和时频局部化方面的优势，在信号处理的各个领域都有着重要的应用和研究价值。但是该方法运用在工程实践中，仍然进度落后于理论，许多实验还是停留在仿真阶段[29]。

3.5. 鲁棒局部加权回归

鲁棒局部加权回归(rloess)首次由Cleveland提出，是一种用于回归分析的统计技术[37]。主要通过对采集的信号值进行多项式加权拟合处理，用最小二乘法进行估计。利用强局部加权回归算法对观测点进行拟合时，多项式的阶数、权值函数、迭代次数以及窗口宽度是很重要的。这种方法的核心思想是，不同数据点对于拟合的贡献应该与其距离目标点的距离成反比。

对每一个观测点$x_i\left(i=1,2,\cdots ,n\right)$ ，以x_i为中心选取窗口宽度f。然后选取合适的权函数W，W需要满足以下性质：

1. $W\left(x\right)>0,\left|x\right|<1$

2. $W\left(x\right)=W\left(-x\right)$

3. 对于$x\ge 0$ ，$W\left(x\right)$ 是非递增函数

4. $\left|x\right|\ge 1,W\left(x\right)=0$

通过权值函数可以得到权值ω_k(x_i)，使用带有权值ω_k(x_i)的加权最小二乘法对x_i进行d阶多项式拟合，得到拟合值y_i [30]。拟合式如下：

$y\left(x_i\right)=a_0\left(x_i\right)+a_1\left(x_i\right)x_i+\cdots +a_p\left(x_i\right)x_i^p,i=1,2,\cdots ,n$ (9)

计算初始拟合值的过程称为局部加权回归。大的残差导致小的权重，小残差会导致较大的权重。新权值和新适应度值的计算现在要重复几次。整个过程，包括初始计算和迭代，被称为鲁棒局部加权回归[38]。其中，残差可以定义为下式：

$e_i=y_i-{\widehat{y}}_i$ (10)

e_i是当前拟合值的残差，假设s是$\left|e_i\right|$ 的中位数。通过下式定义鲁棒性权重：

${\delta }_k=B\left(e_k/6s\right)$ (11)

通过采用鲁棒局部加权回归方法，既避免了复杂的参数调整过程，又不会出现“过拟合”或者“欠拟合”的问题[39]。整个拟合过程会不断地调整局部信号值的权重，使异常信号的影响降到最低。

该滤波技术的计算复杂度较高，因为对每个目标点都需要计算权重和进行拟合。此外，带宽参数的选择对拟合结果有较大影响，需要在应用中谨慎选择。总体而言，鲁棒局部加权回归是一种强大的回归方法，特别适用于处理具有局部信号。

3.6. 基于相关的信号改进

Cui等人在2010年提出的基于相关的信号改进方法(CBSI)，实时矫正信号效果良好，可以有效去除信号中的尖峰型伪影。核心思想是在信号处理过程中，利用信号之间的关联信息，以获得更准确、更稳定的信号结果。该算法需要额外的假设来进行信号校正：(1) 脑活动诱发的氧合血红蛋白(HbO)和脱氧血红蛋白(HbR)的浓度变化呈现负相关。(2) 运动伪影对HbO和HbR的影响相同。(3) 运动伪影与浓度的实际变化无关[40]。真实的HbO与HbR信号可表示为

$x=x_0+\alpha F+\text{Noise}$ (12)

$y=y_0+F+\text{Noise}$ (13)

x₀、y₀表示采集出的HbO与HbR信号，$\alpha$ 是正常数，F为噪声，Noise为fNIRS仪器造成的高频噪声。进一步假设运动相关噪声 F与仪器噪声相互独立，与脑信号无关[41]。为了满足x₀、y₀最大程度的负相关，满足以下关系：

$x_0=-\beta y_0$ (14)

x₀和F应具有最小相关性，因此将这两个量的相关性设为0，由此可以得到：

${\displaystyle \sum _t{x}^2}+\left(\beta -\alpha \right){\displaystyle \sum _{t}xy}-\alpha \beta {\displaystyle \sum _t{y}^2}=0$ (15)

随后，我们假设$\alpha =\beta$ ，可以把式子进行化简：

$\alpha =\frac{\text{std}\left(x\right)}{\text{std}\left(y\right)}$ (16)

之后便能简单得到结果x₀与y₀，其中$\alpha$ 是测量的氧合血红蛋白和脱氧血红蛋白信号的标准偏差的比率。

基于相关的信号改进的一个优点是它可以完全地自动化，并且有效地去除大的尖峰和基线漂移。但是在实际运用该方法的时候，也存在着局限性：由于过滤器的基本假设，CBSI只能应用于浓度变化[42]。在某些病理的状态下，在测量得到的氧合血红蛋白与脱氧血红蛋白并不会总是呈负相关，不满足该方法的基本假设。虽然这种方法改善了具有大量噪声的信号，但在噪声最小的情况下，它不会降低信号质量。

该方法在信号处理领域发挥着重要作用，为提高信号质量和提取有用信息提供了强大的工具。未来的研究方向包括更深入地挖掘信号之间的相关性，设计更复杂、适应性更强的算法，并将这些方法应用于更广泛的领域，以满足日益增长的实际需求。

3.7. 时间导数分布修复

时间导数分布修复(TDDR)是由Fishburn等在2019年提出的运动伪影去除技术。该方法可以同时对尖峰伪影和漂移进行处理，并且不需要调整参数，对数据有最小的假设[43]。虽然一些先前的研究认为大多数运动校正算法产生类似的结果，但其他研究表明时间导数分布修复(TDDR)是消除运动伪影的更好算法[44] [45]。

它基于三种基本假设：1) 非运动相关活动的导数近似正态分布。2) 大部分的波动不包括运动伪影。3) 当存在与运动伪影相关的波动时，比非运动波动的幅度大得多。

步骤一：就是计算测量信号的时间导数。其中，测量的脑信号(x)是时间(t)的函数，如式(16)：

$x_t={\displaystyle \sum _{i=1}^t{y}_t}$ (17)

步骤二：将原始信号对应的观测权值w_t初始化为1。

步骤三：主要通过不同方法更新迭代估计的观测权值。计算出波动的加权平均值$\mu$ 、估计值的绝对残差r_t，可以得到残差标准差的鲁棒估计：

$\sigma =1.4826\times \text{median}\left(r\right)$ (18)

然后使用残差标准差和调整常数计算每个观测点的标度偏差。最后使用Tukey的双权函数计算新的观测权重。

$w_t=\{\begin{array}{l}1-d_t^2\\ 0\end{array}$ (19)

步骤四：时间导数可由上个步骤的鲁棒权重得到。

$y_t=w_t\left(y_t-\mu \right)$ (20)

步骤五：对得到的时间导数进行积分，获得校正后的信号值。

但是，时间导数分布修复也存在一定的局限性。在原始信号中的高频成分无法有效地被矫正，高频成分的存在却会影响该方法的有效性。其次，它对典型的尖峰、漂移是非常有效的，但是对于缓慢漂移等伪影没有很好的效果，需要借助其他滤波器来去除。因为，未来该方法更需要关注是否能够与其他方法进行结合使用，发挥出更好的效果。

3.8. 独立成分分析

独立成分分析(ICA)是盲源信号分离的一种算法，它最重要的作用是通过线性变换，将一个混有多种信号的信息分解为不同的独立信号。独立成分分析可以从多个观测信号中识别出先前未知的、独立的原始信号，即从混合信号中恢复原始信号，而这些信号是不可直接观测到的[46]，如图7、图8展示了源信号的分离过程。源信号代表了产生的运动伪影信号，它们是相互独立的；混合信号为运动伪影与真实的脑信号所形成的采集数据。

Figure 7. Source signal timing diagram

图7. 源信号时序图

Figure 8. Mixed-signal timing diagram

图8. 混合信号时序图

ICA的核心假设是源信号是相互独立的。这是与其他信号分析方法(如主成分分析PCA)的主要区别之一。通过这个假设，我们可以寻找一个变换矩阵(W)，使得(Y = WX)，其中Y矩阵是ICA后的信号矩阵[47]。

首先是要将数据进行中心化处理：为了简化计算，通常首先对观测数据进行中心化，使每个变量的均值为零。然后进行白化：对数据进行线性变换，使得各个分量的协方差矩阵为单位矩阵，从而简化了独立性的估计。通过最大似然估计的准则，找到变换矩阵W，使得独立成分分析后的数据尽可能地相互独立。

ICA适用于光强度、光密度和浓度变化。但是，FastICA等一些常用的ICA算法基于信号的非高斯性，当信号的非高斯性不够显著时，这些算法的性能可能会下降。并且，如果遇到噪声的分布与信号的分布相似，ICA可能会误将噪声当作独立成分进行分离[48]。

3.9. 深度学习

近年来，在算力提升和大数据积累的支持下，基于神经元结构模仿和数据驱动的深度学习技术以其强大的非线性表征的优势，给传统模式分析带来了巨大的影响。与传统的各种去除运动伪影算法不同，深度学习的目标在于训练神经网络模型，直接从大量实例数据中学习如何把采集所得包含噪声的原始信号转变为真正的有效信号并输出结果。深度学习将为fNIRS运动伪影去除带来新的思想和思路[49]。

深度学习模型主要为：卷积神经网络(CNN) [50]、循环神经网络(RNN) [51]与长短期记忆网络(LSTM)、生成对抗网络(GAN) [52]、去噪自编码器(DAE)等。

2023年Gao等人提出了一种基于去噪自动编码器(DAE)的深度学习方法[53]，用于去除功能近红外光谱数据中的运动伪影。该网络采用反向训练，使用最小化损失函数。模拟噪声高频辐射的时间序列数据为输入数据，输出数据为模拟去除伪影的时间序列数据。损失函数是专门针对这一问题设计的。损失函数最小化了预测的fNIRS数据与真实数据之间的差异。

$L_{mes}=\frac{1}{n}{\displaystyle \sum {\left({\widehat{y}}_i-y_i\right)}^2}$ (21)

运动伪影检测是双重的：(1) 标准偏差(STD)超过标准偏差阈值或(2) 幅度变化(AMP)超过幅度阈值。运用了(1)的标准偏差损失函数(L_std)和(2)的幅度损失函数(L_amp)。

$L_{std}=\frac{1}{N_{std}}{\displaystyle \sum _{i=1}^T{\displaystyle \sum _{j\in j|{\Delta }_{i}d{c}_j>m}{\Delta }_{i}d{c}_j}}$ (22)

$L_{amp}=\frac{1}{N_{amp}}{\displaystyle \sum _{i=1}^T{\displaystyle \sum _{j\in j|{\Delta }_{i}d{c}_j>m}{\Delta }_{i}d{c}_j}}$ (23)

其中，

$m={\sigma }_{{\Delta }_1}dc\cdot C_{std}$ (24)

${\Delta }_{i}d{c}_j=d{c}_j-dc\left(j+1\right)$ (25)

在上述方程式中，dc是根据比尔–朗伯定律导出的从预测的血红蛋白值${\stackrel{⌢}{y}}_i$ 得到的投影光强值；C_amp和C_std是波幅和标准差变化的运动检测阈值。如果信号在时间间隔内的变化超过C_amp，则该时间点被标记为运动伪影。最后，训练深度学习网络的损失函数是上述四个个体损失函数的线性组合。

$L_{oss}=L_{mse}+{\theta }_1\cdot L_{var}+{\theta }_2\cdot L_{std}+{\theta }_3\cdot L_{amp}$ (26)

其中，${\theta }_1$ 、${\theta }_2$ 、${\theta }_3$ 分别被设置为1、1和10，以平衡每个项的大小。

将模拟数据集随机划分为训练集、验证集和测试集，比例为8:1:1。模型在训练集上进行训练，在验证集上进行验证，最后在测试集上评估性能。

该方法的不足之处在于模拟数据集可能对其他模型不利，未考虑运动伪影与诱发信号锁定的情况，未模拟硬件设计变化对去噪的影响。

3.10. 混合算法

如今随着时间的推移，越来越多的学者也意识到任何单一流程或单一方法本身都存在着其固有缺陷，比如：基于滤波的方法可能难以消除与血流动力学响应的频率相吻合运动伪影；而基于回归的方法又太过依赖辅助信号的质量与正确性。为了解决此类弊端、寻求更加稳健、更加广泛地去除复杂多变的运动伪影，一种新的混合算法形式开始被广泛运用[54]。目前来看，此类算法已然成为去除运动伪影的有效工具，并逐步地成为主流方向。此思路的核心思想利用“协同”与“互补”的思想，使用两种及以上算法形成一个完整流程，解决信号中不同特征与来源的噪声。

第一阶段为大幅值伪影的检测和剔除：采用对尖峰噪声敏感的方法，例如中值滤波、数学形态学滤波或滑动窗口标准偏差算法对信号进行初步筛查，并进一步修改和删除其内幅值大，易识别的脉冲型运动伪影。

第二阶段为去除残余伪影与基线漂移：用对缓变漂移有较好效果的处理方法：如小波变换、主成分分析或独立成分分析等，去除已经经过预处理的信号的缓慢漂移等残留噪声。

Hayder等在2023年提出的结合小波和基于相关性的信号改进(WCBSI)的算法[55]，WCBSI算法引入小波方法将它与CBSI算法加以组合，对运动伪影校正并应用于多通道数据中。该算法可以应用于所有指标、状态下的HbO数据，并且是一个自动化的算法只需要一个参数。但是WCBSI的运行时间并不是最优，虽然与多数算法是同一量级的，但是仍然可能存在在大规模数据的情况下运行效率不高的问题。

总之，混合算法不是多种方法的堆砌，而是对各种方法及其相关原理的理解和把握基础上的巧妙的拼合。它进一步扩大了可处理的频带范围，增强了适应各种运动伪影类型的能力，改善了去伪影效果的效果稳定性等，是当前fNIRS运动伪影去除领域一个最为重要的具有吸引力的发展趋向。

综合以上所有方法，表1中列举了各个方法的优缺点。

Table 1. The advantages and disadvantages of current motion artifact removal methods

表1. 目前运动伪影去除方法的优缺点

目前运动伪影去除方法	优点	缺点
维纳滤波	1. 在已知噪声特性时效果很好	1. 对非线性和非平稳的伪影信号效果较差
样条插值	1. 方法直观、计算速度快，被广泛集成于各种工具包(Homer3等)。2. 对修复漂移型伪影效果显著。	1. 依赖于运动检测的准确性。如果伪影识别错误，会留下漂移伪影或对正常信号进行插值
主成分分析	1. 无需先验知识。2. 能有效分离方差最大的成分。	1. 假设性过强：假设伪影是信号中方差最大的来源。 2. 需要主观判断选择去除哪些成分。
小波滤波	1. 优秀的时频局部化能力，非常适合处理尖峰类伪影。2. 能同时处理不同类型的伪影(尖峰和基线漂移)	1. 需要选择小波基函数和分解层数，参数设置需要经验。2. 阈值选择不当会过度平滑或残留伪影。
鲁棒局部加权回归	1. 对缓慢的漂移型伪影去除效果很好。2. 由于其“鲁棒性”，使其对尖峰不敏感，拟合曲线不会受尖峰过度影响。	1. 计算成本相对较高、花费时间较长，尤其对于长时间序列数据。2. 需要选择时间窗口参数，窗口大小直接影响去噪效果。
基于相关的信号改进	1. 由于HbO和HbR的负相关性，能很好地保留生理信号。2. 无需检测运动段，整体计算简单。	1. 主要针对基线漂移型伪影，对尖峰型的处理能力有限。2. 需要假设HbO和HbR在运动伪影上高度相关。
时间导数分布修复	1. 无需检测运动段，全自动处理。 2. 通过导数的分布来修复信号去除伪影，能同时处理尖峰和漂移。	1. 方法相对较新，在更广泛数据集上的普适性有待进一步验证。
独立成分分析	1. 利用盲源分离方法，能将信号分解为各个独立的成分，能够区分伪影和生理信号。	1. 成分选择需要依赖主观判断。
深度学习	1. 能处理复杂的伪影模式。2. 无需手动输入参数或阈值。	1. 需要大量的高质量训练数据，计算资源需求高
混合算法	1. 能够灵活地处理尖峰型、漂移型伪影。	1. 取决于使用者具体选择的方法。

4. 评价指标

现如今，对于运动伪影的去除技术已有了很大的进展。在前人的各项研究中，定义了各项指标定量比较伪影去除方案的性能，比较真实的血液动力学反应(HRF)与校正后的HRF进行比较。通常，普遍使用的指标如下：

1) 均方误差(MSE)是表示预测值与真实值之间差异的指标。测量参考信号和被头部运动校正信号之间的非比例平均偏差。MSE的值越小，表示模型的预测值与真实值之间的差异越小，伪影矫正的性能越好[56]。但MSE的也有一定的局限性，是它对异常值比较敏感，因为它对差异值进行了平方操作，使得较大的差异对整体误差的贡献更为显著。

$\text{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\left(y\left(i\right)-\widehat{y}\left(i\right)\right)}^2}}{N}}$ (27)

2) 平均绝对百分比误差(MAPE)：是一种用于衡量预测值与真实值之间相对误差的指标，表示了校正误差与当前参考信号强度的百分比[57]。MAPE的值越小，表示模型的相对误差越小，效果越好。

$\text{MAPE}=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left|\frac{y\left(i\right)-\widehat{y}\left(i\right)}{y\left(i\right)}\right|}$ (28)

3) 皮尔逊相关(R)：是一种统计了两个变量之间线性关系强度的值，它的取值范围在−1到1之间。该值体现了头部运动污染信号的整体相似性，即测量两个信号形状的相似性[58]。皮尔逊相关系数对异常值敏感，比如尖峰伪影等，因此在使用时应注意数据的分布和异常值。

$\begin{array}{c}\rho \left(X,Y\right)=\frac{\mathrm{cov}\left(X,Y\right)}{{\sigma }_x{\sigma }_y}=\frac{E\left[\left(X-EX\right)\left(Y-EY\right)\right]}{{\sigma }_x{\sigma }_y}\\ =\frac{E\left(XY\right)-E\left(X\right)E\left(Y\right)}{\sqrt{E\left(X^2\right)-E^2\left(X\right)}\sqrt{E\left(Y^2\right)-E^2\left(Y\right)}}\end{array}$ (29)

4) 曲线下面积(ΔAUC):是比较两条曲线总体偏离基线的一种全局度量。参考信号和后的信号之间的面积差。如下式(30)，AUC的取值范在0到1之间，值越大，说明了效果越好[59]。

$\Delta \text{AUC}={\text{AUC}}_{\text{HRF}}-{\text{AUC}}_{\widehat{\text{HRF}}}$ (30)

5. 结论与展望

在运动伪影去除技术的研究和应用中，我们可以清晰地看到取得的显著进展和取得的成果。fNIRS在神经科学中占据了一个重要且迅速扩大的空间，提供了巨大的神经成像能力，对参与者行为的限制更少。现实中，fNIRS系统允许对从新生儿到老年人的广泛人群进行神经监测。由此，在fNIRS的使用过程中，由于运动导致的伪影已经是不可避免的。

本文在第二部分系统地阐述了现如今运动伪影去除的主要技术进展，并指出了各项方法的优缺点以及适用范围，帮助读者可以更好地了解该领域的研究成果与发展。第三部分指出了定量比较运动伪影去除方法的定义和公式。

在现有的技术中，每个方法都有自身的特点和可以改进的地方。比如维纳滤波是最优的线性滤波器，但无法处理非线性的信号。样条插值法在去除基线漂移方面有出众的效果，但它依赖于对运动伪影片段的准确检测识别。然而，主成分分析运用正交分解信号来矫正伪影，但有可能存在过度矫正的风险。小波滤波可以显著降低峰值和漂移，但依赖于小波基函数的选择，会很大程度影响滤波的效果。鲁棒局部加权回归可以对信号进行有效地局部平滑，但由于需要每个预测点的权重，因此计算非常复杂。基于相关的信号改进方法也有很强的假设条件，有时候可能并不满足。以上技术虽然是在我们现阶段认为运用较为普遍的fNIRS伪影矫正技术，但这肯定会不断地改进发展。因此越来越多的研究者会使用混合算法来处理运动伪影带来的干扰。

不过，现在该领域的研究也依然存在着一些欠缺，需要更多的研究者来解决。首先，对于复杂运动情况，如非线性运动，现有的运动伪影去除技术可能效果不佳。处理这种复杂运动的情况仍然是一个具有挑战性的问题。其次，当前缺少对伪影去除技术在人的不同状态下稳定性的深入研究，年龄、健康状况等不同，滤波效果是否会存在不同。虽然目前有较多的定量比较指标，但是目前还缺乏一套较为固定的标准评估指标，用于客观评估不同运动伪影去除算法的性能，不同的评价标准会得到不同的比较结果。

期待未来运动伪影去除技术在更多领域得到广泛应用。从实时性要求到复杂运动的处理，从多模态数据的融合到对噪声的鲁棒性，这些挑战将继续激发研究者们寻找创新的解决方案。

NOTES

^*通讯作者。