CN104715241B - 一种基于张量分解的fMRI特征提取与识别方法 - Google Patents

Abstract

本发明公开了一种基于张量分解的fMRI特征提取与识别方法。本发明通过构建fMRI动态功能连接矩阵，并利用多重线性主成分分析方法进行张量分解，提取动态功能连接矩阵的特征张量，最后将特征张量向量化后作为支持向量机的输入去识别不同fMRI数据，整个图像序列的预处理及构建矩阵等步骤都通过Matlab程序实现。本发明不仅可以克服“维数灾难”问题，而且能够挖掘出脑功能连接网络的动态信息。

Description

一种基于张量分解的fMRI特征提取与识别方法

技术领域

本发明涉及一种脑功能核磁共振图像(fMRI)的处理与识别方法，是一种基于张量分解的fMRI特征提取与识别方法；具体涉及数字图像处理、动态功能网络构建、张量分解和模式识别等知识领域与处理方法，本方法对fMRI处理的结果可用对不同类fMRI数据，如成人与儿童数据，醉酒与清醒试验者数据，试验者观看不同情感照片数据等进行分类。属于针对高维图像序列数据的机器学习与模式识别方法。

背景技术

功能磁共振图像(fMRI)技术近年来得到研究者的广泛关注。该技术通过探测大脑中的血氧饱和变化来达到无侵入性观察大脑活动变化的目的。其中静息态功能磁共振图像(Resting State fMRI，R-fMRI)反映的是占大脑活动90％能量的自发性活动，这些自发性活动被认为与人类最基础的大脑活动相关，也可以反映儿童与成人随着年龄的变化，大脑认知状态的变化。

近年来，在针对R-fMRI的研究中，研究者一般关注反映不同脑区间神经生理活动相关的功能连接特性。这些特性能反映不同实验参与者脑功能连接上的区别，然而传统的构建功能网络方法通常使用整个时间序列去构建一个静态的网络。该方法虽然易于检验和解释结果，但却忽略了功能网络状态随时间变化这一特性，而这些随时间变化的网络状态往往是反映无约束大脑活动的非稳定特征，正确发现这类特征对理解脑部活动至关重要。本发明利用时间窗构建的动态功能连接特征，不仅能反映功能网络的时间关联特征，而且可以将该这类非稳定特性作为分类不同实验参与者的重要特征。

动态脑功能连接矩阵一般都是三维矩阵数据，将机器识别技术运用在对这类复杂矩阵特征的提取和分析一直是一个难点问题。传统的方法必须将其中二维矩阵向量化后，再使用诸如主成分分析(Principal Component Analysis，PCA)、独立成分分析(Independent Component Analysis，ICA)等一些降维方法进行处理。该类方法在将90个脑区的功能连接矩阵向量化后，会得到一个(90*89)/2维的向量，然而fMRI数据普遍存在被试数目较小，这样不可避免的在识别中会存在“维数灾难”问题，而且PCA与ICA等算法提取的特征向量也难以描述相邻功能连接矩阵之间的动态信息。

发明内容

1.发明要解决的技术问题

针对传统降维方法难以有效提取fMRI动态功能连接矩阵特征这一问题，本发明提供了一种基于张量分解的fMRI特征提取与识别方法，目的在于克服fMRI数据处理中的“维数灾难”问题，并且挖掘脑功能连接网络的动态信息。

2.技术方案

本发明解决问题所采用技术方案的原理是：利用时间窗提取的脑区时间序列获得多组窗时间序列，并构建动态功能连接矩阵。对该矩阵序列运用张量分解的方法降维，提取特征张量，将特征张量向量化后使用支持向量机分类。

本发明技术方案的具体内容，实现步骤如下：

(1)通过SPM8软件处理fMRI原始数据，首先去除前几个非稳定时间点，随后做时间层校正、头动校正、平滑，并在去除低频漂移影响后对图像做00.1～0.1Hz的滤波；最后通过静息态功能磁共振数据处理助手(DPARSF)软件消除全脑信号、脑脊液信号、头动等协变量参数影响；上述操作均在Windows操作系统下完成，实验硬件如下：Intel Xeon E5-26434C3.3010MB 1600CPU-1、Intel Xeon E5-26434C 3.3010MB 1600CPU-2、32GB DDR3-1866(8x4GB)2CPU RAM、2*1TB 7200RPM SATA 1st HDD、AMD FirePro W70004GB 1st GFX、16XSuperMulti DVDRW SATA 1st ODD、HP 3-3-3Warranty；

