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

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


背景技术：

2.fmri广泛应用于脑功能和神经精神类脑疾病研究。其优点在于高安全性、非侵入式以及毫米级高空间分辨率。多被试fmri数据有5维，包括3维全脑数据、1维全脑扫描次数(即时间点个数)和1维被试个数。在进行盲源分离(blind source separation,bss)时，通常将全脑数据展开成一维体素，此时的多被试fmri数据“空间体素
×
时间
×
被试”也高达3维。
3.张量分解方法因其充分利用数据高维结构信息的优点，非常适合多被试fmri数据的分析。对于体现多被试fmri数据高维空时结构的“空间体素
×
时间
×
被试”形式张量，能够从中提取多被试共享的空间激活成分(spatial maps,sms)及其对应的多被试共享时间过程(timecourses,tcs)。tucker分解是张量分解方法的一种，既能分解出多被试共享sms和共享tcs，还能分解出核张量。核张量中包含了丰富而压缩的原始数据信息，反映了输入高维数据在各个维度上的紧密联系，从中可提取各被试特有的强度信息和特有的空时特征。因此，tucker分解方法是脑功能研究或脑疾病诊断的重要方法论之一，在保持原始高维数据结构的优势下，同时提供多被试群体性与个体性空时脑功能信息。
4.常用tucker分解方法有高阶正交迭代(higher
‑
order orthogonal iteration，hooi)算法和高阶奇异值分解(higher
‑
order svd，hosvd)算法。然而，对于“空间体素
×
时间
×
被试”形式的多被试fmri高维数据张量，这些tucker分解算法均不能提供令人满意的性能。主要原因在于，多被试fmri数据具有高噪性，而且多被试共享sms具有空间稀疏性。这些特性在现有tucker分解模型中均未得到考虑。
5.针对这一问题，本发明采用bahri等人在2019年文章(m.bahri,y.panagakis,and s.zafeiriou,“robust kronecker component analysis,”ieee transactions on pattern analysis and machine intelligence,vol.41,no.10,pp.2365
–
2379,2019)中提出的一种在光学图像降噪中取得优越性能的低秩约束模型robust kronecker component analysis(rkca)用于“空间体素
×
时间
×
被试”fmri张量的降噪；在rkca的低秩约束基础上，继续引入空间稀疏约束，解决rkca算法只适合不具有空间稀疏性的光学图像降噪，而不适于具有空间稀疏性的fmri数据分析问题；提供一种适用于“空间体素
×
时间
×
被试”张量的多被试fmri高维数据tucker分解方法，为脑功能研究或脑疾病诊断同时提供群体性与个体性脑功能信息。


技术实现要素：

