基于稀疏张量分解的动态功能网络分析

葛诗琪，孟祥辉，张浩威，李靓珠

（长安大学 电子与控制工程学院，陕西 西安 710064）

0 引 言

大脑不同功能的实现需要多个脑区的协同作用，这些脑区在空间上可能相隔较远，但在激活模式上却存在着高度的相关性。功能磁共振成像（Functional Magnetic Resonance Imagining, fMRI）是一种基于血氧水平依赖（BOLD）的成像技术，为研究人脑功能提供了一种重要方法[1]。其中静息状态功能磁共振成像（RS-fMRI）由于能够反映人脑的自发神经活动，已成为研究热点[2]。此外，多种脑功能障碍表现为脑区间功能连接之间的显著变化。故对脑区之间的功能连接（Functional Connectivity, FC）进行研究，有利于促进对大脑活动本质的探索[3]。

传统的基于RS-fMRI的功能网络研究通常假设系统在一段时间之内处于平稳状态。虽然RS-fMRI是受试者在无特定任务的情况下得到的磁共振扫描，然而事实上大脑活动是一个复杂的时变过程[4]。Hutchison等人的研究表明，静息态功能连接网络随时间变化具有波动性[5]。因此研究脑功能网络的动态特性具有重要意义。目前最常用的动态功能连接分析方法是滑动窗口法（Sliding Window Analysis）[6]。Sakoğlu等人利用滑动窗口和独立成分分析发现动态连接与受试者的临床指标相关[7]。滑动窗口方法为研究大脑内部状态变化提供了一种有效的方法，但也存在明显的局限性，窗口选择过大无法准确捕捉功能连接的变化，窗口过小会引入不真实的波动。研究表明，30～60 s的窗口能够提供较为稳健的动态连接网络估计[6]。

本文提出一种基于张量分解的动态功能连接网络分析方法。首先根据静息态fMRI不同维度的信息构建张量模型，并根据不同模态的特点选择约束项进行优化求解。然后利用动态检测法对时间改变点进行检测，构建动态网络。为了验证模型的有效性，在仿真数据和真实数据上将本文模型获得的实验结果与滑动窗口方法获得的结果进行了比较。仿真实验结果表明，本文算法能够克服滑动窗口方法对于窗口选择的限制，显著降低了实验结果的相对误差（最多降低38.4%）。真实数据上的实验结果表明，静息状态下不同脑区之间存在功能连接，主要集中在SM、AUD、DMN、VIS和SN之间。

1 基于张量分解的动态功能连接网络

1.1 张量分解模型

张量（tensor）是对高维数组的总称[8]。一维数组称为向量，二维数组称为矩阵，三维及高维数组称为高阶张量。CP分解是一种常用的张量分解方法，可以将高阶张量分解为一系列秩一张量的线性组合。三阶张量T∈RI×J×K的CP分解可以表示为：

式中：“°”表示张量的外积；R为张量的秩。根据张量的标准矩阵表示形式，可通过下列目标函数对公式（1）求解：

式中： ⊙表示 Khatri-Rao积；||·||F表示Frobenius范数；A=[a1,a2, ...,aR]∈RI×R，B=[b1,b2, ...,bR]∈RJ×R，C=[c1,c2, ...,cR]∈RK×R，为因子矩阵；T(k)表示张量以第k模态展开的矩阵。

1.2 基于静息态fMRI的稀疏张量分解

为了对静息态fMRI数据中不同模态的信息进行融合分析，将所有个体的fMRI数据堆叠为一个三阶张量X∈RI×R×K，其中I表示个体数量，J表示ROI节点个数，K表示时间序列长度。

张量分解完成了对高维数据的降维，分解得到的3个矢量分别表示个体、ROI、时间在有效成分中所占的权重。为了使分解的结果包含更多有意义的信息，需要根据不同维度信息对数据进行必要的约束。对A矩阵添加L2,1正则约束项以避免不同个体之间的相互影响；对B矩阵进行正交约束以保证不同ROI之间的线性无关性；对矩阵C添加L1正则项实现稀疏性。添加约束之后的模型如下：

式中：A∈RI×R，air表示个体i在成分r中所占的权重；B∈RJ×R，bjr表示第j个ROI在成分r中所占的权重；C∈RK×R，ckr表示第k个时间点在成分r中所占的权重。λ1和λ3不仅可以控制重建功能性连接网络的特征提取，而且还确定了有效成分的稀疏性和大小，分解的过程中通过调节惩罚参数λ2实现正交性。

1.3 优化求解

使用交替最小二乘（Alternating Least Squares, ALS）算法对式（3）求解，其主要思想是在每次迭代中固定其中2个因子矩阵来更新另一个因子矩阵，重复该步骤直至算法收敛。可得矩阵A、B和C的更新公式分别为：

式中：sgn(·)为符号函数；Σ∈RR×R是对角矩阵，该矩阵中第r个对角元素为1/||ar||2。当目标函数满足|et+1-et|≤10-8时终止迭代，其中et=||Tt-T||2/||T||2。

2 动态功能连接网络

2.1 时间改变点检测

