CN113963349A - 一种提取个体空时特征矢量与被试细分类的方法 - Google Patents

Abstract

一种提取个体空时特征矢量与被试细分类的方法，属于生物医学信号处理领域。其充分挖掘“空间体素×时间×被试”形式多被试fMRI数据经Tucker分解所获取核张量中包含的高维耦合关系，针对空间稀疏约束Tucker分解方法所获取的核张量，提供一种各被试特有空时特征矢量的提取方法，并将这些空时特征矢量用于k‑means被试细分类。在10个健康被试任务态fMRI数据的个体空时特征矢量提取与被试细分类中，根据DMN成分对应的空间特征矩阵，将所有被试分为两组，组1各被试空间激活中IPL区域平均激活体素数比组2多109％；根据任务相关成分所对应的时间特征矩阵，将所有被试分为两组，组1任务态时间过程与参考成分的平均相关系数比组2高55.6％。这些个体空时差异能够为脑功能研究和脑疾病诊断提供新的客观依据。

Description

一种提取个体空时特征矢量与被试细分类的方法

技术领域

本发明属于生物医学信号处理领域，涉及到一种提取个体空时特征矢量与被试细分类的方法，具体是指从多被试功能磁共振成像(functional magnetic resonanceimaging，fMRI)数据的空间稀疏约束Tucker分解核张量中提取个体空时特征矢量与被试细分类的方法。

背景技术

fMRI广泛应用于脑功能和神经精神类脑疾病研究。其优点在于高安全性、非侵入式以及毫米级高空间分辨率。多被试fMRI数据有5维，包括3维全脑数据、1维全脑扫描次数(即时间点个数)和1维被试个数。在进行盲源分离(blind source separation,BSS)时，通常将全脑数据展开成一维体素，此时的多被试fMRI数据“空间体素×时间×被试”也高达3维。

张量分解方法因其充分利用数据高维结构信息的优点，非常适合多被试fMRI数据的分析。Tucker分解是张量分解方法的一种，对于体现多被试fMRI数据高维空时结构的“空间体素×时间×被试”形式张量，既能分解出多被试共享的空间激活成分(spatial maps,SMs)和共享的时间过程成分(timecourses,TCs)，还能分解出蕴含各被试丰富个体信息的核张量。实际上，被试个体间存在认知功能、认知策略、年龄、性别等多方面的差异，在共享的空间激活成分之外，被试间存在着广泛的空间差异性和时间差异性。这些个体空时差异是对被试进行细分类的重要依据。与现有独立矢量分析(independent vector analysis，IVA)算法如IVA-GL(M.Anderson,T.Adali,X.L.Li,“Joint blind source separationwith multivariate Gaussian model:algorithms and performance analysis,”IEEETransactions on Signal Processing,vol.60,no.4,pp.1672-1683,2012)和组独立成分分析(groupindependent component analysis，GICA)方法(V.D.Calhoun,T.Adal₁,et al,“A method for making group inferences from functional MRI data usingindependent component analysis,”Human Brain Mapping,vol.14,no.3,pp.140-151,2001)所提取的单被试特有SM和单被试特有TC相比，Tucker分解核张量所蕴含的个体空间特征矢量和个体时间特征矢量与多被试共享SMs和共享TCs同时获取，完整代表了多被试共享SMs和共享TCs之外的个体空间差异性和时间差异性。而且，因为Tucker分解利用了多被试fMRI数据的高维结构，核张量中的个体空间特征矢量和个体时间特征矢量包括各被试特有SMs成分和特有TCs成分的高度压缩元素，当基于k-means算法对被试进行细分类时，能区分各被试之间空间激活和时间过程的形态差异。

然而，Tucker分解核张量目前很少被利用。目前仅限于在“1维/2维连接×时间×被试”形式(即3维/4维张量)中用于fMRI脑功能连接分析。对于“空间体素×时间×被试”形式的多被试fMRI数据，还没有对Tucker分解核张量进行研究利用的方法。

发明内容