6.本发明提供了一种引入空间稀疏约束的多被试fmri数据tucker分解方法，能有效地从“空间体素
×
时间
×
被试”形式的多被试fmri数据中提取共享sm、共享tc和核张量，显
著提升tucker分解在多被试fmri数据分析中的性能。
7.本发明采用的技术方案如下：
8.在rkca的低秩约束tucker分解模型中引入空间稀疏约束，形成如下模型：
[0009][0010]
其中，是多被试fmri数据，v是脑内体素的个数，t是时间点个数，k是被试个数；是共享sm矩阵，是共享tc矩阵，是核张量，是残差张量，n是共享成分个数，
“×
1”和
“×
2”为模
‑
1乘积和模
‑
2乘积；“||
·
||
f”、“||
·
||
p”、“||
·
||
1”分别为l
f
范数、l
p
范数(p为稀疏参数，且0<p≤1)、和l1范数；δ、λ、γ分别是空间稀疏项、核张量稀疏项和残差张量稀疏项参数；式(1)中，空间稀疏约束通过s的l
p
范数实现，s和b的低秩约束由l
f
范数实现，g和e的稀疏约束通过l1范数实现。
[0011]
由式(1)，得到增广拉格朗日函数如下：
[0012][0013]
式中，是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，“:”表示取张量对应维的所有元素。
[0014]
根据式(2)，利用交替方向乘子法(alternating direction method of multipliers，admm)和半二次分裂法对共享sm、共享tc以及核张量进行更新；具体实现步骤如下：
[0015]
第一步：输入数据。以“空间体素
×
时间
×
被试”形式输入多被试fmri数据张量
[0016]
第二步：参数设置。设置共享成分个数10≤n≤t；设置式(1)中p、δ、λ、γ四个稀疏项参数为(0,1]范围内的值；设置式(4)半二次分裂法中分裂变量的稀疏项参数ξ＞0，以及牛顿法求解分裂变量的最大迭代次数iter_y
max
(见第五步)；设置admm最大迭代次数iter
max
、最小迭代误差ε
min
、最小相对误差δε
min
；
[0017]
第三步：初始化。用hosvd算法对x的分解结果对共享sm矩阵共享tc矩阵以及核张量进行初始化；求解残差张量e＝x
‑
g
×1s
×2b；令g的分
裂变量r＝g，s的分裂变量y＝s(见式(4))，拉格朗日乘子u＝0、w＝0、q＝0(见式(4))，惩罚参数α0＝k/||x||
f
、β0＝k/||r||
f
；令迭代次数iter＝1，迭代误差ε0＝1，相对误差δε0＝1；
[0018]
第四步：更新共享tc矩阵b。由式(2)的增广拉格朗日函数，得到b的更新如式(3)：
[0019][0020]
式中，i表示单位阵。
[0021]
第五步：更新共享sm矩阵s、分裂变量y、y的一阶导数y
d
和y的二阶导数y
dd
。利用半二次分裂法，在式(2)中引入分裂变量y，则增广拉格朗日函数可以写为：
[0022][0023]
其中，l1(b,g,e,r,u,w,α,β)是不包括s的增广拉格朗日项；ξ为稀疏项参数，q是拉格朗日乘子。由式(4)，推导s的计算公式如下：
[0024][0025]
再令iter_y＝1，采用牛顿法，每次迭代时利用式(6)
‑
(8)迭代更新分裂变量y、y的一阶导数y
d
和y的二阶导数y
dd
：
[0026]
y＝y
‑
y
d
./y
dd
ꢀꢀ
(6)
[0027][0028]
y
dd
＝ξp(p
‑
1)|y|
p
‑2‑
δ1
ꢀꢀ
(8)
[0029]
直到iter_y＝iter_y
max
；其中“./”为矩阵点除运算，p为l
p
范数的稀疏参数，“sgn(
·
)”为符号函数，为矩阵点乘，“|
·
|”表示取绝对值，1是与y相同大小的全1矩阵。
[0030]
第六步：更新核张量g和分裂变量r。利用软阈值方法，根据式(9)更新核张量g：
[0031][0032]
式中k＝1,2,...,k。对于核张量g的分裂变量r，利用离散李雅普诺夫方程求解：
[0033][0034]
其中
[0035]
第七步：更新残差张量e。利用软阈值方法，根据式(11)更新e：
[0036][0037]
其中
[0038]
第八步：根据式(12)
‑
(14)，更新拉格朗日乘子u、w、q：
[0039]
u
←
u α(x
‑
r
×1s
×2b
‑
e)
ꢀꢀ
(12)
[0040]
w
←
w β(g
‑
r)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0041]
q
←
q δ(y
‑
s)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0042]
第九步：根据式(15)和(16)，更新惩罚参数α和β：
[0043]
α
←
ηα
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0044]
β
←
ηβ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0045]
式中惩罚参数α和β更新的增长率η大于1。
[0046]
第十步：根据式(17)和(18)，更新迭代误差ε
iter
和相对误差δε
iter
：
[0047]
ε
iter
＝||x
‑
g
×1s
×2b
‑
e||
f
/||x||
f
ꢀꢀꢀꢀꢀꢀꢀ
(17)
[0048]
δε
iter
＝|ε
iter
‑1‑
ε
iter
|/ε
iter
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0049]
第十一步：若迭代误差ε
iter
小于预设误差阈值ε
min
，或者相对误差δε
iter
小于预设误差阈值δε
min
，或者iter大于预设最大迭代次数iter
max
，则跳转到第十二步，否则执行iter＝iter 1并跳转到第四步。
[0050]
第十二步：输出共享sm矩阵s，共享tc矩阵b以及核张量g。
[0051]
本发明聚焦多被试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％。
附图说明
[0052]
图1为本发明的实现流程图。
[0053]
图2为本发明与rkca方法所提取两个共享sm和共享tc成分的性能比较。
具体实施方式
[0054]
下面结合技术方案详细叙述本发明的一个具体实施例。
[0055]
现有k＝10个被试的任务态fmri数据，包括10个健康人。每个被试含有t＝165次全脑扫描，每次全脑扫描图共有153594体素，其中脑内体素v＝59610。
[0056]
第一步：输入数据。以“空间体素
×
时间
×
被试”形式输入多被试fmri数据张量
[0057]
第二步：参数设置。设置共享成分个数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；
[0058]
第三步：初始化。用hosvd算法对x的分解结果对共享sm矩阵共享tc
矩阵以及核张量进行初始化；求解残差张量e＝x
‑
g
×1s
×2b；令g的分裂变量r＝g，s的分裂变量y＝s，拉格朗日乘子u＝0、w＝0、q＝0，惩罚参数α0＝k/||x||
f
、β0＝k/||r||
f
；令迭代次数iter＝1，迭代误差ε0＝1，相对误差δε0＝1；
[0059]
第四步：应用式(3)更新共享tc矩阵b；
[0060]
第五步：应用式(5)更新共享sm矩阵s；应用式(6)
‑
(8)更新分裂变量y、y的一阶导数y
d
和y的二阶导数y
dd
；
[0061]
第六步：应用式(9)更新核张量g；应用式(10)更新分裂变量r；
[0062]
第七步：应用式(11)更新残差张量e；
[0063]
第八步：根据式(12)
‑
(14)，更新拉格朗日乘子u、w、q；
[0064]
第九步：根据式(15)和(16)，更新惩罚参数α和β；
[0065]
第十步：根据式(17)和(18)，更新迭代误差ε
iter
和相对误差δε
iter
；
[0066]
第十一步：若迭代误差ε
iter
小于预设误差阈值ε
min
，或者相对误差δε
iter
小于预设误差阈值δε
min
，或者iter大于预设最大迭代次数iter
max
，则跳转到第十二步，否则执行iter＝iter 1并跳转到第四步。
[0067]
第十二步：输出共享sm矩阵s，共享tc矩阵b以及核张量g。
[0068]
第十三步：基于与参考成分相关系数最大的原则，从共享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 during activation and rest,proceedings of the national academy of sciences of the united 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)。