近年来, 脑研究已经引起了人们的广泛关注.神经影像技术能够以图像形式表示大脑的解剖结构与功能活动, 极大地推动了脑研究的发展.而在非侵入式脑成像技术中, 功能磁共振成像(Functional magnetic resonance imaging, fMRI)因其高空间、时间分辨率被广泛应用于脑功能研究与脑疾病诊断等领域^[1].目前, fMRI主要的成像机制是血氧水平依赖(Blood oxygenation level dependent, BOLD)对比:神经元细胞激活后的供氧需求导致邻近区域内相关的血液动力学状态发生变化, 而BOLD信号测量的则是其中脱氧血红蛋白浓度的变化^[2].并且, 该过程在不同被试甚至同一被试的不同脑区上都会有所区别^[3].这意味着BOLD信号只是神经元活动的一种间接度量, 不能完全等价表示神经元活动.而事实上, 大多数脑功能探究本质上需要的是神经元层面上的直接信息.例如, 脑激活定位是通过神经元活动与实验任务刺激的匹配程度来判定脑区的激活程度^[4], 而脑效应连接则需要通过神经元活动之间的关系来探索它们对应脑区之间的因果关系^[5].因此, 从BOLD信号中还原出真正的神经元活动能够为脑功能研究提供更为准确的大脑活动信息, 该问题也引起了脑功能领域研究者的广泛关注^[6-7].

目前, 大部分还原方法本质上都是对血液动力学模型(Hemodynamic model, HDM)^[8]的反演计算.血液动力学模型是由Balloon模型^[9]派生而来, 用于描述血液动力学过程的数学模型.血液动力学模型使用一组非线性微分方程定量地表达了一系列血液动力学状态之间的关系, 并通过测量函数将相关的血液动力学状态映射成BOLD信号.该模型不仅准确地对血液动力学过程进行了建模, 还有效地解释了BOLD信号中存在的非线性特征^[10].另外, 血液动力学模型中的参数都具有实际的生理意义, 参数值在不同被试和同一被试不同脑区上的变化能够解释血液动力学过程的差异性.血液动力学模型可以被看作一个非线性动态系统, 其输入-状态-输出形式表示如下:

其中下标$t$表示时间, $u$表示系统输入, 即神经元活动, $\pmb{x}$表示一系列血液动力学状态, $\dot{\pmb{x}}$表示其关于时间的导数. $\boldsymbol{\theta}$表示血液动力学模型参数, $y$表示系统输出, 即BOLD信号. $ h $与$g$均为非线性函数.

血液动力学模型使人们对血液动力学过程有了更直观的理解, 并试图利用BOLD信号估计相应的模型参数、血液动力学状态以及神经元活动. 2000年, Friston等首次尝试解决该问题, 他们使用Volterra核函数拟合BOLD信号, 并分析了Volterra核函数与血液动力学模型参数之间的数值关系^[8].随后, 他们又提出一种贝叶斯估计方法, 通过计算血液动力学模型参数在高斯分布假设下的后验分布来对其进行估计^[11]. Riera等提出了一种基于局部线性滤波器(Local linearization filter, LLF)的方法, 对模型的参数与血液动力学状态进行估计^[12].该方法不仅使用了Wiener过程对血液动力学状态的噪声进行建模, 使其具有更强的抗噪能力, 还使用了径向基函数表示神经元活动, 使其能够对任意的神经元活动进行估计.然而, 已有方法大部分都使用线性方法拟合血液动力学过程, 而血液动力学模型中的非线性特征会使它们受到一定限制.针对该问题, Johnston等提出一种基于粒子滤波器(Particle filter)的方法, 该方法能够保留血液动力学模型中的非线性特征, 并对状态以及模型参数进行估计^[13].类似方法还包括基于粒子滤波与平滑(Particle filter and smoothing)的算法^[14], 基于无迹卡尔曼滤波器(Unscented Kalman filter, UKF)的算法^[15].除Riera等^[12]之外, 上述方法在未通过实验设计得知神经元活动的情况下, 都无法同时对神经元活动进行估计.因此, Friston等提出了一种基于变分贝叶斯(Variational Bayesian)的动态系统反演方法, 称为动态期望最大化(Dynamic expectation maximisation, DEM)算法, 该算法不仅能够对模型参数与血液动力学状态进行估计, 还能够估计潜在的神经元活动^[16].随后, Havlicek等提出了一种结合平方根容积卡尔曼滤波器(Square-root cubature Kalman filter, SCKF)与Rauch-Tang-Striebel平滑器的方法SCKS, 用于血液动力学状态、模型参数以及神经元活动的估计^[17]. 2013年, Laleg-Kirati等提出了一种基于观测器(Observer-based)的方法, 并将其应用于血液动力学模型中, 对血液动力学状态以及神经元活动进行估计^[18].类似地, Zayane-Aissa等提出了一种高阶滑模(High order sliding mode, HOSM)观测器, 来通过fMRI数据估计血容积与脱氧血红蛋白含量的变化^[19]. Khorama等提出了一种结合Tikhnov正则牛顿法与CKF的方法TNM-CKF, 在一定程度上解决了SCKS对模型参数先验值的需求^[20]. Aslan等使用粒子平滑期望最大化(Particle smoother expectation maximization, PSEM)算法对血液动力学模型进行反演计算, 以达到对血液动力学状态与模型参数的联合估计^[21].

