CN109700462B - 多被试复数fMRI数据移不变CPD分析方法 - Google Patents

Abstract

一种引入空间源相位稀疏约束的多被试复数fMRI数据移不变CPD分析方法，属于医学信号处理领域。在移不变CPD算法的基础上，利用交替最小二乘法更新被试共享SM成分、共享TC成分和各被试强度；将基于实数数据的被试时延估计扩展到复数数据的时延估计；对共享SM成分采用空间源相位稀疏约束方式再次更新，具体包括两步：首先对共享SM成分进行相位校正；然后对校正后的空间源相位的大相位体素，利用平滑L0范数逼近函数，对共享SM成分更新。本发明能有效提取出多被试复数fMRI数据中更为全面的共享复数脑功能信息，这些信息在今后脑认知和脑疾病研究都有非常好的应用前景。

Description

多被试复数fMRI数据移不变CPD分析方法

技术领域

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

背景技术

fMRI是脑科学研究的重要技术之一，具有空间分辨高且安全无侵入等优势。由磁共振扫描仪获取的fMRI数据本质是复数的，包括幅值和相位。虽然目前大部分研究集中在幅值fMRI数据分析，但越来越多文献表明，只有利用复数fMRI数据，才能提取完整的脑功能信息。尽管复数fMRI数据因为引入相位数据而具有高噪声性，但已经有了有效的消噪方案，即发明专利“林秋华,于谋川,龚晓峰,丛丰裕.一种对复数fMRI数据进行ICA分析的后处理消噪方法.中国,CN201410191416.6”和文献“M.C.Yu,Q.H.Lin,L.D.Kuang,X.F.Gong,F.Cong,and V.D.Calhoun.ICA of full complex-valued fMRI data using phaseinformation of spatial maps.Journal of Neuroscience Methods,vol.249,pp.75-91,2015”提出的空间源相位消噪方法。该方法表明，在复数SM成分中，感兴趣体素的空间源相位集中在小相位范围[-π/4,π/4]，而噪声体素分散在大相位范围[-π,-π/4)∪(π/4,π]，称之为空间源相位的小相位特性。

多个被试的fMRI数据分析能够获取群体性特征，在疾病研究方面比单被试分析更有意义。因此，研究多被试复数fMRI数据的分析方法变得越来越重要。多被试fMRI数据易于表示为一个三维张量(空间维×时间维×被试维)。因此，张量分解算法适于分析多个被试的fMRI数据，但需要解决多被试fMRI数据与张量模型的失配问题。

首先，由于各被试脑响应的快慢和先后不可避免会存在差异，多被试fMRI数据存在时间差异性。针对该问题，可采用

等人在2008年文章“M.

L.K.Hansen,S.M.Arnfred,L.H.Lim,and K.H.Madsen.Shift-invariant multilinear decompositionof neuroimaging data.Neuroimage,vol.42,pp.1439–1450,2008.”中提出的移不变CPD(canonical polyadic decomposition)算法，该算法物理意义明确，将多被试幅值fMRI数据分解为各被试间共享脑空间激活成分(spatial map,SM)、共享时间过程成分(timecourse,TC)、各被试时延信息、以及各被试的强度差异信息。这些信息能够为脑功能研究和脑疾病诊断提供重要特征。

其次，多被试fMRI数据还具有空间差异性，主要体现为各被试的脑空间激活大小、位置、角度差异。目前，对于多被试幅值fMRI数据，空间差异性的解决方案有施加独立性约束，而多被试复数fMRI数据的空间差异性问题的解决方案还未见报道。

发明内容

本发明的目的在于，提供一种引入空间源相位稀疏约束的多被试复数fMRI数据移不变CPD分析方法，通过增加对共享SM成分的空间源相位稀疏约束，解决多被试复数fMRI数据的空间差异性问题。与此同时，因为空间源相位具有消噪特性，该方法能起到一定的消噪作用，进一步提升算法的整体性能。

