CN105069307A - 一种结合ICA与移不变CPD的多被试fMRI数据分析方法 - Google Patents

Abstract

一种结合ICA与移不变CPD的多被试fMRI数据分析方法，属于医学信号处理领域。采用以CPD为中心，以ICA为预处理环节的主辅结合方式；ICA预处理为移不变CPD提供联合混合矩阵；移不变CPD采用秩一估计方法，从联合混合矩阵中逐一分解出多被试共享TC、与共享TC相对应的各被试时延以及各被试强度；利用移不变CPD的输出重构联合混合矩阵，采用最小二乘法估计多被试共享SM。本发明能够对多被试任务态fMRI数据进行稳定而更为有效的分析，而且当被试间SM差异性和TC差异性较大时，优势更为明显，所估计的共享SM成分和共享TC成分与先验参考信号有着更高的相关性，TC大时延的估计精度高，且计算量小。

Description

一种结合ICA与移不变CPD的多被试fMRI数据分析方法

技术领域

本发明涉及医学信号处理领域，特别是涉及一种多被试功能磁共振成像(functionalmagneticresonanceimaging,fMRI)数据的分析方法。

背景技术

fMRI数据是由磁共振成像扫描仪对被试(健康人或病人)大脑进行扫描所采集到的脑功能数据，具有无损伤和空间分辨率高等优势。通过利用多被试fMRI数据分析方法，人们能够提取各被试的共享脑空间激活区(spatialmap,SM)、共享时间过程(timecourse,TC)，以及各被试的强度信息。这些信息对于脑功能研究和临床诊断具有重要价值。

多被试fMRI数据是有着空间、时间和被试三个维度的高维数据。最为有效的分析方法应能同时利用fMRI数据的统计信息和结构信息。因此，张量分解方法近年来得到了广泛关注。然而，fMRI数据的信噪比很低，而且即使在相同的任务刺激下，多被试间也存在着SM差异性和TC差异性，所以现有的张量分解方法尚存在鲁棒性差、分解成分性能低等问题。为此，Beckmann和Smith于2005年在文章“Beckmann,C.F.,Smith,S.M.TensorialextensionsofindependentcomponentanalysisformultisubjectfMRIanalysis.NeuroImage25,294-311”中提出了一种ICA与CPD(canonicalpolyadicdecomposition，一种典型的张量分解算法)相结合的多被试fMRI数据分析方法。该方法首先利用ICA获取多被试的共享SM成分，再利用CPD获取多被试的共享TC成分和各被试的强度信息。该方法的优点是，利用ICA为CPD加入了SM空间独立性约束，解决了因fMRI低信噪比及被试间SM差异性所导致的算法鲁棒性问题，并在一定程度上提高了共享SM成分的性能。

然而，Beckmann和Smith方法仍然存在两个问题：第一，没有考虑多被试间的TC差异性。因为不同的被试在响应时间或血液动力学时延方面存在差异，所以各被试的共享TC会存在时延差异。像SM差异性一样，TC差异性也影响张量分解方法的鲁棒性和分解成分性能。等人2008年在文章“M.,Hansen,L.K.,Arnfred,S.M.,Lim,L.H.,Madsen,K.H.Shift-invariantmultilineardecompositionofneuroimagingdata.NeuroImage42,1439-1450”中提出了一种移不变CPD方法，将TC的时延纳入CPD。然而，因为没有考虑SM空间差异性，直接应用移不变CPD方法分析多被试fMRI数据时，其性能甚至劣于Beckmann和Smith方法。

第二，Beckmann和Smith方法没有一个整体的目标函数，在ICA与CPD的结合过程中，ICA和CPD的两个目标函数完全独立并行，不分主次。其结果是，ICA和CPD之间的迭代运算易于发散，进而导致分解失败。因此，Beckmann和Smith方法对多被试fMRI数据的分析结果尚不能令人满意。

发明内容

本发明的目的在于，提供一种结合ICA与移不变CPD的多被试fMRI数据分析方法，同时考虑被试间的SM差异性和TC差异性，并将ICA用作移不变CPD的预处理环节，解决Beckmann和Smith方法未考虑TC差异性以及缺乏整体目标函数的问题，显著改善张量分解方法的鲁棒性和分解成分性能。

本发明的技术方案是，在ICA与CPD结合过程中，采用以CPD为中心、以ICA为预处理环节的主辅结合方式；采用移不变CPD嵌入各被试的TC时延，实现对被试间SM差异性问题和TC差异性问题的同时处理。ICA预处理为移不变CPD提供联合混合矩阵；移不变CPD采用秩一估计方法，从联合混合矩阵中逐一分解出多被试共享TC、与共享TC相对应的各被试时延，以及各被试强度；最后，利用移不变CPD的输出重构联合混合矩阵，采用最小二乘法估计多被试共享SM。具体实现步骤如下：