本发明充分挖掘“空间体素×时间×被试”形式多被试fMRI数据经Tucker分解所获取核张量中包含的高维耦合关系，针对性能最优的空间稀疏约束Tucker分解方法所获取的核张量，提供一种提取各被试特有空间特征矢量和时间特征矢量的方法，并将这些空、时特征矢量用于k-means被试细分类，为个体空时差异性研究提供更为有效的方法。

本发明采用的技术方案如下：

首先，在Tucker分解模型中引入空间稀疏约束，形成如下模型：

其中，

是多被试fMRI数据，V是脑内体素的个数，T是时间点个数，K是被试个数；

是共享SM矩阵，

是共享TC矩阵，

是核张量，

是残差张量，N是共享成分个数；“×₁”和“×₂”为模-1乘积和模-2乘积；“||·||_F”、“||·||₁”、“||·||_p”分别为l_F范数、l₁范数和l_p范数(p为稀疏参数，且0<p≤1)；δ、λ、γ分别是空间稀疏项、核张量稀疏项和残差张量稀疏项参数，取值都在(0,1]范围内；式(1)中，空间稀疏约束通过S的l_p范数实现，S和B的低秩约束由l_F范数实现，G和E的稀疏约束通过l₁范数实现。

由式(1)，得到增广拉格朗日函数如下：

式中，

是G的分裂变量，

是拉格朗日乘子，V是脑内体素的个数，T是时间点个数，K是被试个数，N是共享成分个数；α、β是惩罚参数，“<·>”是矩阵内积；

和

分别为张量R、X、E、U、W和G的第k个正面切片，分别满足R _k＝R(:,:,k)、X _k＝X(:,:,k)、E _k＝E(:,:,k)、U _k＝U(:,:,k)、W _k＝W(:,:,k)和G _k＝G(:,:,k)，其中k＝1,2,...,K，“:”表示取张量对应维的所有元素。根据式(2)，利用交替方向乘子法(AlternatingDirection Method of Multipliers，ADMM)和半二次分裂法对共享SM、共享TC以及核张量进行更新。

接着，根据Tucker分解得到的多被试共享TC矩阵和共享SM矩阵，分别提取感兴趣成分的索引；然后，根据感兴趣成分在共享TC矩阵和共享SM矩阵中的索引，分别提取包含所有被试N个特征的空间特征矩阵和时间特征矩阵；最后，从空间特征矩阵和时间特征矩阵中分别提取每个被试的空间特征矢量和时间特征矢量，利用k-means算法分别根据空间特征和时间特征对被试进行细分类；具体实现步骤如下(见图1)：

第一步：输入数据。以“空间体素×时间×被试”形式输入多被试fMRI数据张量

第二步：根据式(1)和式(2)，对多被试fMRI数据X通过ADMM和半二次分裂法实现Tucker分解，得到核张量

共享SM矩阵

和共享TC矩阵

求解步骤具体如下：

步骤1：参数设置。设置共享成分个数10≤N≤T；设置式(1)中p、δ、λ、γ四个稀疏项参数为(0,1]范围内的值；设置式(4)半二次分裂法中分裂变量的稀疏项参数ξ＞0，以及牛顿法求解分裂变量的最大迭代次数iter_y_max(见步骤4)；设置ADMM最大迭代次数iter_max、最小迭代误差ε_min、最小相对误差Δε_min；

步骤2：初始化。用HOSVD算法对X的分解结果对共享SM矩阵

共享TC矩阵

以及核张量

进行初始化；求解残差张量E＝X-G×₁S×₂B；令G的分裂变量R＝G，S的分裂变量Y＝S(见式(4))，拉格朗日乘子U＝0、W＝0、Q＝0(见式(4))，惩罚参数α₀＝K/||X||_F、β₀＝K/||R||_F；令迭代次数iter＝1，迭代误差ε₀＝1，相对误差Δε₀＝1；

步骤3：更新共享TC矩阵B。由式(2)的增广拉格朗日函数，得到B的更新如式(3)：

式中，I表示单位阵。

步骤4：更新共享SM矩阵S、分裂变量Y、Y的一阶导数Y_d和Y的二阶导数Y_dd。利用半二次分裂法，在式(2)中引入分裂变量Y，则增广拉格朗日函数可以写为：