为了对BOLD信号与神经元活动之间的非线性关系进行有效的建模, 以上方法大都不可避免地需要使用模型参数, 而这些参数的先验信息是无法通过测量得到的, 这就使这些方法在实际应用中受到了一定的限制.人工神经网络是一种强大的监督学习模型, 它能够通过训练自适应地学习数据中的非线性关系, 且已经成功应用在了时序预测^[22]以及模型反演^[23]等问题中. Karam等首次将人工神经网络方法应用至血液动力学状态估计问题中, 提出了一种基于NARX (Nonlinear auto-regressive with exogenous input)网络的血液动力学状态估计方法^[24], 并且不需要任何模型参数的先验信息.然而, 该方法只考虑了BOLD信号与血液动力学状态之间的关系, 而没有考虑血液动力学状态之间逐层递进的非线性关系.另外, BOLD信号是一种时间序列, 它在当前时刻的值包含了前面时刻值的一些信息, 这种时序关系是血液动力学模型内在的一种特性, 而上述大多数非线性方法并没有充分利用这种特性.为此, 本文提出一种基于循环神经网络(Recurrent neural network, RNN)的方法, 能够充分利用血液动力学模型中的时序关系, 对血液动力学状态进行逐层映射与估计.具体来说, 首先利用血液动力学模型中非线性函数的反函数表示BOLD信号与血液动力学状态之间的映射关系, 构成了模型的反演过程.然后, 采用一种新的神经网络结构, 用于拟合上述非线性映射关系, 使其能够通过BOLD信号估计血液动力学状态.该结构包含三个RNN模块, 并通过堆栈方式连接, 称为栈式循环神经网络(Stacked recurrent neural network, SRNN).实验结果表明:与代表算法相比, 新算法能够更准确、快速地对血液动力学状态进行估计, 更具有实际应用价值.

1.   背景知识

1.1   血液动力学模型

血液动力学模型描述了神经元活动、血液动力学状态以及BOLD信号之间的动态关系.该模型首先使用一组微分方程计算血液动力学状态的变化趋势, 然后通过一个非线性函数将血液动力学状态映射成为BOLD信号.血液动力学状态在$ t $时刻的变化趋势计算如下:

其中$ u $表示神经元活动, $ \pmb{x} = \left[s, f_{{\rm in}}, v, q \right] ^{\rm T} $表示归一化的各个血液动力学状态, 包括血脉舒张信号(Vasodilation signal)、血流量(Cerebral blood flow, CBF)、血容积(Cerebral blood volume, CBV)以及脱氧血红蛋白(Deoxyhemoglobin, dHb), $ \dot{\pmb{x}} $表示其关于时刻$ t $的导数.另外, $ \boldsymbol{\theta} = \left[\epsilon, \kappa, \gamma, \tau, \alpha, E_0, V_0 \right]^{\rm T} $表示血液动力学模型涉及到的生理参数, 这些参数的含义与作用将在下面进行详细说明.

具体来说, 当大脑中某个脑区被激活后, 该脑区中神经元活动$ u $的产生会导致血脉舒张信号$ s $的增强, 增强的比例由功效系数$ \epsilon $决定.血脉舒张信号的大小则直接决定了血液流入血管的速度$ \dot{f_{{\rm in}}} $.同时, 血脉舒张信号也会受到其自身以及血流量$ f_{{\rm in}} $的抑制, 分别由对应的衰减系数$ \kappa $与$ \gamma $控制. $ \tau $为常数, 与血液穿过血管的时间相关.伴随着血管中血液的流入与流出, 血容积$ v $也会增大或减小.其中血液流出的速度由血容积与Grubb常数$ \alpha $决定:

血液中脱氧血红蛋白的含量$ q $会在神经元供氧后增加, 同时也会伴随着血液流出而降低, 其中流入的血液中的氧提取率可由下式计算:

其中$ E_0 $表示静息态下氧提取率.

最后, BOLD信号由血容积与脱氧血红蛋白的非线性组合计算得到, 可由函数$ g $表示如下:

其中$ V_0 $为静息态血容积比率, 其他参数如下所示^[9]:

血液动力学模型(2), (5)的连续形式可通过积分计算:

更近一步, 该模型可以转换为离散形式, 表示如下:

其中$ {\pmb f} $为函数向量, 表示一系列状态函数, 由式(2)可知:

