近红外光谱技术是一种利用近红外波段(700~1 000 nm)的光来测定大脑局部区域、 氧合、 血流和代谢状态的一种非侵入性的脑功能神经成像技术, 测量由刺激期间相关神经元细胞代谢活动引起的局部毛细血管血流的变化。 当相应脑区处于活跃状态并参与执行某项任务时, 大脑对氧和葡萄糖的代谢需求就会相应的增加, 导致局部脑血流量(cerebral blood flow, CBF)供过于求, 以此来满足大脑增加的代谢需求。 具体表现为: 在任务刺激的作用下, 相应的大脑神经区域脑血流携带的氧会增加, 而氧是通过血液中的血红蛋白进行运输的, 所以在该大脑活动区域血液中的氧合血红蛋白浓度会上升, 而脱氧血红蛋白浓度则会下降, 而我们可以通过测量fNIRS的光衰减变化进行估算, 刺激期间的血液动力学变化如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 刺激期间血液动力学变化 (a): 无激活状态下脑血流; (b)激活后脑血流Fig.1 Hemodynamic changes during stimulation (a): Cerebral blood flow in the absence of stimulation; (b): Stimulation-induced cerebral blood flow

大脑组织活动时, 血氧水平会影响其光学特性。 生理组织对光有两种反应: 一是吸收, 二是散射。 血液的主要成分(水、 氧合血红蛋白和脱氧血红蛋白)对700~900 nm的近红外光的吸收非常小, 具有良好的散射性, 因此700~900 nm的近红外光也被称为光谱窗, 如图2所示。 在这个光谱窗760和850 nm左右的近红外光对氧合血红蛋白(HbO₂)和脱氧血红蛋白(HbR)的吸收率是不同的, 760 nm时脱氧血红蛋白的吸收率更高, 850 nm时氧合血红蛋白的吸收率更高。 由于头部不同皮层内散射的数量不会因为神经活动而变化, 皮层组织散射导致的衰减被认为是恒定的。 因此认知活动中测量到的衰减的变化就近似被认为是吸收的变化, 而这个变化又是大脑活动组织中氧合血红蛋白和脱氧血红蛋白的变化所致。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 HbO2和HbR在近红外段的吸收谱线Fig.2 Absorption spectra of HbO2 and HbR in the near-infrared range

比尔朗伯(Beer-Lamber)定律是1729年由布格(Bouguer)和朗伯(Lambert)提出的, 即: 对于一个均匀吸收粒子分布的媒介, 平行光在媒介中传播的衰减程度依赖于媒介的吸收物质、 光波长与媒介的厚度, 对媒介微元厚度积分, 便可求得透射过的光强衰减。 1852年比尔提出在介质中的吸收物质, 介质吸收光强的能力是吸收物质的浓度与光在介质中的传播的光程长度的指数函数关系, 提出光密度OD(optical density, OD)的计算公式如式(1)

OD=lgI0I=εcl(1)