其中，L₁(B,G,E,R,U,W,α,β)是不包括S的增广拉格朗日项；ξ为稀疏项参数，Q是拉格朗日乘子。由式(4)，推导S的计算公式如下：

再令iter_y＝1，采用牛顿法，每次迭代时利用式(6)-(8)迭代更新分裂变量Y、Y的一阶导数Y_d和Y的二阶导数Y_dd：

Y＝Y-Y_d./Y_dd (6)

Y_d＝ξp·sgn(Y)ο|Y|^p-1+δ(S-Y)-Q (7)

Y_dd＝ξp(p-1)|Y|^p-2-δ1 (8)

直到iter_y＝iter_y_max；其中“./”为矩阵点除运算，p为l_p范数的稀疏参数，“sgn(·)”为符号函数，“ο”为矩阵点乘，“|·|”表示取绝对值，1是与Y相同大小的全1矩阵。

步骤5：更新核张量G和分裂变量R。利用软阈值方法，根据式(9)更新核张量G：

式中

k＝1,2,...,K。对于核张量G的分裂变量R，利用离散李雅普诺夫方程求解：

其中

步骤6：更新残差张量E。利用软阈值方法，根据式(11)更新E：

其中

步骤7：根据式(12)-(14)，更新拉格朗日乘子U、W、Q：

U←U+α(X-R×₁S×₂B-E) (12)

W←W+β(G-R) (13)

Q←Q+δ(Y-S) (14)

步骤8：根据式(15)和(16)，更新惩罚参数α和β：

α←ηα (15)

β←ηβ (16)

式中惩罚参数α和β更新的增长率η大于1。

步骤9：根据式(17)和(18)，更新迭代误差ε_iter和相对误差Δε_iter：

ε_iter＝||X-G×₁S×₂B-E||_F/||X||_F (17)

Δε_iter＝|ε^iter-1-ε^iter|/ε^iter-1 (18)

步骤10：若迭代误差ε_iter小于预设误差阈值ε_min，或者相对误差Δε_iter小于预设误差阈值Δε_min，或者iter大于预设最大迭代次数iter_max，则跳转到步骤11，否则执行iter＝iter+1并跳转到步骤3。

步骤11：输出共享SM矩阵S，共享TC矩阵B以及核张量G。

第三步：分别构建共享TC和共享SM感兴趣成分的参考成分。当构建共享TC感兴趣成分的参考成分时，对于多被试任务态fMRI数据，任务相关参考成分由实验范式与血液动力学反应函数(hemodynamic response function,HRF)线性卷积生成，默认网络(defaultmode network，DMN)参考成分为任务相关参考成分的反相成分；对于任务态fMRI数据其他(非任务相关)成分以及多被试静息态fMRI数据，参考成分采用GICA方法(V.D.Calhoun,T.Adal₁,et al,“A method for making group inferences from functional MRI datausing independent component analysis,”Human Brain Mapping,vol.14,no.3,pp.140-151,2001)获取，具体如下：采用GICA对时间维串联多被试fMRI数据进行分离，得到每被试特有TC，再对所有被试的特有TC取平均，得到共享TC感兴趣成分的参考成分。当构建共享SM感兴趣成分的参考成分时，对于多被试任务态fMRI数据，任务相关参考成分通过组广义线性模型(generalized linear model,GLM)方法构建(L.D.Kuang,Q.H.Lin,etal.,“Shift-invariant canonical polyadic decomposition of complex-valued multi-SubjectfMRI data with a phase sparsity constraint,”IEEE Transactions on MedicalImaging,vol.39,no.4,pp.844-853,2019)；对于任务态fMRI数据其他(非任务相关)成分以及多被试静息态fMRI数据，参考成分选自文献“S.M.Smith,P.T.Fox et al.,“Correspondence of the brain's functional architecture during activation andrest,”Proceedings of the National Academy of Sciences of the United States ofAmerica,vol.106,no.31,pp.13040-13045,2009”的结果。

第四步：分别提取共享TC和共享SM感兴趣成分的索引。基于与感兴趣成分的参考成分相关系数最大的原则，从B＝[b₁,b₂,...,b_N]中提取共享TC的感兴趣成分