1.2   循环神经网络

循环神经网络是人工神经网络的一种, 相较于传统人工神经网络, 其相邻时刻隐层神经元之间的连接使得它能够有效地记忆前面时刻的信息, 已经成功应用于手写体识别^[25-26], 语音识别^[27-28]等领域.一种常见的RNN结构如图 1所示, 由输入层、隐藏层与输出层组成.在$ t $时刻, 隐层状态$ {\pmb h}_t = [h_t^1, h_t^2, \cdots, h_t^{N_H}]^{\rm T} $由输入层$ {\pmb x}_t = [x_t^1, x_t^2, \cdots, x_t^{N_I}]^{\rm T} $与前一时刻隐层状态$ {\pmb h}_{t-1} $共同决定:

图 1  基本循环神经网络结构

Fig. 1  A simple recurrent neural network

下载: 全尺寸图片 幻灯片

其中$ N_H $和$ N_I $分别表示隐层与输入层的神经元节点数, $ U $表示输入层与隐层之间的连接矩阵, $ W $表示相邻时刻隐层之间的连接矩阵, $ \sigma_s $为隐层神经元激活函数, 常见的激活函数有sigmoid函数与tanh函数. RNN的输出由隐藏状态变量计算得到:

其中$ V $表示隐层与输出层之间的连接矩阵, $ \sigma_o $表示输出神经元的激活函数.

激活函数使得RNN能够通过建立输入与输出之间的关系$ {F^{\rm RNN}} $, 拟合任意非线性映射函数:

因此, 本文使用RNN对血液动力学模型中非线性函数的反函数进行拟合.

2.   基于SRNN的血液动力学状态估计

为了利用BOLD信号准确地估计血液动力学状态, 本文提出了一种基于栈式循环神经网络的方法.该方法分为两个步骤, 首先建立BOLD信号至血液动力学状态的非线性映射关系, 然后使用神经网络拟合该映射关系.具体来说, 如式(9)所示, 血液动力学状态与BOLD信号之间的关系可以通过一系列非线性映射函数表示, 通过这些函数的反函数, 便可以建立反向映射关系.针对式(9)中的非线性映射函数$ f_2 $, $ f_4 $以及$ g $, 本文采用一种新的神经网络结构, 使用不同的RNN模块对它们的反函数进行拟合, 并将它们通过栈式堆叠的方式连接起来, 该结构称为SRNN.

2.1   基于反函数的血液动力学模型反演

从BOLD信号中估计血液动力学状态的问题可以通过对血液动力学模型进行反演完成, 即通过求解式(8)的反函数, 可以表示为:

其中$ \tilde{{\pmb x}}_t $表示在$ t $时刻, 通过映射函数与BOLD信号得到的血液动力学状态估计值. $ N $表示时间序列的长度.

因此, 血液动力学状态的估计可以通过求解映射函数的反函数完成.映射函数如式(9)所示, 通过求解其中三个映射函数的反函数, 便可以完成由BOLD信号对所有血液动力学状态的估计, 其中包括测量函数$ g $与两个状态函数$ f_2 $与$ f_4 $.反演过程如图 2所示, 实线部分表示反演过程中的非线性映射关系.接下来对这三个映射函数及其反函数进行介绍.

图 2  血液动力学模型反演过程

Fig. 2  The inversion process of the hemodynamic model

下载: 全尺寸图片 幻灯片

2.1.1   测量函数

测量函数$ {g} $是将血容积$ v $与脱氧血红蛋白含量$ q $映射至BOLD信号的非线性函数, 其反函数$ {g}^{-1} $则是由BOLD信号到血容积与脱氧血红蛋白含量的映射:

然而, 求解二元函数的反函数是一个不适定问题, 通常不存在唯一解, 为了取得令人满意的结果, 需要在求解过程中利用其他信息.

将测量函数$ y = g \left(v, q, \boldsymbol{\theta} \right) $在$ \left(v_{t}, q_{t} \right) $进行泰勒展开:

其中$ \frac{\partial ^{n_1 + n_2} g}{\partial v^{n_1} \cdot \partial q^{n_2}} $表示函数$ g $对$ v $和$ q $的$ n_1, n_2 $阶偏导数.由测量函数的泰勒展开可知, 当前时刻的BOLD信号包含了当前时刻与前一时刻血容积与脱氧血红蛋白含量的信息.因此, 利用这些信息, 我们可以建立映射关系$ \bar{g}^{-1} $:

其中, $ v_t $、$ q_t $、$ v_{t-1} $与$ q_{t-1} $分别表示当前时刻与前一时刻$ v $和$ q $的值.

因此, 根据映射函数$ \bar{g}^{-1} $, 我们可以根据当前时刻的BOLD信号与前一时刻血容积与脱氧血红蛋白含量的估计值, 计算当前时刻的血容积与脱氧血红蛋白含量值.

2.1.2   状态函数