本发明的技术方案是，在移不变CPD算法的基础上，利用交替最小二乘法更新被试共享SM成分、共享TC成分和各被试强度；将基于实数数据的被试时延估计扩展到复数数据的时延估计；对共享SM成分采用空间源相位稀疏约束方式再次更新，具体包括两步：首先对共享SM成分进行相位校正；然后对校正后的空间源相位的大相位体素，利用平滑L0范数逼近函数，对共享SM成分更新。具体实现步骤如下：

第一步：输入多被试复数fMRI数据

其中V、J和K分别表示脑内体素数目(空间维)、扫描次数(时间维)和被试数目(被试维)。

第二步：初始化。设成分个数为N(N为大于0的正整数)。随机初始化共享SM成分

共享TC成分

和被试强度

初始化被试时延

为零矩阵。令迭代次数iter＝0，相对误差Δε_iter＝1，计算迭代误差ε_iter：

其中，式(1)也是移不变CPD算法模型，τ_k,n表示第k被试第n个成分的时延，在此令其为整数。b_n(j-τ_k,n)表示为b_j,n时移了τ_k,n个点，具体地，若τ_k,n＞0，第k被试第n个TC成分

相对共享TC成分b_n循环左移τ_k,n个点，否则若τ_k,n＜0，则循环右移|τ_k,n|个点。

第三步：更新共享TC成分B。采用

等人提出的移不变CPD算法中共享TC成分B的更新方法(交替最小二乘法)对B进行更新。

第四步：对共享SM成分S和X进行降维。采用发明专利“邝利丹,林秋华,龚晓峰,丛丰裕.一种适于多被试fMRI数据分析的快速移不变CPD方法.中国,CN201811510882.0”，对S和X降维成

和

从而保证第五步的被试时延快速估计。其中

为X的1模展开形式，将

张量化为

第五步：更新被试时延

由于多被试复数fMRI数据是复数的，不能直接采用

等人提出的基于实数数据的被试时延估计方法。下面对复数数据的被试时延进行估计。对于第k被试第n'(n'＝1,…,N)个成分的时延τ_k,n'估计，首先进行如下一些定义。令

为

进行3模展开矩阵

的第k行向量，定义向量

满足

为第n'个成分的信号，其中c_k,n为C的元素；

为

的第n列向量。

的元素表示为

(j＝1,…,J；n＝1,…,N)。定义Re{·}和Im{·}分别表示取实部和虚部，并做如下运算：

其中，j＝1,…,J，

为

的第n'个元素。接着将g_RR(j)、g_IR(j)、g_RI(j)和g_II(j)(j＝1,…,J)转换到频域形式

和

f＝1,…,F。继续做如下运算：

其中，上标^*表示取共轭，

和

分别表示对共享TC成分b_n'的实部和虚部转换到频域形式后的第f个元素。再将

(f＝1,…,F)傅里叶反变换到时域形式φ_RRR(j)、φ_RRI(j)、φ_IRR(j)、φ_IRI(j)、φ_RIR(j)、φ_RII(j)、φ_IIR(j)和φ_III(j)，j＝1,…,J。

令

(j＝1,…,J)满足

的绝对值最大化实际上等价于均方误差最小化，那么第k个被试第n'个成分的时延τ_k,n'可通过最大化式(5)来获得，具体实现是：

在此，被试时延τ_k,n'为整数。对被试时延

的所有元素τ_k,n'(k＝1,…,K，n'＝1,…,N)按照上述方法进行更新。

第六步：更新共享SM成分S。采用

等人提出的移不变CPD算法中更新S的方法(交替最小二乘法)对S更新。

第七步：对共享SM成分S进行相位校正。首先求取联合混合矩阵

其元素满足z_j+k(J-1),n＝c_k,nb_n(j-τ_k,n)。对每个成分n(n＝1,…,N)的共享SM成分s_n，根据发明专利“林秋华,于谋川,龚晓峰,丛丰裕.一种对复数fMRI数据的ICA估计成分进行相位校正的方法.中国,CN201410189199.7”，将专利中a_n替换成联合混合向量z_n，即通过对exp(-iθ)z_n(0≤θ≤π，exp{·}为指数函数)进行实部能量最大化，获取相位校正旋转角度θ_n，得到相位校正后的共享SM成分