得到其索引j∈{1,2,...,N}，其中b_j＝B(:,j)；从S＝[s₁,s₂,...,s_N]中提取共享SM感兴趣成分

得到其索引i∈{1,2,...,N}，其中s_i＝S(:,i)；“:”表示取矩阵对应维的所有元素。

第五步：根据共享TC和共享SM感兴趣成分的索引j和i，提取包含所有被试N个特征的空间特征矩阵

和时间特征矩阵

G_S＝[G(:,j,1),G(:,j,2),...,G(:,j,K)] (19)

G_B＝[G(i,:,1),G(i,:,2),...,G(i,:,K)] (20)

其中“:”表示取张量对应维的所有元素；G_S和G_B的每一列元素分别对应每个被试的N个空间特征矢量和时间特征矢量。

第六步：根据空间特征矩阵G_S和时间特征矩阵G_B，利用k-means算法对K个被试进行细分类。重写式(19)和(20)如下：

G_S＝[g_s1,g_s2,...,g_sK] (21)

G_B＝[g_b1,g_b2,...,g_bK] (22)

其中，

和

分别为被试k对应的空间特征矢量和时间特征矢量，k＝1,2,…,K；利用k-means算法，分别对所有K个被试的空间特征矢量g_s1,g_s2,...,g_sK和时间特征矢量g_b1,g_b2,...,g_bK进行聚类；聚类数目由肘部法则决定，聚类的结果即为被试空间特征细分类结果和时间特征细分类结果。

第七步：输出空间特征矩阵G_S、时间特征矩阵G_B、被试空间特征细分类和时间特征细分类结果。

本发明聚焦多被试fMRI数据的个体空时差异性挖掘问题，提供了一种基于空间稀疏约束Tucker分解核张量的各被试空时特征矢量提取和被试细分类方法，从“空间体素×时间×被试”形式的多被试fMRI数据核张量中有效提取单被试的空间特征矢量和时间特征矢量，体现各被试空间激活和时间过程的形态差异。利用个体空时特征矢量通过k-means算法对被试进行细分类，能为脑认知和脑疾病研究提供新方法和新依据。以10个健康被试的任务态fMRI数据为例，本发明提取的感兴趣DMN成分所对应空间特征矩阵G_S捕捉了每被试特有空间激活形态的不同；利用k-means算法将所有被试分为两组，组1各被试SMs的顶下小叶(inferior parietal lobule，IPL)区域平均激活体素数目比组2多109％(组1：677，组2：324)；感兴趣任务相关成分所对应时间特征矩阵G_B捕捉了每被试特有时间过程形态的不同；利用k-means算法将所有的被试分为两组，组1任务相关时间过程与参考成分的平均相关系数比组2高55.6％(组1：0.84，组2：0.54)。

附图说明

图1为本发明的实现流程图。

图2为DMN成分的多被试空间特征矩阵。

图3为任务相关成分的多被试时间特征矩阵。

图4为IVA-GL所获取每被试特有DMN空间成分及其性能。

图5为DMN空间激活的参考信号。

图6为IVA-GL所获取每被试特有任务相关时间过程成分及其性能。

图7为任务相关时间过程的参考信号。

具体实施方式

下面结合技术方案详细叙述本发明的一个具体实施例。

现有K＝10被试的敲击手指任务态fMRI数据，包括10个健康人。每个被试含有T＝165次全脑扫描，每次全脑扫描图共有153594体素，其中脑内体素V＝59610。选取DMN为共享TC的感兴趣成分、任务相关成分为共享SM的感兴趣成分，进行各被试空时特征矢量提取和k-means聚类，且与IVA-GL输出的各被试DMN空间成分和任务相关时间过程成分进行比较。

第一步：输入数据。以“空间体素×时间×被试”形式输入多被试fMRI数据张量

第二步：根据式(1)和式(2)，对多被试fMRI数据X通过ADMM和半二次分裂法实现Tucker分解，得到核张量

共享SM矩阵

和共享TC矩阵

具体步骤如下：