如式(9)所示, 状态函数$ f_3 $与$ { f}_4 $的反函数分别表示如下:

事实上, 两者在本质上并没有区别, 这是由于在使用神经网络对这些映射关系进行拟合时, 为了充分利用血液动力学模型中的时序特征, 我们使用上一个模块的隐含层作为当前模块的输入.隐含层作为输出层的一种特征表示形式, 包含了输出层的所有信息.这里使用$ f_4 $的反函数表示对血流量的估计, 因此, 根据映射函数$ f_4^{-1} $可知, 我们可以使用血容积与脱氧血红蛋白含量的估计结果$ \tilde{v} $与$ \tilde{q} $对血流量$ f_{{\rm in}} $进行估计.

需要注意的是上述过程中涉及到血液动力学状态在时间上的不同.其中反函数(14)的输出是$ t $时刻的$ \left[v_t, q_t \right] $, 而在式(15)中, 通过输入$ t $时刻的结果后, 可以得到$ t-1 $时刻的$ f_{{\rm in}} $估计值.

同理, 通过状态函数$ f_2 $的反函数, $ t-1 $时刻的$ f_{{\rm in}} $估计值将被用于对$ t-2 $时刻的$ s $进行估计:

至此, 通过使用血液动力学模型中映射函数的反函数, 我们构建了从BOLD信号到所有血液动力学状态的映射关系, 从而建立了模型的反演过程, 接下来我们将使用神经网络结构对该过程进行拟合.

2.2   基于栈式循环神经网络的血液动力学模型反演

为了拟合上述反演过程, 本节将介绍一种新的神经网络结构.如图 3所示, 该结构以RNN为基础模块, 以时间为单元进行展开.而BOLD信号与血液动力学状态也是时间序列, 因此, 该结构每个时刻的输入为BOLD信号在各时刻的值$ y_1, y_2, \cdots, y_N $.输出为各血液动力学状态的估计值$ {\tilde{\pmb x}}_1, {\tilde{\pmb x}}_2, \cdots, {\tilde{\pmb x}}_N $.

图 3  栈式循环神经网络结构

Fig. 3  An illustration of the proposed SRNN

下载: 全尺寸图片 幻灯片

该结构由三个RNN模块组成, 如结构图 3右侧所示, 每个模块可以表示为一个如式(10)所示的非线性映射, 用于拟合反演过程中涉及到的三个反函数(14) $ \sim $ (16).

首先, 模块1通过拟合式(14), 建立由BOLD信号到血容积与脱氧血红蛋白含量的非线性映射, 完成对这两个血液动力学状态的估计.即模块1的输入为BOLD信号$ y $, 输出为血容积$ v $与脱氧血红蛋白含量$ q $.在时刻$ t $, 模块1的隐层状态计算如下:

模块1的输出如下所示:

其中, 上标$ (1) $表示该元素属于模块1, $ \tilde{v_t} $与$ \tilde{q_t} $为由模块1得到的血容积与脱氧血红蛋白含量的估计值.

然后, 如式(15)所示, 血流量的计算需要的是血容积与脱氧血红蛋白含量的估计值, 即模块1的输出.然而, 从本质上来说, 单隐层神经网络结构中隐层节点表示的是输入层与输出层的一种编码形式, 它包含了输出层的信息, 因此能够用于代替输出层作为模块2的输入, 即通过栈式堆叠连接两个模块, 具体计算方式如下:

血流量$ {f_{{\rm in}}}_{t-1} $的估计值可以通过下式得到:

与模块2类似, 根据式(16), 模块3对血液动力学状态$ s_{t-2} $的估计可表示为:

最后, 通过最小化对应的均方差(Mean square error, MSE), 分别完成三个模块的训练.

其中$ \tilde{v}_t $, $ \tilde{q}_t $, $ \tilde{f_{{\rm in}}}_t $和$ {\tilde{s}}_t $分别表示各RNN模块的输出.在训练神经网络的过程中, 本文使用的是BPTT (Backpropagation through time)算法, 该算法通过将误差反向传播到前面的时刻, 根据梯度下降的原则调整每个模块中的网络权重.

当神经网络完成训练后, 它便可以通过输入BOLD时间序列, 得到各个时刻血液动力学状态的估计值.

3.   实验结果与讨论

为了对新方法的性能进行验证和分析, 本文选取了两种代表性的滤波器算法与一种神经网络方法进行比较, 分别是DEM^[16]、SCKS^[17]与NARX网络^[24].本章首先介绍了数据的生成, 然后介绍了SRNN的训练过程, 最后从各方面对比了本文算法与代表算法的性能.

3.1   仿真数据