静息态fMRI数据随时间而波动，故实验的一个重点就是对时间改变点进行检测，并根据时间改变点对状态进行划分。假设时间权重矩阵C可以被划分为τ+1个状态，时间改变点分别位于{t1,t2, ...,tτ}，使用公式（7）对连续2个时间点的距离进行度量[9]。

式中：C(t,k)表示时间点t在第k个有效成分中所占的权重；R表示有效成分的总数。通过式（7）确定d(t)是否为时间改变点。

式中，T为时间点个数，由于距离d(t)服从高斯分布，故选择μ+2σ作为阈值。

2.2 功能连接网络构建

在得到时间改变点之后，提取各状态的主导ROI，计算各状态下的动态功能连接网络。基于稀疏张量CP分解的动态功能连接网络构建主要分为以下步骤。

（1）张量分解：对原始数据构成的张量模型根据式（3）进行分解，得到关于个体、ROI及时间的矩阵A、B、C。

（2）时间改变点检测：根据式（8）基于时间权重矩阵C确定时间改变点{t1,t2, ...,tτ}，并进行状态划分。

（3）提取各状态的主导ROI：2个时间改变点(ti,ti+1)之间为一段稳定的状态。通过矩阵C求出每个稳定状态下的主要成分，提取B矩阵中的对应列，并保留其中权重较大的部分，提取出其中主要的ROI。

（4）构建dFC网络：根据原始数据的连接网络及步骤（3）得到的主导ROI的权重矩阵，建立每个状态下的动态连接网络。

3 实验及分析

3.1 模拟实验

为了评估算法的性能，首先使用模拟数据集进行验证，并与滑动窗口方法进行比较。使用SimTB[10]工具箱进行BOLD时间序列信号生成。产生的信号包括每个受试者的60 ROI，每个ROI包含120个时间点，重复时间为2 s。时间序列被分为4个不同的状态，每个状态包含8～16个不同的ROI，共100个受试者。用二进制表示信号的激活状态，1表示激活，0表示未激活[11]。功能区划分及激活状态如图1所示。

滑动窗口方法中我们测试了两种不同的窗口长度（l=30，40）。本文算法的稀疏性参数λ1、λ2和λ3使用网格搜索法确定，分别为0.9、2和0.1×10-1，R设置为60，此时拟合率相对较高，可达到85.2%。使用真实连通矩阵Cg与算法估计所得连通矩阵CE的相对误差作为评价指标：

图1 四个功能区和四个功能区的激活状态

表1显示了滑动窗口法（SW）和本文方法（STensor）在模拟数据集上的相对误差。可以看出，四种状态下本文算法对连接矩阵的拟合均优于滑动窗口方法，相对误差均明显降低。在状态1中，当滑动窗口的窗口长度取40时，相对误差最多降低了38.4%。

表1 相对误差

3.2 动态功能连接网络分析

本文使用的数据来源于费城神经发育队列（Philadelphia Neurodevelopmental Cohort, PNC）[12]。实验中选取其中878名受试者的RS-fMRI数据，时间序列为124，ROI数量为264。实验中三种稀疏性参数使用网格搜索法进行选择，λ1，λ2和λ3分别取0.1，6×103和0.1×10-3。有效成分个数R取100。在完成稀疏张量分解之后，首先使用动态方法检测得到时间改变点为{3，40，75，90}。由于初始时机器及人体心跳呼吸等噪音造成的波动较大，故舍弃时间改变点3，于是时间序列被分为四段相对稳定的状态[0，40]，[41，75]，[76，90]和[91，120]。找出每个状态的主要成分及主导ROI，计算动态功能连接矩阵，这些ROI可以被分为不同的功能区域。

图2分别为四种状态下每种功能区域之间的连接数，连接数越大表示连接强度越高。从图2中可以看出静息状态下大脑并不是静息的，各个脑区之间存在着许多自发性的神经活动，尤其是视觉网络（VIS）、听觉网络（AUD）、默认网络（DMN）等之间存在着高强度的功能连接，其中DMN的连通性强于其他脑区，这与文献[13]提出的DMN自身具有相对较高的连通性的结论一致。文献[14]指出DMN和VIS之间有着显著的脑行为关联。

图2 dFC网络

图2中4个状态的dFC网络中，较高强度的功能连接均集中在SM、AUD、DMN、VIS和SN几个脑区中，说明大部分功能连接所对应的脑区具有相似性，且在神经解剖学上有重叠部分[2]。

4 结 语

针对静息态fMRI数据的特性，本文提出了一种基于稀疏张量CP分解的动态脑功能网络分析方法。首先利用动态连接检测法对张量分解的结果进行状态划分，该方法打破了滑动窗口方法对于窗口长度选择的限制，显著降低了实验的相对误差。然后对不同状态下的功能连接进行建网分析。实验结果表明，静息状态下各个脑区之间依然存在联系，这些脑区可能在解剖学上相互分离，但由于自发性的神经活动形成了功能连接网络。在静息状态下，功能连接所对应的脑区具有相似性，且大部分集中在SM、AUD、DMN、VIS和SN几个脑区当中，这一结论也在其他多个文献中得到了验证。