第一步：输入多被试fMRI数据k＝1,…,K。其中K表示被试数目；J表示时间维的全脑扫描次数；I表示空间维的脑内体素数目。

第二步：进行两级PCA(principlecomponentanalysis)压缩。采用文章“Erhardt,E.B.,Rachakonda,S.,Bedrick,E.J.,Allen,E.A.,Adali,T.,Calhoun,V.D.,2011.Comparisonofmulti-subjectICAmethodsforanalysisoffMRIdata.HumanBrainMapping32,2075-2095”中的两级PCA方法对多被试fMRI数据进行压缩。设多被试共享SM和共享TC的成分数是N。将原始的多被试时间维串联数据压缩为Z＝VX，其中，第一级PCA对每个被试数据进行PCA压缩，使时间维由J降为N₁，然后再做第二级PCA压缩，将多被试时间维串联数据由N₁K×I压缩到N×I。

第三步：ICA预处理。对Z进行ICA，得到分离矩阵进而计算得到联合混合矩阵上标表示伪逆。

第四步：移不变CPD分解。M₀的每一列包含一个共享TC成分以及各被试的时延和强度信息，所以本发明采用移不变秩一估计方法逐一分解M₀的每一列。设M₀的第n列(n＝1,…,N)包含第r个(r＝1,…,N)成分：

$m_n=\left[\begin{array}{c}{b}_{r\left(1\right)}{c}_{1r}+{\mathrm{\ϵ}}_{r\left(1\right)}\\ b_{r\left(2\right)}{c}_{2r}+{\mathrm{\ϵ}}_{r\left(2\right)}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ b_{r\left(K\right)}{c}_{Kr}+{\mathrm{\ϵ}}_{r\left(K\right)}\end{array}\right]---\left(1\right)$

其中，b_r(k)包含多被试的第r个共享TC成分，以及第k个被试的时延信息τ_kr。当τ_kn＞0时，b_r(k)由b_r左移τ_kn点获取；当τ_kn＜0时，b_r(k)由b_r右移τ_kn点获取。另外，c_kr表示第k个被试的强度信息，ε_r(k)表示其他成分的串扰。

定义第k列矩阵化式(1)如下：

对做离散傅里叶变换，得到其频域形式且

${\tilde{m}}_{nk}={\tilde{b}}_r{c}_{kr}\exp \[-\mathrm{\ι}2\mathrm{\π}\frac{f-1}{J}{\mathrm{\τ}}_{kr}\]+{\widetilde{\mathrm{\ϵ}}}_{r\left(k\right)}---\left(3\right)$

式中exp[·]为指数函数，令各被试b_r成分的时延为τ_r＝[τ_1r,…,τ_Kr]^T、各被试的强度为c_r＝[c_1r,…,c_Kr]^T，在随机初始化b_r、τ_r和c_r之后，对其进行迭代估计，直至收敛或达到最大迭代次数。首先，在频域上采用最小二乘法估计b_r，并在时域上进行归一化处理：

然后，采用等人于2007年在文章“M.,Madsen,K.H.,Hansen,L.K.Shiftednon-negativematrixfactorization.ProceedingsoftheIEEESignalProcessingSocietyWorkshop,139-144”和文章“M.,Madsen,K.H.,Hansen,L.K.Shiftedindependentcomponentanalysis.InternationalConferenceonIndependentComponentAnalysisandSourceSeparation,89-96”中提出的算法估计时延τ_r。如果是整数时延，最大化如果是非整数时延，最大化最后，在时域上采用最小二乘法估计c_r并进行归一化：

第五步：联合混合矩阵重构。利用上述移不变CPD的分解结果b_r、τ_r和c_r，重构联合混合矩阵如下：

第六步：采用最小二乘法估计共享SM成分：

包含N个共享SM成分。

第七步：输出多被试共享SM成分s₁,…,s_N、共享TC成分b₁,…,b_N、与b₁,…,b_N相对应的各被试时延τ₁,…,τ_N，以及各被试的强度信息c₁,…,c_N。

本发明所达到的效果和益处是，与Beckmann和Smith方法相比，本发明能够对多被试任务态fMRI数据进行稳定而更为有效的分析，而且当被试间的TC时延差异性及SM空间差异性较大时，优势更为明显。例如，DMN(defaultmodenetwork)是一种重要的非任务相关成分，广泛用于脑功能和脑疾病研究。与任务相关成分相比，不同被试DMN成分的SM差异性和TC差异性更大。所以，在针对16被试敲击手指的fMRI数据分析中，Beckmann和Smith方法对DMN成分的分解性能不高。相比之下，本发明分解的共享SM和共享TC与先验参考信号有着更高的相关性，且能以高精度估计出TC的大时延。此外，本发明计算量小，在数据量巨大的多被试fMRI数据分析中具有非常突出的性能优势。因此，本发明在脑功能研究和脑疾病诊断方面具有良好的应用前景。