(2)预处理后的fMRI数据，对照生理学模板(本发明选择Anatomical AutomaticLabeling模版)将除去小脑部分的区域划分成90个功能脑区，并求取单个脑区的平均体素时间序列h_a＝{h_a，1，h_a，2，...h_a，N}1≤a≤90，其中N为时间长度，h_a，N为第a个体素在时间N时的灰度值，于是有全脑90个时间序列H＝{h₁，h₂，....h₉₀}；

(3)利用时间窗思想，将原始长度为N的时间序列，通过设定窗口长度W及相对应的更新步长划分为Z个有重叠的窗口时间序列其中为原全脑时间序列的第k个片段，为第k个片段的单脑区的平均体素时间序列。

(4)通过公式计算脑区窗口时间序列间的皮尔逊相关系数，其中 为同一实验参与者在第k个片段窗口下的第i脑区和第j脑区的时间序列。其中皮尔逊系数常用来表示两个连续变量的线性相关性；反映了第i脑区和第j脑区的功能连接情况，即对应脑区活动的相似程度；共有90个脑区，因此单个窗口的功能连接矩阵为

因为共有Z个有重叠的窗口时间序列，由此得到Z个90*90的功能连接矩阵构成动态功能矩阵。将该矩阵记为一个三阶张量其中模1，2为单个90*90的功能连接矩阵，模3为窗口滑动的时间序列，b指第b个被试；

(5)随机抽取实验数据中的一半数据作为训练集，并将训练集中数据张量记作χ＝ (χ₁，χ₂，...，χ_M)，其中M为训练样本个数。对其计算训练数据张量的均值计 算每个样本的张量与平均张量的差值得到中心化的训练样本张量集

(6)构造训练集的协方差矩阵其中是张量在n(n＝1，2，3)模方向展开的矩阵；对该协方差矩阵进行特征分解，将特征值按大到小进行排序为其所对应的特征向量为使用前P_n个向量初始化特征向量U_n，即由于目标的保留维数{P_n，n＝1，2，3}是未知的；在本发明中，使用一种定义保留比率的方法确定保留维数：该方法首先定义比率其中表示为第i_n个投影到模n方向的特征值，其中I_n为张量n-模的维数。比率Q⁽ⁿ⁾的意义为模n方向保留P_n个特征向量后的保留部分占总张量离散度的比例。在针对本发明数据使用时，在三个模方向上设定Q⁽¹⁾＝Q⁽²⁾＝Q⁽³⁾＝Q，即可在每个模方向上均保留前P_n个特征向量。

(7)确定P_n的值之后，在训练集上使用多重线性主成分分析方法(Multilinear Principal Component Analysis，MPCA)去寻找投影矩阵，该算法即满足其目标函数为其中为张量的总散度，而目标函数的物理意义在于让训练集在张量子空间的距离 最远。

(8)根据目标函数，确定最佳的三个投影矩阵U₁，U₂，U₃；将训练集数据和测试集数据通过三个投影矩阵投影到张量子空间，得到训练样本的特征张量记为γ^tranin＝{γ₁ ^t，γ₂ ^t，...，γ_m ^t}，测试集样本张量为γ^prediction＝{γ₁ ^p，γ₂ ^p，...，γ_m ^p}，其中测试集数据也经过去均值处理；

(9)将全部特征张量按照模1方向展开，向量化后的训练样本特征集记为Z^train＝{z^t ₁，z^t ₂.....z^t _M}，测试样本特征集为Z^prediction＝{z^p ₁，z^p ₂.....z^p _M}，其中z_i为第i个实验数据向量化后的一维向量；将这个特征向量作为支持向量机(SVM)的特征输入，构建识别模型。

3.有益效果