步骤1：参数设置。设置共享成分个数N＝50；设置式(1)中的p＝0.3、δ＝0.4、λ＝0.4、γ＝0.6四个稀疏项参数；设置式(4)半二次分裂法中分裂变量的稀疏项参数ξ＝0.4，牛顿法求解分裂变量的最大迭代次数iter_y_max＝10；设置ADMM最大迭代次数iter_max＝300、最小迭代误差ε_min＝10^-7、最小相对误差Δε_min＝10^-4；

步骤2：初始化。用HOSVD算法对X的分解结果对共享SM矩阵

共享TC矩阵

以及核张量

进行初始化；求解残差张量E＝X-G×₁S×₂B；令G的分裂变量R＝G，S的分裂变量Y＝S，拉格朗日乘子U＝0、W＝0、Q＝0，惩罚参数α₀＝K/||X||_F、β₀＝K/||R||_F；令迭代次数iter＝1，迭代误差ε₀＝1，相对误差Δε₀＝1；

步骤3：应用式(3)更新共享TC矩阵B；

步骤4：应用式(5)更新共享SM矩阵S；应用式(6)-(8)更新分裂变量Y、Y的一阶导数Y_d和Y的二阶导数Y_dd；

步骤5：应用式(9)更新核张量G；应用式(10)更新分裂变量R；

步骤6：应用式(11)更新残差张量E；

步骤7：根据式(12)-(14)，更新拉格朗日乘子U、W、Q；

步骤8：根据式(15)和(16)，更新惩罚参数α和β；

步骤9：根据式(17)和(18)，更新迭代误差ε_iter和相对误差Δε_iter；

步骤10：若迭代误差ε_iter小于预设误差阈值ε_min，或者相对误差Δε_iter小于预设误差阈值Δε_min，或者iter大于预设最大迭代次数iter_max，则跳转到步骤11，否则执行iter＝iter+1并跳转到步骤3。

步骤11：输出共享SM矩阵S，共享TC矩阵B以及核张量G。

第三步：分别构建共享TC中DMN成分的参考成分和共享SM中任务相关成分的参考成分。DMN参考成分由实验范式与血液动力学反应函数线性卷积生成，并进行反相；任务相关参考成分通过组GLM得到。

第四步：分别提取共享TC中DMN成分的索引和共享SM中任务相关成分的索引。基于与DMN参考成分相关系数最大的原则，从B＝[b₁,b₂,...,b₅₀]中提取DMN成分

得到其索引j＝33，其中b₃₃＝B(:,33)，“:”表示取矩阵对应维的所有元素；基于与任务相关参考成分相关系数最大的原则，从S＝[s₁,s₂,...,s₅₀]中提取任务相关成分

得到其索引i＝5，其中s₅＝S(:,5)。

第五步：根据DMN成分的索引j＝33和任务相关成分的索引i＝5，分别应用式(19)和(20)提取包含每个被试50个空间特征矢量的空间特征矩阵

以及包含每个被试50个时间特征矢量的时间特征矩阵

第六步：根据空间特征矩阵G_S和时间特征矩阵G_B，利用k-means算法对10个被试进行细分类。分别提取空间特征矩阵和时间特征矩阵中10个被试对应的空间特征矢量

和时间特征矢量

k＝1,2,…,10；利用k-means算法，分别对所有10个被试的空间特征矢量g_s1,g_s2,...,g_s10和时间特征矢量g_b1,g_b2,...,g_b10进行聚类；由肘部法则求解空间特征矢量和时间特征矢量的聚类数目，都等于2。