第八步：对共享SM成分S进行空间源相位稀疏约束更新。将S所有体素相位值从大到小排序，记

为第V/3个体素的相位值。对S的V/3大相位体素值，采用L0范数平滑函数(详见文章“H.Mohimani,M.Babaie-Zadeh,and C.Jutten.A fast approach forovercomplete sparse decomposition based on smoothed l0 norm.IEEE Transactionson Signal Processing,vol.57,no.1,pp.289-301,2009”)进行稀疏约束：

其中，s_v,n为S的元素，v＝1,…,V,n＝1,...,N，F_σ(|s_v,n|)＝f_σ(|s_v,n|)，且

其中参数σ充分小时，f_σ(|s_v,n|)接近L0范数；σ越大，f_σ(|s_v,n|)越平滑。具体采用最速下降法对S进行更新：

其中λ为正定步长，ΔS的每个元素Δs_v,n(v＝1,…,V,n＝1,…,N)满足

Δs_v,n＝exp{θ(s_v,n)}f′_σ(|s_v,n|) (10)

在此，θ(s_v,n)为s_v,n(v＝1,…,V；n＝1,…,N)的相位值，f′_σ(|s_v,n|)为f_σ(|s_v,n|)的一阶导数，满足

第九步：更新被试强度C。根据

等人提出的移不变CPD算法中更新C的方法(交替最小二乘法)对C更新。

第十步：计算误差。令iter←iter+1；σ←0.999σ，使σ缓慢变化；根据式(1)，计算本次迭代误差ε_iter，以及相对误差Δε_iter：

Δε_iter＝|(ε_iter-1-ε_iter)/ε_iter-1| (12)

第十一步：若ε_iter小于预设误差阈值ε_{iter_min}，跳转到第十四步，否则执行第十二步。

第十二步：若Δε_iter小于预设相对误差阈值Δε_{iter_min}，跳转到第十四步，否则执行第十三步。

第十三步：若iter大于预设最大迭代次数iter_max，跳转到第十四步，否则执行第三步。

第十四步：对共享SM成分S进行相位消噪。对每个成分n(n＝1,…,N)的共享SM成分s_n，采用发明专利“林秋华,于谋川,龚晓峰,丛丰裕.一种对复数fMRI数据进行ICA分析的后处理消噪方法.中国,CN201410191416.6”的方法进行相位消噪。

第十五步：输出相位消噪的共享SM成分S、共享TC成分B、被试时延

与被试强度C。

本发明所达到的效果和益处是，能够对任务态多被试复数fMRI数据的任务相关成分进行有效估计。在10被试敲击手指任务的复数fMRI数据分析实验中，与移不变CPD算法(被试时延估计采用本发明的时延估计方式)和tensor ICA算法相比，本发明估计的共享任务相关SM成分和共享任务相关TC成分的性能分别提升21.7％～36.6％和6.25％～16.44％。另外，与幅值多被试fMRI分析(ICA和移不变CPD相结合的方法，详见发明专利“林秋华，邝利丹，龚晓峰，丛丰裕.一种结合独立成分分析与移不变规范多元分解的多被试功能核磁共振成像数据分析方法，CN201510510622.3”)相比，本发明提取的共享任务相关SM成分在任务相关区域(初级运动区和辅助运动区)中连续激活体素数目多了178.7％。因此，本发明能有效提取出多被试复数fMRI数据中更为全面的共享复数脑功能信息，这些信息在今后脑认知和脑疾病研究都有非常好的应用前景。

附图说明

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

具体实施方式

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

现有10被试执行敲击手指任务下采集的复数fMRI数据，即K＝10。每个被试都进行了J＝165次扫描，每次扫描采集了53×63×46的脑图像，去掉脑外数据体素，保留脑内数据体素V＝59610。假设共享SM和TC成分的成分个数N＝50，采用本发明进行多被试复数fMRI数据分析的步骤如附图所示。

第一步：输入多被试复数fMRI数据

第二步：初始化。随机初始化共享TC成分

共享SM成分

被试强度

初始化被试时延

为零矩阵，步长λ＝4，σ＝2，iter＝0，相对误差Δε_iter＝1，根据式(1)计算迭代误差ε_iter。

