CN113792254B

CN113792254B - 一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法

Info

Publication number: CN113792254B
Application number: CN202111302720.XA
Authority: CN
Inventors: 林秋华; 韩越
Original assignee: Dalian University of Technology
Current assignee: Dalian University of Technology
Priority date: 2021-08-17
Filing date: 2021-11-05
Publication date: 2024-05-28
Anticipated expiration: 2041-11-05
Also published as: CN113792254A

Abstract

一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法，属于生物医学信号处理领域。其在RKCA低秩约束Tucker分解模型中引入了空间稀疏约束，解决了“空间体素×时间×被试”形式多被试fMRI数据的高噪性问题，匹配了多被试共享空间成分的稀疏性，从中有效提取了多被试的共享空间成分、共享时间成分，以及含有丰富被试个体信息的核张量。在10个健康被试任务态fMRI数据分解中，与RKCA方法相比，所提取任务态和默认网络的共享空时成分与参考成分的相关系数分别提升了58.3％和29.6％以及31.7％和31.9％，在空间参考内的激活体素数目分别增加了51.1％和26.0％，为脑认知和脑疾病研究提供了新方法。

Description

一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法

技术领域

本发明属于生物医学信号处理领域，涉及到一种引入空间稀疏约束的多被试功能磁共振成像(functional magnetic resonance imaging，fMRI)数据Tucker分解方法。

背景技术

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

张量分解方法因其充分利用数据高维结构信息的优点，非常适合多被试fMRI数据的分析。对于体现多被试fMRI数据高维空时结构的“空间体素×时间×被试”形式张量，能够从中提取多被试共享的空间激活成分(spatial maps,SMs)及其对应的多被试共享时间过程(timecourses,TCs)。Tucker分解是张量分解方法的一种，既能分解出多被试共享SMs和共享TCs，还能分解出核张量。核张量中包含了丰富而压缩的原始数据信息，反映了输入高维数据在各个维度上的紧密联系，从中可提取各被试特有的强度信息和特有的空时特征。因此，Tucker分解方法是脑功能研究或脑疾病诊断的重要方法论之一，在保持原始高维数据结构的优势下，同时提供多被试群体性与个体性空时脑功能信息。

常用Tucker分解方法有高阶正交迭代(higher-order orthogonal iteration，HOOI)算法和高阶奇异值分解(higher-order SVD，HOSVD)算法。然而，对于“空间体素×时间×被试”形式的多被试fMRI高维数据张量，这些Tucker分解算法均不能提供令人满意的性能。主要原因在于，多被试fMRI数据具有高噪性，而且多被试共享SMs具有空间稀疏性。这些特性在现有Tucker分解模型中均未得到考虑。

针对这一问题，本发明采用Bahri等人在2019年文章(M.Bahri,Y.Panagakis,andS.Zafeiriou,“Robust Kronecker component analysis,”IEEE Transactions onPattern Analysis and Machine Intelligence,vol.41,no.10,pp.2365–2379,2019)中提出的一种在光学图像降噪中取得优越性能的低秩约束模型Robust Kronecker componentanalysis(RKCA)用于“空间体素×时间×被试”fMRI张量的降噪；在RKCA的低秩约束基础上，继续引入空间稀疏约束，解决RKCA算法只适合不具有空间稀疏性的光学图像降噪，而不适于具有空间稀疏性的fMRI数据分析问题；提供一种适用于“空间体素×时间×被试”张量的多被试fMRI高维数据Tucker分解方法，为脑功能研究或脑疾病诊断同时提供群体性与个体性脑功能信息。

发明内容

本发明提供了一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法，能有效地从“空间体素×时间×被试”形式的多被试fMRI数据中提取共享SM、共享TC和核张量，显著提升Tucker分解在多被试fMRI数据分析中的性能。

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

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

其中，是多被试fMRI数据，V是脑内体素的个数，T是时间点个数，K是被试个数；/>是共享SM矩阵，/>是共享TC矩阵，/>是核张量，/>是残差张量，N是共享成分个数，“×₁”和“×₂”为模-1乘积和模-2乘积；“||·||_F”、“||·||_p”、“||·||₁”分别为l_F范数、l_p范数(p为稀疏参数，且0<p≤1)、和l₁范数；δ、λ、γ分别是空间稀疏项、核张量稀疏项和残差张量稀疏项参数；式(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)，利用交替方向乘子法(Alternating Direction Method ofMultipliers，ADMM)和半二次分裂法对共享SM、共享TC以及核张量进行更新；具体实现步骤如下：

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

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

第三步：初始化。用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；

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

式中，I表示单位阵。