本发明通过时间窗函数将初始的时间序列划分为若干个时间窗序列，这些时间窗序列能够较好的反映大脑的非平稳特性，克服在传统方法下，使用整个时间序列无法反映全部非平稳特征的缺点。

本发明选择用张量表示动态功能连接矩阵，并使用张量分解的方法去提取特征张量。该方法不仅能较为全面的反映动态功能连接矩阵的特性，而且更为重要的是提取的特征具有一定的可解释性。此外从机器学习算法角度出发，使用张量分解不仅解决了“维数灾难”问题，而且提高了识别模型的鲁棒性。

附图说明：

图1为本发明的流程图；

图2为窗口化时间序列的示意图。(图2中N为时间序列长度，W为窗口长度，A_z为单个窗口功能连接矩阵)

图3为模-1方向展开示意图

具体实施例

假设目前有健康成年人30人，10岁以下健康儿童30人，所有被试均是右利手，且通过同一台仪器扫描脑功能图像数据，扫描序列数为224。为测试算法有效性，将成人与儿童被试数据分别取出一半作为训练集数据，余下一半作为测试集。具体算法如下，流程如图1所示。

第一步，对所有数据进行预处理，去除前十个非稳定状态扫描序列，通过SPM8软件对fMRI原始数据做时间层校正、头动校正、平滑，并在去除低频漂移影响后对图像做00.1～0.1Hz的滤波；最后通过DPARSF软件消除全脑信号、脑脊液信号、头动协变量参数影响；

第二步，使用生理学模版(本发明选择Anatomical Automatic Labeling模版)将大脑划分为90个独立的脑区，并对属于同一时间点的各个脑区体素值求取其平均，由此得到单个脑区i的平均体素时间序列h_i＝{h_i，1，h_i，2....h_i，214}，于是有全脑90脑区的平均体素时间序列集H＝{h₁，h₂，....h₉₀}。该步骤目的在于借鉴同一脑区体素具有很强的一致性的先验知识，大大减少特征数据。

第三步，将获取的全脑时间序列窗口化(如图所示)，划分为若干个窗口时间序列。每个窗口长度固定为30个扫描序列，每次滑动步长更新2个扫描序列，最终得到92个互有重叠的窗口时间序列其中1≤k≤92为原全脑时间序列的一个片段，如图2所示。

第四步，利用皮尔逊相关公式构建功能连接矩阵。皮尔逊相关系数体现了在同一窗时间下两脑区活动的一致性，一致性越高，相关系数越高。ρ^k _i，j反映了在第k窗口下，第i脑区和第j脑区的功能连接情况，即对应脑区活动的相似程度；共有90个脑区，因此单个窗口的对角连接矩阵为

因为共有92个有重叠的窗口时间序列，由此得到92个90*90的功能连接矩阵构成动态功能矩阵，将求得的被试的动态功能矩阵记一个三阶张量为

第五步，随机抽取实验数据中的一半数据作为训练集，于是有训练集数据χ^train＝(χ^t ₁，χ^t ₂，...，χ^t ₃₀)，测试集数据χ^prediction＝(χ^p ₁，χ^p ₂，...，χ^p ₃₀)，上述数据集中前15个均为儿童数据，余下为成人数据。对其计算训练数据张量的均值计算每个样本的张量与平均张量的差值得到实验训练集数据和测试集

第六步，在训练集上通过构造训练集的协方差矩阵其中是张量在n(n＝1，2，3)模方向展开的矩阵；对该协方差矩阵进行特征分解，将特征值按大到小进行排序为其所对应的特征向量为因为需要使用前P_n个向量初始化特征向量U_n，即U_n＝{u₁，u₂，...，u_Pn}，由于目标的保留维数{P_n，n＝1，2，3}是未知的；在本发明中，使用一种定义保留比率的方法确定保留维数：该方法首先定义比率其中表示为第i_n个投影到模n方向的特征值，其中I_n为张量n-模的维数。比率Q⁽ⁿ⁾的意义为模n方向保留P_n个特征向量后的保留部分占总张量离散度的比例。在针对本发明数据使用时，在三个模方向上设定Q⁽¹⁾＝Q⁽²⁾＝Q⁽³⁾＝Q，即可在每个模方向上均保留前P_n个特征向量。本发明设定比率为90％，由此在三个模方向上均保留5个特征。