附图说明

附图是本发明分析多被试fMRI数据的工作流程图。

具体实施方式

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

现有16被试在执行敲打手指任务下采集的fMRI数据，即K＝16。每个被试都进行了J＝165次扫描，每次扫描都获得了53×63×46的全脑数据，脑内体素数I＝59610。假设16被试共享SM和共享TC的成分数N＝35，采用本发明进行fMRI数据分析的步骤如附图所示。

第一步：输入多被试fMRI数据k＝1,…,16，

第二步：进行两级PCA压缩。第一级PCA将每个被试数据由165×59610压缩为50×59610，即N₁＝50；第二级PCA将16被试时间维串联数据由800×59610压缩到35×59610；获取压缩阵

第三步：ICA预处理。采用fastICA算法对35×59610的多被试fMRI数据进行ICA，非线性函数选择为tanh，算法停止准则为|ΔW|＜10^-6或者达到最大迭代次数1000。得到分离矩阵进而计算得到联合混合矩阵

第四步：移不变CPD分解。采用移不变秩一估计方法逐一分解M₀的每一列，得到16被试的35个共享TC成分b₁,…,b₃₅、与b₁,…,b₃₅相对应的各被试时延τ₁,…,τ₃₅，以及各被试的强度信息c₁,…,c₃₅。

第五步：联合混合矩阵重构。利用上述移不变CPD的分解结果b₁,…,b₃₅、τ₁,…,τ₃₅和c₁,…,c₃₅，按照公式(6)重构联合混合矩阵

第六步：按照公式(7)，采用最小二乘法估计35个共享SM成分s₁,…,s₃₅。

第七步：输出多被试共享SM成分s₁,…,s₃₅、共享TC成分b₁,…,b₃₅、与b₁,…,b₃₅相对应的各被试时延τ₁,…,τ₃₅，以及各被试的强度信息c₁,…,c₃₅。

Claims (2)

1.一种结合ICA与移不变CPD的多被试fMRI数据分析方法，其特征是，在ICA与CPD结合过程中，采用以CPD为中心、以ICA为预处理环节的主辅结合方式；采用移不变CPD嵌入各被试的TC时延，实现对被试间SM差异性问题和TC差异性问题的同时处理；ICA预处理为移不变CPD提供联合混合矩阵；移不变CPD采用秩一估计方法，从联合混合矩阵中逐一分解出多被试共享TC、与共享TC相对应的各被试时延，以及各被试强度；最后，利用移不变CPD的输出重构联合混合矩阵，采用最小二乘法估计多被试共享SM。

2.根据权利要求1所述的一种结合ICA与移不变CPD的多被试fMRI数据分析方法，其特征包括以下步骤：

第一步：输入多被试fMRI数据k＝1,...,K；其中K表示被试数目；J表示时间维的全脑扫描次数；I表示空间维的脑内体素数目；

第二步：进行两级PCA压缩；采用两级PCA方法对多被试fMRI数据进行压缩；设多被试共享SM和共享TC的成分数是N，将原始的多被试时间维串联数据压缩为Z＝VX，其中，第一级PCA将每个被试数据的时间维由J降为N₁，第二级PCA将多被试时间维串联数据由N₁K×I压缩到N×I；

第三步：ICA预处理；对Z进行ICA，得到分离矩阵进而计算得到联合混合矩阵上标表示伪逆；

第四步：移不变CPD分解；采用移不变秩一估计方法逐一分解M₀的每一列，得到共享TC成分b₁,...,b_N、与b₁,...,b_N相对应的各被试时延τ₁,...,τ_N，以及各被试的强度信息c₁,...,c_N；其中，在频域上采用最小二乘法估计b_r，并在时域上进行归一化处理：

b_r＝b_r/||b_r||

然后，采用的算法估计时延τ_r；最后，在时域上采用最小二乘法估计c_r并进行归一化：

c_r＝c_r/||c_r||

第五步：联合混合矩阵重构；利用上述移不变CPD的分解结果重构联合混合矩阵如下：

第六步：采用最小二乘法估计N个共享SM成分

第七步：输出多被试共享SM成分s₁,...,s_N、共享TC成分b₁,...,b_N、与b₁,...,b_N相对应的各被试时延τ₁,...,τ_N，以及各被试的强度信息c₁,...,c_N。