第五步：更新共享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_dd＝ξp(p-1)|Y|^p-2-δ1 (8)

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

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

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

其中

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

其中

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

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

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

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

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

α←ηα (15)

β←ηβ (16)

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

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

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

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

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

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

本发明聚焦多被试fMRI高维数据的盲源分离问题，提供了一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法，可以对“空间体素×时间×被试”形式的三维多被试fMRI数据进行共享SM、共享TC以及含有被试丰富个体信息的核张量进行有效估计，为脑认知和脑疾病研究提供方法论。以10个健康被试的任务态fMRI数据为例，本发明与RKCA方法相比，所提取任务态(task)和默认网络(default mode network，DMN)共享SM成分与空间参考的相关系数分别提升了58.3％和29.6％，在空间参考内的激活体素数目分别增加了51.1％和26.0％，共享TC成分与时间参考的相关系数也分别提升了31.7％和31.9％。

附图说明

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

图2为本发明与RKCA方法所提取两个共享SM和共享TC成分的性能比较。

具体实施方式

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

现有K＝10个被试的任务态fMRI数据，包括10个健康人。每个被试含有T＝165次全脑扫描，每次全脑扫描图共有153594体素，其中脑内体素V＝59610。

第二步：参数设置。设置共享成分个数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；

第三步：初始化。用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)更新共享TC矩阵B；

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

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

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

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

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

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

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

第十三步：基于与参考成分相关系数最大的原则，从共享SM矩阵S、共享TC矩阵B中分别提取task成分和DMN成分的共享SM和共享TC。其中Task成分的共享SM参考通过广义线性模型(generalized linear model,GLM)得到，DMN的共享SM参考选自文献“S.M.Smith,P.T.Fox et al.,Correspondence of the brain's functional architecture duringactivation and rest,Proceedings of the National Academy of Sciences of theUnited States of America,vol.106,no.31,pp.13040-13045,2009”的结果。对于共享TC，task成分的参考由实验范式与血液动力学反应函数(hemodynamic response function,HRF)线性卷积生成，DMN成分的参考由task参考成分取反得到。本发明与RKCA方法得到的共享SM以及共享TC的对比结果见图2。与RKCA方法相比，本发明所提取task和DMN共享SM成分与空间参考的相关系数分别提升了58.3％和29.6％(本发明：0.57，0.70；RKCA：0.36，0.54)，共享TC成分的相关系数分别提升了31.7％和31.9％(本发明：0.79，0.62；RKCA：0.60，0.47)，共享SM在参考成分内的激活体素数目分别增加了51.1％和26.0％(本发明：2045，2534；RKCA：1353，2011)。

Claims

1.一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法，其特征在于，在RKCA的低秩约束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以及核张量进行更新；

所述利用ADMM和半二次分裂法对共享SM、共享TC以及核张量进行更新；具体实现步骤如下：

第一步：输入数据；以“空间体素×时间×被试”形式输入多被试fMRI数据张量第二步：参数设置；设置共享成分个数N、式(1)中p、δ、λ、γ四个稀疏项参数、式(4)半二次分裂法中分裂变量的稀疏项参数ξ、牛顿法求解分裂变量的最大迭代次数iter_y_max、ADMM最大迭代次数iter_max、最小迭代误差ε_min和最小相对误差Δε_min；

第三步：初始化；用HOSVD算法对X的分解结果对共享SM矩阵S、共享TC矩阵B以及核张量G进行初始化；求解残差张量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；

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

式中，I表示单位阵；

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

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

Y＝Y-Y_d./Y_dd(6)

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

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

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

式中对于核张量G的分裂变量R，利用离散李雅普诺夫方程求解：/>其中第七步：更新残差张量E；利用软阈值方法，根据式(11)更新E：/>其中/>第八步：根据式(12)-(14)，更新拉格朗日乘子U、W、Q：

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

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

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

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

α←ηα (15)

β←ηβ (16)

式中η为惩罚参数α和β更新的增长率；

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

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

第十一步：若迭代误差ε_iter小于预设误差阈值ε_min，或者相对误差Δε_iter小于预设误差阈值Δε_min，或者iter大于预设最大迭代次数iter_max，则跳转到第十二步，否则执行iter＝iter+1并跳转到第四步；

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

2.根据权利要求1所述的一种引入空间稀疏约束的多被试fMRI数据Tucker分解方法，其特征在于：共享成分个数10≤N≤T；式(1)中p、δ、λ、γ四个稀疏项参数为(0,1]范围内的值；式(4)半二次分裂法中分裂变量的稀疏项参数ξ＞0；惩罚参数α和β更新的增长率η>1。