在真实fMRI数据中, 其对应的神经元活动与血液动力学状态是不可测量的.因此, 通过生成仿真数据来训练SRNN以及验证其有效性是很有必要的.同时, 为了保证SRNN能够在真实fMRI数据中具有实际应用价值, 需要使仿真数据能够尽可能真实地反映实际fMRI数据的特性.具体来说, 生成仿真数据主要包含以下两个步骤:首先生成神经元活动时间序列, 然后将神经元活动作为血液动力学模型的输入, 并对模型进行求解, 最后得到对应的血液动力学状态以及BOLD信号.因此, 本文从神经元活动与血液动力学模型两个部分来介绍如何在生成仿真数据的同时保持其对真实fMRI数据的还原程度.

3.1.1   神经元活动

对于血液动力学模型来说, 神经元活动$ u(t) $是其唯一的外部输入, 在血液动力学模型参数保持不变以及不考虑系统振动的情况下, 神经元活动序列对血液动力学状态以及BOLD信号的变化起到了关键性的作用.因此, 根据fMRI实验范式来生成对应的神经元活动是仿真数据拟合真实fMRI数据的关键步骤之一.目前fMRI实验中, 神经元激活模式主要包括组块设计模式(Blocked-design)与事件相关(Event-related)模式, 分别表示对神经元的持续刺激与瞬时刺激.其中, 事件相关模式具有更高的灵活性, 能够更真实地反映大脑的活动模式, 因此, 本文根据事件相关范式生成神经元活动.

我们使用高斯函数来表示事件相关模式中神经元的激活^[17].首先介绍单个刺激的情况, 假设神经元在时刻$ t = a $受到刺激, 其响应$ u(t) $可表示为如下形式:

其中$ m(t) $为高斯函数, 其中$ \sigma $是用于调整最终生成的BOLD信号的强度的参数, 能使其与实际采集到的信号强度相近, 根据实际fMRI信号, 本文令$ \sigma = 8 $. $ v(t)*m(t) $表示$ v(t) $与$ m(t) $的卷积, 用于描述神经元对刺激的响应与恢复静息态的过程.

为了表示更复杂的情况, 我们采用如下方式定义神经元在受到多个刺激下的时间序列:

其中$ N_a $表示刺激的个数, $ N_a \sim {\rm U}(3, 5) $并对其取整. $ a_i $表示每个刺激出现的时间, $ a_1, a_2, \cdots, a_{N_a} $相互独立且有$ a_i \sim {\rm U}(0, T) $, $ T $为时间序列长度, 在本文中$ T = 64 $. $ \varepsilon_i $表示每个刺激的强度, $ \varepsilon_1, \varepsilon_2, \cdots, \varepsilon_{N_u} $相互独立且有$ \varepsilon_i \sim {\rm U}(0, 1) $.通过对上述参数$ N_a, a_i, \sigma_i $进行随机采样, 便可以得到相互独立且同分布的神经元时间序列.例如, 随机确定$ N_u = 4 $, 它们的位置为[7; 25; 34; 56], 大小为[1; 0.7; 0.9; 0.2], 卷积后的神经元活动如图 4所示.

图 4  随机生成的神经元活动

Fig. 4  The random generated neural activity

下载: 全尺寸图片 幻灯片

通过重复采样, 便可得到各种刺激模式下神经元响应的时间序列.

3.1.2   血液动力学模型

在生成了神经元活动时间序列后, 需要选择合适的模型来生成对应的血液动力学状态以及BOLD信号, 用于SRNN的训练与验证.血液动力学模型是描述生成过程的生理物理机制的模型, 其有效性是保证血液动力学状态及BOLD信号与真实fMRI数据的匹配程度的另一个关键因素.血液动力学模型的具体原理已经在第1.1节中进行了详细的介绍, 事实上, 其对真实生理物理过程建模的有效性已经在文献[8]中进行了验证, 并且该模型已经在血液动力学系统反演的问题中受到了广泛的应用^[11-18].因此, 本文通过求解以神经元活动作为外部输入的血液动力学模型, 来生成仿真数据.另一方面, 在实际情况中, 一些血液动力学状态的值应该始终为正(如血流量, 血容积与脱氧血红蛋白含量)^[16], 所以, 在生成数据的过程中, 我们对相应的血液动力学状态进行如下转换:

因此, 血液动力学模型便转换为:

根据式(25), 血液动力学模型可以分为两个部分:状态转换函数与观测函数.其中状态转换函数可以看成一个一阶非线性微分方程组, 通过给定血液动力学状态初始值与外部输入(神经元活动)并对其进行求解, 便可得到血液动力学状态时间序列.模型的求解使用局部线性(Local linearization, LL)方法^[29], 在求解血液动力学模型的过程中, 模型参数设置如表 1所示, 最后得到的血液动力学状态如图 5 (a)所示, 由于各状态采用不同单位, 因此表示为arbitrary unit (a.u.).

图 5  随机生成的血液动力学状态及BOLD信号

Fig. 5  The random generated hemodynamic states and the BOLD signal

下载: 全尺寸图片 幻灯片