第三步：更新共享TC成分B。采用

等人提出的移不变CPD算法中共享TC成分B的更新方法(交替最小二乘法)对B更新，其中令B转换到频域形式的频点数F＝J。

第四步：对共享SM成分S和X进行降维。采用发明专利“邝利丹,林秋华,龚晓峰,丛丰裕.一种适于多被试fMRI数据分析的快速移不变CPD方法.中国,CN201811510882.0”的方法，将S和X降维成

和

保证第五步的被试时延快速估计。

第五步：更新被试时延

对于每个被试每个成分的时延τ_k,n'，首先，根据公式(2)，计算向量

然后，根据公式(3)，求解g_RR(j)、g_IR(j)、g_RI(j)和g_II(j)，j＝1,…,165，并转换到频域形式

f＝1,…,165。根据公式(4)，求解

f＝1,…,165，再傅里叶反变换到时域形式φ_RRR(j)、φ_RRI(j)、φ_IRR(j)、φ_IRI(j)、φ_RIR(j)、φ_RII(j)、φ_IIR(j)和φ_III(j)，j＝1,…,165。最后，根据式(5)和(6)求解τ_k,n'，k＝1,…,10，n'＝1,…,50。

第六步：更新共享SM成分S。采用

等人提出的移不变CPD算法中更新S的方法对S更新。

第七步：对共享SM成分S进行相位校正。对每个成分n(n＝1,…,N)的共享SM成分s_n，采用发明专利“林秋华,于谋川,龚晓峰,丛丰裕.一种对复数fMRI数据的ICA估计成分进行相位校正的方法.中国,CN201410189199.7”进行相位校正，并将该发明专利中的TC成分b_n替换成联合混合向量z_n，得到相位校正后的共享SM成分

第八步：对共享SM成分S进行空间源相位稀疏约束更新。根据式(9)～(11)，对共享SM成分S进行更新。

第九步：更新被试强度C。采用

等人提出的移不变CPD算法中更新C的方法对C更新。

第十步：计算误差。令iter←iter+1，σ←0.999σ，根据式(1)和式(12)，分别计算本次迭代误差ε_iter和相对误差Δε_iter。

第十一步：预设误差阈值ε_{iter_min}＝10^-4。若ε_iter＜ε_{iter_min}，跳转到第十四步，否则执行第十二步。

第十二步：预设相对误差阈值Δε_{iter_min}＝10^-6。若Δε_iter＜Δε_{iter_min}，跳转到第十四步，否则执行第十三步。

第十三步：预设最大迭代次数iter_max＝500。若iter＞iter_max，跳转到第十四步，否则执行第三步。

第十四步：对共享SM成分S进行相位消噪。对每个成分n(n＝1,…,N)的共享SM成分s_n，采用发明专利“林秋华,于谋川,龚晓峰,丛丰裕.一种对复数fMRI数据进行ICA分析的后处理消噪方法.中国,CN201410191416.6”的方法进行相位消噪。

第十五步：输出消噪后的共享SM成分S、共享TC成分B、被试时延

与被试强度C。

Claims (1)

1.一种多被试复数fMRI数据移不变CPD分析方法，其特征包括以下步骤：

第一步：输入多被试复数fMRI数据

其中V、J和K分别表示脑内体素数目、扫描次数和被试数目；

第二步：初始化；设成分个数为N，N为大于0的正整数；随机初始化共享SM成分

共享TC成分

和被试强度

初始化被试时延

为零矩阵；令迭代次数iter＝0，相对误差Δε_iter＝1，计算迭代误差ε_iter：

其中，式(1)也是移不变CPD算法模型，τ_k,n表示第k被试第n个成分的时延，令其为整数；b_n(j-τ_k,n)表示为b_j,n时移了τ_k,n个点，具体地，若τ_k,n＞0，第k被试第n个TC成分

相对共享TC成分b_n循环左移τ_k,n个点，否则若τ_k,n＜0，则循环右移|τ_k,n|个点；

第三步：更新共享TC成分B；采用移不变CPD算法中共享TC成分B的更新方法对B进行更新；