第七步，确定P_n的值为5后，在训练集上使用多重线性主成分分析方法 (Multilinear Principal Component Analysis，MPCA)去寻找投影矩阵，即满足其目标函数为其中为张量的总散度，而目标函数的物理意义在于让训练集在张量子空间的距离 最远，其中MPCA的算法流程如下：

第八步，求得最佳的三个投影矩阵U₁，U₂，U₃；将训练集数据和测试集数据通过三个投影矩阵投影到张量子空间，在经由公式计算γ_m，因为P_n的值为5，所以γ_m为一个5*5*5的张量，由此得到训练样本的特征张量记为γ^train＝{γ₁ ^t，γ₂ ^t，...，γ₃₀ ^t}，测试集样本张量为γ^prediction＝{γ₁ ^p，γ₂ ^p，...，γ₃₀ ^p}；训练集及测试集中均为前15个数据为成人，余下15个为儿童。

第九步，将γ^test按照模1方向展开成向量(如图3所示)，由此每个数据转换为一个1*125的特征向量，记为将训练集特征向量记为Z^train＝{z^t ₁，z^t ₂.....z^t ₃₀}，该训练集中前15个数据为成人，余下15个为儿童。最后将训练集特征向量及对应的特征标签作为支持向量机的输入，训练识别模型。

第十步，根据训练好的模型，使用测试集数据得到根据模型的预测识别标签，根据真实标签分析识别准确率、特异性等性能指标，并分析特征重要度。

本发明首先通过时间窗构建动态功能连接矩阵，该方法克服了传统方法只能获取静态功能连接的缺陷，能够更加全面的反映fMRI数据特点。在将动态功能连接矩阵记为三阶张量后，利用张量分解的方法提取特征张量，完成对动态功能连接矩阵数据的降维。该降维方法相较传统的降维方法，可以很好的保留动态功能连接矩阵前后时间窗间的相关特征，并提供更好的数据分类特征。

Claims (1)

1.一种基于张量分解的fMRI特征提取与识别方法，包括如下步骤：

(1)对获取到的静息态fMRI实验数据进行预处理，消除外界干扰信号并将所有采集到的实验数据标准化到统一的时间和空间域；所述的预处理过程通过SPM8软件实现，包括时间层校准，头动校正，空间标准化和空间平滑；

(2)从预处理后的fMRI数据中提取出相关信息，对照生理学模板，将每个实验者大脑脑区内的所有体素的时间序列进行平均，获得90个脑区的平均体素时间序列；

(3)利用时间窗思想，将单脑区平均体素时间序列划分为多个时间序列，得到多个窗口平均体素时间序列；

(4)计算90脑区在同一窗口下的时间序列间的皮尔逊相关系数，在单一窗口下得到一个功能矩阵，由此将所有窗口下得到的功能矩阵组合，得到动态功能连接矩阵；

(5)随机抽取实验数据中的一半数据作为训练集，并按上述步骤得到动态功能连接矩阵，并将动态功能矩阵作为一个三阶张量，构建张量训练集；

(6)构造训练集的协方差矩阵，设定Q值，确定保留的特征张量维数；

(7)根据特征张量维数，在满足投影的张量总散度最小的情况下，确定最佳投影矩阵；

(8)根据最佳投影矩阵，将训练集及测试集原始张量投影，并求取训练集和测试集的特征张量；

(9)将所有特征张量展开为一维向量；

(10)将该一维向量作为支持向量机的输入，根据训练集构建识别模型，并用测试集检验该模型；

所述步骤(2)中所使用的生理学模版是指生理解剖学模版；步骤(3)中利用时间窗思想，是指将原始长度为N的时间序列，通过设定窗口长度W及相对应的更新步长划分为Z个有重叠的窗口时间序列；步骤(7)中确定最佳投影矩阵的方法为多重线性主成分分析方法；步骤(9)中将特征张量展开为一维向量，指的是将张量按照模-1方向展开成向量。