表 1  血液动力学模型参数默认值

Table 1  The default value of hemodynamic model parameters

Parameters	Description	Default
$ \epsilon $	Neural efficiency	0.50
$ \kappa $	Rate of signal decay	0.65
$ \gamma $	Rate of flow dependent elimination	0.41
$ \tau $	Hemodynamic transit time	0.98
$ \alpha $	Grubb's exponent	0.32
$ E_0 $	Resting oxygen extraction fraction	0.34
$ V_0 $	Resting blood volume fraction	0.08

下载: 导出CSV 

| 显示表格

在得到了血液动力学状态时间序列后, 通过观测函数计算对应的BOLD信号.然而, 在实际fMRI数据采集过程中, 往往存在许多噪声, 如热效应、磁场扰动等.因此, 根据对真实数据中噪声的估计, 在BOLD信号中加入了方差为0.0025的噪声^[12].

最后, 我们便得到了仿真数据集$ u^{(k)}(t) $, $ x_i^{(k)}(t) $, $ y^{(k)}(t) $, 其中$ t = 1, 2, \cdots, 64 $表示时刻, $ i = 1, $ 2, 3, 4分别对应各个血液动力学状态变量, $ k = 1, 2, \cdots, $ 10 000表示生成了10 000组样本, 其中6 000组作为训练集, 用于SRNN的训练, 2 000组作为验证集, 用于对SRNN的参数调整, 2 000组作为测试集, 用于对SRNN的性能进行测试, 并与各代表算法进行比较.

3.2   训练神经网络

我们使用谷歌开源框架Tensorflow对SRNN进行构建与训练, 训练过程中的参数设置为:训练批量大小为256, 学习率为0.1, 衰减系数为0.96, 每迭代200次进行1次衰减, 训练过程中使用Dropout技术防止过拟合.训练过程包括对三个模块的单独训练与结构调整, 其中每个模块包括参数调整与模型训练两个步骤: 1)使用训练集进行训练, 使用验证集误差验证不同网络结构的有效性, 选择最合适的模型.其中, 模块1的隐层节点数调整区间为15至30个, 模块2与3的隐层节点数调整区间为10至20个, 在保证精度接近的情况下尽可能地使用更少的节点数, 以得到更简单的模型. 2)使用训练集对选定的模型进行训练, 直到验证集误差小于一定阈值或不再下降, 则认为该模块已经收敛.在完成前一个模块的参数调整与训练后, 再使用前一个模块的训练结果根据第2.2中的规则进行下一个模块的训练.

通过验证集对模型进行调整后, 我们得到了如下网络结构:模块1的隐层节点为25, 模块2与模块3的隐层节点数为15, 并对其进行训练直至收敛.最后, 使用训练得到的神经网络进行测试, 并与对比算法进行比较.对比算法DEM、SCKS的参数与SPM12中DEM工具箱函数(DEM\_demo\_hdm\_SCK.m)相同, NARX网络的结构与参数与原文相同^[24].

3.3   评估方式

为了对不同算法进行比较, 我们使用平方误差损失(Squared error loss, SEL)函数对算法性能进行评估, 具体如下:

其中$ K $表示测试集样本数量. $ \tilde{x}_i^{(k)} (t) $表示模型对第$ k $个样本中第$ i $个血液动力学状态在$ t $时刻的估计值.我们使用平方误差损失的对数形式lg(SEL)展示各算法的误差, lg(SEL)的值越小则表示算法的结果越准确.

3.4   实验结果

3.4.1   时间特性提取

在本节中, 我们测试不同的循环神经网络结构的非线性拟合能力以及对时间特性的提取能力, 包括如图 1所示的基本RNN结构与长短期记忆(Long-short term memory, LSTM)网络^[30]. RNN与LSTM都属于循环神经网络, 相较于传统神经网络, 它们的优势在于能够将之前的信息应用于当前时刻的任务中.然而, RNN与LSTM的区别在于, 由于循环神经网络是以时间为单位进行展开, 误差的反向传播由当前时刻传递至前一时刻, 梯度消失现象使得RNN无法学习到多个时刻之前的信息. LSTM则通过刻意的结构设计来避免这一问题, 使其能够在不需要付出额外的代价的情况下, 自动地学习到多个时刻之前的信息.

我们分别以RNN与LSTM结构为SRNN中的基础单元, 并使用生成的仿真数据训练集对两种SRNN进行训练.在训练过程中, 两种SRNN对血容积$ v(t) $估计结果的lg(SEL)随迭代过程的变化如图 6 (a)所示, 收敛后结果如图 6 (b)所示.

图 6  RNN模块与LSTM模块结果对比

Fig. 6  The comparison of RNN module and the LSTM module

下载: 全尺寸图片 幻灯片