第七步：输出空间特征矩阵G_S(如图2所示)、时间特征矩阵G_B(如图3所示)、被试空间特征细分类结果：组1＝{被试#1,#4,#5,#7}，组2＝{被试#2,#3,#6,#8,#9,#10}、被试时间特征细分类结果：组1＝{被试#1,#2,#3,#5,#6,#7,#8,#9,#10}，组2＝{被试#4}。

图4为IVA-GL所获取每被试特有DMN空间成分及其性能；DMN提取方法具体如下：采用IVA-GL方法以模型阶数N＝50分离10被试fMRI数据，得到每被试的50个空间成分；从“S.M.Smith,P.T.Fox et al.,“Correspondence of the brain's functionalarchitecture during activation and rest,”Proceedings of the National Academyof Sciences of the United States of America,vol.106,no.31,pp.13040-13045,2009”结果中选取DMN空间参考(见图5)，基于与该DMN参考成分相关系数最大的原则，从每被试的50个空间成分中选取DMN成分。

以图4为对照，可以验证本发明基于空间特征矢量进行被试细分类的有效性。由图可见，两组被试(组1＝{被试#1,#4,#5,#7}，组2＝{被试#2,#3,#6,#8,#9,#10})的空间激活形态在IPL区域明显不同。其中，组1被试的IPL激活体素数目均大于等于542，而组2的IPL激活体素数目均小于等于463。组1各被试IPL区域的平均激活体素数比组2多109％(组1：677，组2：324)。然而，若以相关系数为分类指标，则被试#1和#2将属于同一组(与DMN空间参考的相关系数都是0.55)，被试#7和#10也将属于同一组(与DMN空间参考的相关系数都是0.58)，但这些被试在IPL区域存在明显的激活形态差异。所以，本发明能够捕捉每被试DMN成分特有空间激活形态的差异，解决相关系数无法区分各被试空间激活形态差异的问题。

图6为IVA-GL所获取每被试特有任务相关时间过程成分及其性能；任务相关成分提取方法具体如下：采用IVA-GL方法以模型阶数N＝50分离10被试fMRI数据，得到每被试的50个时间过程成分；由实验范式与HRF线性卷积生成的时间过程作为参考成分(见图7)，基于与参考成分相关系数最大的原则，从每被试50个时间过程成分中选取任务相关成分。

以图6为对照，可以验证本发明基于时间特征矢量进行被试细分类的有效性。由图可见，两组被试(组1＝{被试#1,#2,#3,#5,#6,#7,#8,#9,#10}，组2＝{被试#4})的时间过程形态和性能有明显不同。其中，组1各被试的任务相关时间过程噪声小，与参考成分相关系数达到0.70～0.91；组2只有被试#4，相比于组1各被试，被试#4的时间过程所含噪声大，与参考成分的相关系数只有0.54。组1与参考成分的平均相关系数比组2高55.6％(组1：0.84，组2：0.54)。所以，本发明也能够捕捉每被试特有时间过程的形态差异。

综合空间特征矢量和时间特征矢量的被试细分类结果，可以发现，被试4与其他被试在空、时特征上均有所不同。因此，本发明能够为被试细分类提供空间特征差异、时间特征差异、空时特征综合差异等多维度、多角度特征，为脑功能研究和脑疾病诊断提供新的客观依据。

Claims (2)

1.一种提取个体空时特征矢量与被试细分类的方法，是从多被试fMRI数据的空间稀疏约束Tucker分解核张量中提取个体空时特征矢量与被试细分类的方法，其特征在于，

在Tucker分解模型中引入空间稀疏约束，形成如下模型：

其中，

是多被试fMRI数据，V是脑内体素的个数，T是时间点个数，K是被试个数；

是共享SM矩阵，

是共享TC矩阵，

是核张量，

是残差张量，N是共享成分个数；“×₁”和“×₂”为模-1乘积和模-2乘积；“||·||_F”、“||·||₁”、“||·||_p”分别为l_F范数、l₁范数和l_p范数，p为稀疏参数；δ、λ、γ分别是空间稀疏项、核张量稀疏项和残差张量稀疏项参数；式(1)中，空间稀疏约束通过S的l_p范数实现，S和B的低秩约束由l_F范数实现，G和E的稀疏约束通过l₁范数实现；

由式(1)，得到增广拉格朗日函数如下：

式中，

是G的分裂变量，

是拉格朗日乘子，V是脑内体素的个数，T是时间点个数，K是被试个数，N是共享成分个数；α、β是惩罚参数，“<·>”是矩阵内积；

和

分别为张量R、X、E、U、W和G的第k个正面切片，分别满足R _k＝R(:,:,k)、X _k＝X(:,:,k)、E _k＝E(:,:,k)、U _k＝U(:,:,k)、W _k＝W(:,:,k)和G _k＝G(:,:,k)，其中k＝1,2,...,K，“:”表示取张量对应维的所有元素；

根据式(2)，利用ADMM和半二次分裂法对共享SM、共享TC以及核张量进行更新；

根据Tucker分解得到的多被试共享TC矩阵和共享SM矩阵，分别提取感兴趣成分的索引；然后，根据感兴趣成分在共享TC矩阵和共享SM矩阵中的索引，分别提取包含所有被试N个特征的空间特征矩阵和时间特征矩阵；最后，从空间特征矩阵和时间特征矩阵中分别提取每个被试的空间特征矢量和时间特征矢量，利用k-means算法分别根据空间特征和时间特征对被试进行细分类。

2.根据权利要求1所述的一种提取个体空时特征矢量与被试细分类的方法，其特征在于，具体实现步骤如下：

第一步：输入数据；以“空间体素×时间×被试”形式输入多被试fMRI数据张量

第二步：根据式(1)和式(2)，对多被试fMRI数据X通过ADMM和半二次分裂法实现Tucker分解，得到核张量

共享SM矩阵

和共享TC矩阵

第三步：分别构建共享TC和共享SM感兴趣成分的参考成分；当构建共享TC感兴趣成分的参考成分时，对于多被试任务态fMRI数据，任务相关参考成分由实验范式与HRF线性卷积生成，DMN参考成分为任务相关参考成分的反相成分；对于任务态fMRI数据其他成分以及多被试静息态fMRI数据，参考成分采用GICA方法获取，具体如下：采用GICA对时间维串联多被试fMRI数据进行分离，得到每被试特有TC，再对所有被试的特有TC取平均，得到共享TC感兴趣成分的参考成分；当构建共享SM感兴趣成分的参考成分时，对于多被试任务态fMRI数据，任务相关参考成分通过组GLM方法构建；对于任务态fMRI数据其他成分以及多被试静息态fMRI数据，参考成分选自文献“S.M.Smith,P.T.Fox et al.,“Correspondence of thebrain's functional architecture during activation and rest,”Proceedings ofthe National Academy of Sciences of the United States of America,vol.106,no.31,pp.13040-13045,2009”的结果；

第四步：分别提取共享TC和共享SM感兴趣成分的索引；基于与感兴趣成分的参考成分相关系数最大的原则，从B＝[b₁,b₂,...,b_N]中提取共享TC的感兴趣成分

得到其索引j∈{1,2,...,N}，其中b_j＝B(:,j)；从S＝[s₁,s₂,...,s_N]中提取共享SM感兴趣成分

得到其索引i∈{1,2,...,N}，其中s_i＝S(:,i)；“:”表示取矩阵对应维的所有元素；

第五步：根据共享TC和共享SM感兴趣成分的索引j和i，提取包含所有被试N个特征的空间特征矩阵

和时间特征矩阵

G_S＝[G(:,j,1),G(:,j,2),...,G(:,j,K)] (3)

G_B＝[G(i,:,1),G(i,:,2),...,G(i,:,K)] (4)

其中“:”表示取张量对应维的所有元素；G_S和G_B的每一列元素分别对应每个被试的N个空间特征矢量和时间特征矢量；

第六步：根据空间特征矩阵G_S和时间特征矩阵G_B，利用k-means算法对K个被试进行细分类；重写式(3)和式(4)如下：

G_S＝[g_s1,g_s2,...,g_sK] (5)

G_B＝[g_b1,g_b2,...,g_bK] (6)

其中，

和

分别为被试k对应的空间特征矢量和时间特征矢量，k＝1,2,…,K；利用k-means算法，分别对所有K个被试的空间特征矢量g_s1,g_s2,...,g_sK和时间特征矢量g_b1,g_b2,...,g_bK进行聚类；聚类数目由肘部法则决定，聚类的结果即为被试空间特征细分类结果和时间特征细分类结果；

第七步：输出空间特征矩阵G_S、时间特征矩阵G_B、被试空间特征细分类和时间特征细分类结果。

Images