第四步：对共享SM成分S和X进行降维；对S和X降维成

和

从而保证第五步的被试时延快速估计；其中

为X的1模展开形式，将

张量化为

第五步：更新被试时延

对复数数据的被试时延进行估计；对于第k被试第n'(n'＝1,…,N)个成分的时延τ_k,n'估计，首先进行如下一些定义；令

为

进行3模展开矩阵

的第k行向量，定义向量

满足

为第n'个成分的信号，其中c_k,n为C的元素；

为

的第n列向量；

的元素表示为

定义Re{·}和Im{·}分别表示取实部和虚部，并做如下运算：

其中，j＝1,…,J，

为

的第n'个元素；接着将g_RR(j)、g_IR(j)、g_RI(j)和g_II(j)(j＝1,…,J)转换到频域形式

和

f＝1,…,F；继续做如下运算：

其中，上标^*表示取共轭，

和

分别表示对共享TC成分b_n'的实部和虚部转换到频域形式后的第f个元素；再将

和

傅里叶反变换到时域形式φ_RRR(j)、φ_RRI(j)、φ_IRR(j)、φ_IRI(j)、φ_RIR(j)、φ_RII(j)、φ_IIR(j)和φ_III(j)，j＝1,…,J；令

满足

的绝对值最大化实际上等价于均方误差最小化，那么第k个被试第n'个成分的时延τ_k,n'通过最大化式(5)来获得，具体实现是：

被试时延τ_k,n'为整数；对被试时延

的所有元素τ_k,n'，k＝1,…,K，n'＝1,…,N；按照上述方法进行更新；

第六步：更新共享SM成分S；采用移不变CPD算法中更新S的方法对S更新；

第七步：对共享SM成分S进行相位校正；首先求取联合混合矩阵

其元素满足z_j+k(J-1),n＝c_k,nb_n(j-τ_k,n)；对每个成分n的共享SM成分s_n，n＝1,…,N，将TC成分b_n替换成联合混合向量z_n，即通过对exp(-iθ)z_n进行实部能量最大化，获取相位校正旋转角度θ_n，得到相位校正后的共享SM成分

0≤θ≤π，exp{·}为指数函数；

第八步：对共享SM成分S进行空间源相位稀疏约束更新；将S所有体素相位值从大到小排序，记

为第V/3个体素的相位值；对S的V/3大相位体素值，采用L0范数平滑函数进行稀疏约束：

其中，s_v,n为S的元素，v＝1,…,V,n＝1,...,N，F_σ(|s_v,n|)＝f_σ(|s_v,n|)，且

其中参数σ充分小时，f_σ(|s_v,n|)接近L0范数；σ越大，f_σ(|s_v,n|)越平滑；具体采用最速下降法对S进行更新：

其中λ为正定步长，ΔS的每个元素Δs_v,n(v＝1,…,V,n＝1,…,N)满足

Δs_v,n＝exp{θ(s_v,n)}f_σ′(|s_v,n|) (10)

在此，θ(s_v,n)为s_v,n(v＝1,…,V；n＝1,…,N)的相位值，f_σ′(|s_v,n|)为f_σ(|s_v,n|)的一阶导数，满足

第九步：更新被试强度C；根据移不变CPD算法中更新C的方法对C更新；

第十步：计算误差；令iter←iter+1；σ←0.999σ，使σ缓慢变化；根据式(1)，计算本次迭代误差ε_iter，以及相对误差Δε_iter：

Δε_iter＝|(ε_iter-1-ε_iter)/ε_iter-1| (12)

第十一步：若ε_iter小于预设误差阈值ε_{iter_min}，跳转到第十四步，否则执行第十二步；

第十二步：若Δε_iter小于预设相对误差阈值Δε_{iter_min}，跳转到第十四步，否则执行第十三步；

第十三步：若iter大于预设最大迭代次数iter_max，跳转到第十四步，否则执行第三步；

第十四步：对共享SM成分S进行相位消噪；对每个成分n的共享SM成分s_n，进行相位消噪；

第十五步：输出相位消噪的共享SM成分S、共享TC成分B、被试时延

与被试强度C。

Images