从收敛后的状态估计结果也可以看出, RNN得到的结果与真实值之间的误差比LSTM更大.该现象可能由以下原因导致: 1)尽管两者使用相同的网络参数, 但是由于LSTM在隐层的操作更为复杂, 使其具有比RNN更好地非线性拟合能力, 能更好地提取血液动力学模型中的非线性特征. 2)由血液动力学模型可知, 血液动力学状态在前后时刻存在依赖关系, BOLD信号作为血液动力学状态的一种表现形式, 其各个时刻之间的值也存在相应的时序关系.事实上, BOLD响应往往在神经元被激活的数秒后到达峰值, 并且神经元的影响还会持续很长的一段时间. RNN无法很好地处理这种长期依赖关系, 而LSTM是针对这种长期依赖关系刻意设计的结构, 因此能够取得更好的效果.可以看出, 针对血液动力学模型的反演, LSTM比RNN具有更好的非线性拟合能力, 并且能够更好地提取BOLD信号中存在的时间特性.因此, 在接下来的实验中, 我们使用LSTM作为基础单元, 训练并测试SRNN的性能, 并与代表算法进行对比.

3.4.2   血液动力学状态估计

我们使用训练集对SRNN进行训练, 并用测试集对比各算法的性能.图 7给出了SRNN与代表算法对各血液动力学状态在100个测试集样本上的平均结果, 在这里使用-lg(SEL)表示各算法的结果.如图 8所示, 为了更直观地展示各算法的结果, 我们在测试集的结果中, 选取了一组最接近平均值的结果.其中图 8 (a)、8 (b)、8 (c)、8 (d)分别是各算法对血脉舒张信号、血流量、血容积以及脱氧血红蛋白含量的估计.图 8 (e)表示由测量函数以及SRNN对血容积与脱氧血红蛋白含量的估计值计算得到的BOLD信号及其误差.

图 7  各算法对比结果

Fig. 7  The result of different approaches

下载: 全尺寸图片 幻灯片

图 8  DEM、SCKS以及SRNN的实验结果对比

Fig. 8  The results of different approaches for estimation of hemodynamic states

下载: 全尺寸图片 幻灯片

在对血脉舒张信号$ s(t) $, 血流量$ f_{{\rm in}}(t) $, 血容积$ v(t) $以及脱氧血红蛋白含量$ q(t) $的估计中, SRNN均取得了最好的效果, 这表明了SRNN能够动态地拟合血液动力学模型中的非线性特征, 而LSTM特殊的网络结构也能够很好地提取BOLD时间序列中存在的时间特性, 更准确地估计血液动力学状态.

从NARX网络的结果可以看出, 其对各个血液动力学状态的估计结果随着模型反演过程的递进而呈现出降低的趋势, 这是由NARX网络结构特性导致的.具体来说, 它使用同一种网络用于预测不同的血液动力学状态, 输入均为BOLD信号.然而, 由血液动力学模型可知, 从BOLD信号到不同血液动力学状态之间的非线性程度存在一定差异, 这就导致了NARX对血脉舒张信号与血流量的估计相对较差. SRNN通过栈式堆叠各个模块, 用于拟合各个血液动力学状态内部的非线性关系.从实验结果可以看出, 其在血脉舒张信号与血流量上也取得了较好的效果, 有效地解决了这个问题.

另外, 虽然SRNN在血液动力学状态的估计上取得了优于代表算法的结果, 但也存在一定的误差.由图 8可知, 对于所有的血液动力学状态估计, 在前几个时刻的结果往往并不理想, 这是由于在训练以及测试时, 网络的隐层状态均初始化为零, 导致SRNN的隐层状态在前几个时刻不能得到充分的信息.而随着时间的推移, SRNN中隐层状态得到了足够的信息, 神经网络也趋于稳定.这也说明了新方法能够提取血液动力学模型中的时序特征, 有助于对血液动力学状态的估计.

尽管本文并不是对神经元活动进行估计, 但准确地对相应的血液动力学状态进行估计能够为神经元活动估计提供帮助.由状态转换函数$ f_1 $可知, 对神经元活动的估计可以由血脉舒张信号与血流量得到, 而新方法对这两个血液动力学状态的估计优于代表算法, 因此, 新方法能够得到更准确的接近神经元活动层面的信息, 为脑功能研究提供更准确的大脑活动信息, 进一步推动脑认知的发展.

3.4.3   时间性能比较

在本节中, 我们将对各个算法的时间性能进行比较, 并讨论它们在实际情况下的时间成本与应用价值.

首先需要说明的是, SRNN与NARX算法属于监督学习算法, 需要使用大量的数据进行训练, 使神经网络能够从数据中学习出其内在的映射规律.而DEM与SCKS则不需要进行预训练, 可以直接输入BOLD信号, 对血液动力学状态进行估计.因此, 将神经网络在训练上的时间复杂度、它们的训练时间, 以及所有算法在测试集样本上运行的时间进行了比较, 如表 2所示.其中, $ K_D, K_{SC}, K_N, K_{SR} $分别表示各算法的迭代次数, $ T $表示BOLD时间序列的长度, $ n_s, n_p $分别表示算法需要估计的血液动力学状态与模型参数的数量, 本文中$ n_s = 4, n_p = 7 $.在神经网络训练时间复杂度中, $ m $表示训练批量大小, $ n_{hn} $与$ n_{hs} $表示神经网络的最大隐层数.