基本的朗伯比尔定律原则上只适用均一且无散射的媒介, 实际上生物组织体是强散射体, 光的路程并非沿着直线传播, 而是随机改变传播方向。 Delpy等提出微分光路长并总结出散射对光能量损耗的影响, 称为修正的朗伯-比尔定律(modified Lambert-Beer's law, MLBL)表达如式(2)

OD(t, λ)=-log10I(t, λ)I0(t, λ)=∑i=1nεi(λ)ci(t)B(λ)l+G(λ)(2)

式(2)中, n为组织中的多个吸收物质, G为光子背景引起的损耗包括骨骼组织等的吸收, 近红外光照射组织前后其值基本保持不变, 可定为一个常数, B为差分路径因子(differential path-length factor, DPF), 描述的是因散射作用引起的传播路径增长, DPF是光波长的函数, 可由蒙特卡罗模拟或者分辨光谱的实验获得, 在不同组织中DPF也不同, 其值通常大于1, 人体生物组织DPF通常取1^[8]。

在人体脑组织中HbO₂、 HbR是主要的生色团^[9], 人体血流动力学主要检测HbO₂和HbR的浓度变化量, 但是除待测组织外还包括骨骼、 头皮、 皮下脂肪等表层介质, 这些组织将不可避免的对测量产生一些影响。 由上述修正的朗伯比尔定律可知, HbO₂和HbR作为吸收物质, 因为散射衰减因子G的存在, 在脑组织的强散射条件下无法定量确定吸收物质的浓度, 但可从光密度的差量表征HbO₂和HbR的变化量。 初始状态下, 近红外光I入射颅脑表面, 在参考状态下测得透射光I₀, 改变状态下测得透射光I_t, 根据公式可求得参考点和改变状态下的光密度

OD(0)=log10I0I=BLc0+GOD(t)=log10ItI=BLct+G(3)

两个公式相减消去G, 改变状态点相较于参考点浓度该变量为

ΔOD=log10I0It=BL(εHbO2εHb)(ΔcHbO2ΔcHb)(4)

采用双波长测量, 在双波长的光密度差变化量

ΔOD(λ1)=B(λ1)L(εHbO2λ1ΔcHbO2+ελHb1ΔcHb)ΔOD(λ2)=B(λ2)L(ελ2HbO2ΔcHbO2+ελ2HbΔcHb)(5)

式(5)中, εHbO2和ε ^Hb为HbO₂和HbR在λ ₁, λ ₂处的摩尔消光系数, 式(5)可视为一组二元一次方程, 求解可得Δ cHbO2, Δ c_Hb

$\left[\begin{array}{c} \Delta c_{\mathrm{HbO}_{2}} \\ \Delta c_{\mathrm{Hb}} \end{array}\right]=\left(\boldsymbol{E}^{\mathrm{T}} \boldsymbol{E}\right)^{-1} \boldsymbol{E}^{\mathrm{T}} \times \frac{1}{L} \times\left[\begin{array}{c} \Delta O D^{\lambda_{1}} / B^{\lambda_{1}} \\ \Delta O D^{\lambda_{2}} / B^{\lambda_{2}} \end{array}\right]$(6)

其中

E=ελ1HbO2ελ1Hbελ2HbO2ελ2Hb(7)

在生物组织的光学窗口中, HbO₂和HbR的吸收系数在700~900 nm波段处吸收较大, 选用合适的近红外光波长以及光源检测距离便可求得血氧浓度变化量。

本次实验共招募10名健康被试人员(所有被试均为在校学生, 男生8名, 女生2名, 年龄范围为23岁左右)。 被试者双眼矫正视力处于正常范围, 没有任何身体、 精神或者生理方面的问题, 为确保良好的实验效果, 被试者被要求坐在舒服的靠背座椅上, 距离计算机屏幕保持在60~80 cm左右, 实验环境光线偏暗, 以避免环境光的干扰。 所有受试者均同意进行实验, 实验前均已经了解实验范式。

实验使用的是35通道fNIRS成像采集装置(NirSmart, 中国丹阳慧创), 该设备主要采集受试者在进行脑活动实验时大脑的HbO₂和HbR浓度的变化。 设备根据10/20国际标准导联系统进行定位, 由14个光源探头、 14个探测探头构成35个有效通道, 覆盖了受试者双侧前额和运动皮质, 选取的光源和探头(SD)排布方案如图3所示。 光源和探测器用头帽固定, 两者之间的距离为30 mm, 每个通道的中点(即光源和探测器连线的中点)下方脑区为该通道的主要探测区域, 并以此为标准进行通道的脑区标定。 全通道采样频率为11 Hz, 光源探头所用波长分别为730和850 nm, 输出功率大于等于50 mW, 检测探头动态范围大于等于100 dB。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 HbO2和HbR测量SD排布图Fig.3 SD layout of the brain HbO2 and HbR measurements

词语流畅性(Verbal fluency)测验, 是1962年由Thurstone提出, 并用于检测大脑优势半球额叶和颞叶, 主要对被试者的语言表达与沟通、 语义、 记忆和语言理解与执行等功能进行评估。 具体测验的要求为, 让被试者在规定时间内尽可能多地说出实验规定的某一类词语。 通常分为语义流畅性表达、 语音流畅性表示和动作连贯性等。 根据词语流畅性测验的定义以及实验的要求和目的, 为被试者设计并制定实验范式。 具体内容为: 给定一个字母, 要求被试者在规定的时间内尽可能的回忆出以给定字母开头的单词, 例如给定字母“ C” , 被试者此时应该去联想以字母“ C” 开头的单词, 如“ Class” 、 “ Cat” 、 “ Collect” 等。 当电脑屏幕出现“ +” 时, 提醒被试者此时需要处于休息状态, 不进行任何大脑活动想象。 整个实验过程一共设计了 8个任务块, 每个任务块中任务持续时间为10 s, 休息持续时间为20 s。

具体实验范式示意图如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 实验范式示意图Fig.4 Schematic diagram of the experimental paradigm

数据处理流程使用的是Homer2工具箱对fNIRS数据进行预处理, 其中包括数据伪迹识别与校正、 滤波、 计算血氧浓度等, 数据处理流程所用到的函数以及其中参数的设置如表1所示。

表1

Table 1

表1(Table 1)

表1 Homer2工具箱有关函数及参数 Table 1 Functions and parameters of the Homer2 toolbox 函数名 参数名 参数值 hmrIntensity2OD - - hmrMotionArtifactByChannel tMotion tMask STDEVthresh AMPthresh 0.5 1.0 20.0 2.00 hmrMotionCorrectSpline p turnon 0.99 1 hmrBandpassFilt hpf lpf 0.010 0.1 hmrOD2Conc ppf 6.0 6.0 hmrBlockAvg trange -5.0 25

表1 Homer2工具箱有关函数及参数 Table 1 Functions and parameters of the Homer2 toolbox

在表1中, 函数hmrIntensity2OD的作用是将光强度转化为光密度; 函数hmrMotionArtifactByChannel的作用是识别运动伪迹。 以标准差20、 波幅2的标准, 以0.5 s为一个识别窗口进行一个数据的扫描识别。 当检测到数据中的伪迹时, 以该伪迹为开始, 标记1 s时长的数据作为该伪迹的范围。 函数hmrMotionCorrectSpline是Homer2工具箱中提供的一种基于样条插值(spline interpolation)的运动伪影校正函数, 用于对数据进行平滑处理, 以此来校正由运动引起的信号的失真和偏移。 其基本原理是将每个通道的时间序列数据分成多个小段, 对每个小段进行插值运算, 得到运动矫正后的数据。 该函数还可以自动检测运动幅度过大的数据段, 并将这些数据标记为不可靠的数据。 本工作设置的参数如表1所示, 通过观察修正前后的信号变化可知, 样条插值函数有效地改善了数据的质量, 提高了数据的准确性和可靠性; 函数hmrBandpassFilt是滤波函数, 参数0.010~0.1是根据近红外常见噪声所设置的参数值; 函数hmrOD2Conc将光密度转化为血氧浓度, ppf代表的是成人的差分路径因子(老人、 孩子的差分路径因子与成人略有不同); 函数hmrBlockAvg是计算时间窗口为任务开始前5~25 s的血氧浓度的均值叠加。 本文实验对任务状态和静息状态的采样点分别打上数据标签, 在特定的时间窗口内提取特征值。 根据前人的研究和文献参考, 选取了fNIRS常用的特征, 均值、 斜率、 峰峰值以及它们相互组合分别作为分类模型的输入, 以找出分类效果最佳的情况。 分类算法采用的是随机森林模型进行分类识别。

随机森林(random forest, RF)是一种通过集成学习的思想将多棵树集成的一种算法^[11]。 其基本单元是决策树, RF模型实现过程为: 采用随机且有放回的方式随机抽样n个训练样本作为训练子集; 随机地从特征集中选取m个特征子集, 每次树进行分裂时, 根据Gini指数给每次分裂效果量化评分, 不断选择出最佳特征和分割点值作为分裂依据, 直到分裂停止; 之后采用投票策略对N棵决策树的预测结果统计。 伪代码如表2。

表2

Table 2

表2(Table 2)

表2 随机森林算法模型伪代码 Table 2 Pseudocode of random forest algorithm model Algorithm 1: 随机森林模型 1 Input: 训练集D, 特征集A, 决策树个数T Output: 随机森林模型, 由T棵决策树组成 2 begin 3 对于每棵树t=1, 2, …T; 4 从训练集D中随机有放回地抽取n个样本, 构成新的训练D’ (n为样本数量); 5 从特征集A中随机选择m个特征, 构成新的特征集A’ (m为特征数量); 6 构建决策树, 选择最优特征划分节点, 使用Gini系数作为衡量指标; 7 if当前节点样本属于同一类别 then 8 将当前节点标记为叶子节点, 类别为该类; 9 else 10 if 当前节点样本的特征集为空集 then 11 将当前节点标记为叶子节点, 类别为该节点样本中最多的类别; 12 else 13 计算每个特征的Gini系数, 选择Gini系数最小的特征作为最优特征划分节点; 14 end if 15 end if 16 递归地对子节点进行划分, 直接满足停止条件; 17 将构建好的决策树t添加到随机森林模型中; 18 返回随机森林模型; 19 end

表2 随机森林算法模型伪代码 Table 2 Pseudocode of random forest algorithm model

本次实验使用MATLAB编写代码, 主要包括数据训练集测试集的划分, 数据归一化等操作进行模型的构建。 使用70%的数据作为训练集, 而剩下的30%作为测试集, 提取出来的均值、 斜率、 峰峰值以及它们的排列组合作为模型的输入集合训练模型, 其中最大决策树数量设置为60, 模型中的最小叶子数设置为1, 通过袋外数据(out of bag, OOB)对RF分类模型的误差率进行评价, 并将训练20次的输出结果的均值作为分类的最终结果。

涉及前额叶脑区一共有19个通道, 考虑到由于前额中央是脑沟回, 通道信息较弱, 为了使分类效果更加明显, 仅选取左前额叶6个通道和右前额叶6个通道, 共12个通道。 分别提取这12个通道的均值、 斜率、 峰峰值作为模型的输入值进行训练。 同时将被试者进行任务状态的脑活动信号数据设置标签为1, 静息状态设置标签为2, 并对这两种状态进行二分类的识别。 分类结果表明多特征与单特征相较可以提高最终的分类结果。 词语流畅性实验中三特征融合的情况状态分类结果最佳为93.84%, 比单特征均值、 斜率、 峰值分别提高了2.32%、 4.25%、 5.33%。 由实验数据结果可知: 多特征融合在一定程度上可以提高近红外脑功能成像分类的精度。 图5和表3显示了单个特征以及组合特征所对应的运行20次分类结果的平均值以及10名被试者在不同特征下的分类准确率的平均值。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 10名被试者在不同特征下的平均分类准确率Fig.5 The average classification accuracy of 10 participants under different characteristics

表3

Table 3

表3(Table 3)

表3 10名被试者词语流畅性测验的平均准确率/% Table 3 The average accuracy of the word fluency test among 10 participants/% 被试 特征 Mean Slope Peak M+S M+P S+P M+S+P Sub1 90.35 88.94 85.97 90.11 91.30 90.12 91.95 Sub 2 92.85 89.58 94.30 95.11 95.71 92.97 95.59 Sub 3 91.19 92.08 89.58 96.07 93.57 93.33 94.4 Sub 4 90.95 85.83 87.77 90.47 90.38 91.07 92.73 Sub 5 91.21 89.25 88.33 90.35 90.59 91.3 93.09 Sub 6 90.71 92.08 90.13 95.71 93.33 92.73 94.04 Sub 7 90.95 89.09 86.66 92.50 92.26 90.11 93.45 Sub 8 95.11 94.02 86.80 96.90 95.71 94.28 96.54 Sub 9 90.04 89.16 89.72 92.85 93.09 92.50 94.04 Sub 10 91.88 85.83 85.83 93.92 91.90 90.22 92.61 平均值 91.52 89.59 88.51 93.40 92.78 91.86 93.84

表3 10名被试者词语流畅性测验的平均准确率/% Table 3 The average accuracy of the word fluency test among 10 participants/%

由表3和图5可以发现, 单个特征时, 均值作为特征值输入时, 10名被试者的平均分类结果达到91.25%, 斜率作为特征值输入时, 10名被试者的平均分类结果达到89.59%, 峰峰值作为特征值输入时, 10名被试者的平均分类结果达到88.51%, 均值作为单个特征值输入时表现效果最佳。 特征组合时三个特征值融合的分类效果最好, 10名被试者的平均分类结果达到93.84%。 因此, 综合分析可得, 应用特征融合的方法在一定程度上可以提高分类准确率。

以每1 s的时间窗数据作为一个样本, 每个被试者共提取出120个样本特征数据, 每个被试者数据均用于进行随机森林模型的评估。 每个被试者有12个通道, 每个通道均提取3种特征, 共36个特征值。 将36个特征值进行随机森林模型的评估, 选用特征重要性和误差分析作为评估指标。 如图6所示为其中某个被试者的模型评估结果。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 随机森林模型评估指标 (a): 特征重要性; (b): 误差曲线图Fig.6 Random forest model evaluation index (a): Feature importance; (b): Error curve

特征重要性图显示了每个特征对预测结果的影响程度, 由柱状图可以清楚的看出不同特征对于模型性能的重要程度, 其中特征12对模型的影响程度最大, 特征20对模型的贡献程度最小; 误差曲线显示: 随着决策数目的增加, 误差急剧下降并逐渐趋于稳定, 说明模型所设置的决策数目参数60是合适的。