表 2  时间性能比较

Table 2  The comparison of time performance

方法	训练复杂度	运算复杂度	训练时间(s)	运算时间(s)
DEM[16]	$ - $	$ {\rm O}(K_D T n_s^3) $	$ - $	$ 16.3 $
SCKS[17]	$ - $	$ {\rm O}(2K_{SC} T n_p^3) $	$ - $	$ 35.7 $
NARX[24]	$ {\rm O}(K_N T m n_{hn}) $	$ {\rm O}(K_N T n_{hn}) $	11 367	0.09
SRNN	$ {\rm O}(4K_{SR} T m n_{hs}) $	$ {\rm O}(4K_{SR} T n_{hs}) $	1 831	0.2

下载: 导出CSV 

| 显示表格

由表 2中结果可以看出, SRNN与NARX网络均需要花费一定的时间对网络进行训练, 其中SRNN能够以更快的速度收敛.另一方面, 由各算法在单个样本上运行所需要的平均时间可看出, 神经网络方法具有较大的优势.这是由于对于每一个BOLD时间序列, DEM与SCKS方法都需要对模型参数以及血液动力学状态进行迭代估计, 每次迭代的时间与时间序列的长度以及需要估计的模型参数数量有关.而神经网络方法在经过训练后, 将BOLD信号输入训练好的模型, 进行网络的前向计算便可得到其对应的血液动力学状态估计, 其运算时间则取决于神经网络的层数、每层的隐层节点数以及输入的时间序列长度等参数.

事实上, fMRI数据是典型的高维数据, 一个被试的全脑中包含的BOLD时间序列个数可能高达几十万个, 经过挑选后用于功能分析的数量也仍会是一个很大的数目.因此, 当样本数量很大时, 尽管神经网络方法需要一定的训练成本, 但是其在总体时间性能上仍然优于传统算法.另外, 相较于NARX网络, SRNN能够更快地收敛, 因此, 本文算法在时间性能上优于其他对比算法, 具有较高的实际应用价值.

3.4.4   模型参数先验对代表算法的影响

本文提出的方法不需要对血液动力学模型中的模型参数进行确定, 而代表算法则会受到一定的影响.为了探究模型参数的先验值对代表算法的限制, 我们使用不同的模型参数先验值来测试代表算法的性能.以$ \boldsymbol{\theta} $表示模型参数的真值, 以$ \boldsymbol{\theta}_{D} $与$ \boldsymbol{\theta}_{S} $分别表示代表算法DEM与SCKS中模型参数的先验值.我们对$ \boldsymbol{\theta}_{D} $与$ \boldsymbol{\theta}_{S} $在区间$ \left[0.8\boldsymbol{\theta}, 1.2\boldsymbol{\theta} \right] $上均匀取100个值, 在其他参数相同的情况下运行算法, 对各血液动力学状态估计的结果如图 9所示.从结果可以看出, 两个代表算法均会受到其先验信息的影响, 算法的参数先验值与参数的真实值的偏差越大, 算法的误差则越大.事实上, 不同的被试或者同一被试不同脑区上的参数都会有所区别, 这使得代表算法在实际应用中受到了一定的限制.另一方面, 由于本文算法直接通过神经网络对BOLD信号与血液动力学状态之间的非线性关系进行拟合, 而不需要对血液动力学模型参数进行估计, 因此能够避免先验信息带来的限制, 更具有实际应用的价值.

图 9  血液动力学模型参数先验值对算法DEM与SCKS性能的影响

Fig. 9  The influence of prior knowledge on model parameters on DEM and SCKS

下载: 全尺寸图片 幻灯片

4.   结论

本文提出了一种基于栈式循环神经网络的方法, 使用BOLD信号对血液动力学状态进行估计.具体来说, 首先根据血液动力学模型中映射函数的反函数建立模型的反演过程, 然后使用循环神经网络对反函数进行拟合, 并通过栈式堆叠的方式连接各RNN模块, 完成血液动力学状态的估计.循环神经网络结构能够充分利用血液动力学模型中的时序关系, 有效地估计血液动力学状态.另外, 作为人工神经网络的一种, RNN具有强大的非线性拟合能力, 能够学习BOLD信号与血液动力学状态之间的非线性关系.实验结果表明:与代表算法相比, 新方法能够更快速、准确地估计血液动力学状态.另外, 由于新方法不需要对血液动力学模型参数设置先验值, 因此具有更高的实际应用价值.