CN104434109A - 一种功能核磁共振时间序列匹配方法 - Google Patents

Abstract

本发明公开了一种功能核磁共振时间序列匹配方法，包括以下步骤：计算系数矩阵β的初始值；计算观测信号时间序列Y与设计时间序列X间的范式距离F_dist；随机获取新的记录点；计算新的范式距离F_dist；输出β，完成Y与X的β匹配。本发明在fMRI时间序列匹配中将时域信号经快速离散傅立叶模变换后得到频域序列，从而完全消除相位信息，以达到消除要匹配的fMRI时间序列间相位差的目的。本发明与当前相位校正的方法相比，更为简单，消耗的计算量少。本发明对离散傅里叶变换进行加权约束，以降低高频部分的影响，优先考虑了低频有效信号，更为明确地确定“最感兴趣频率”的位置，从而降低激活检验结果出现“伪阳性”体素的可能性。

Description

一种功能核磁共振时间序列匹配方法

技术领域

本发明涉及一种数据处理技术，特别是一种功能核磁共振时间序列同频不同相的匹配方法。

背景技术

功能磁共振(fMRI)技术，是无损地研究人脑行为活动的重要影像技术，具有广泛而重要的应用前景。fMRI信号是一种难以直观解读的观测型信号，因此针对fMRI信号的合理分析具有重要意义，这既是热点但也是难点问题。

fMRI技术的通过大脑的血氧水平依赖效应呈现大脑的活动情况，通过该效应，可以实时地检测大脑活动的脑血流变化情况。fMRI数据有以下几个特点：

1)延迟特性。一个脑功能刺激，fMRI信号往往并不能做出及时的反馈，而是在几秒后达到峰值。

2)混淆特性。在以往的信号分析研究中，fMRI信号通常被看作一种由不同刺激经过线性叠加而形成的观测信号。

3)采样点多。fMRI的扫描精度比较高，若以2秒每帧的速度进行扫描，可以获得64*64*32分辨率的大脑图像。因此对维数比较敏感的算法，在应用到fMRI数据中往往会遇到困难。

3)信噪比低。大脑是一个复杂的整体，在进行某些特定认知任务的时候，大脑还肩负着人体的调节工作，比如呼吸、心跳等一系列维持生命所必须的生理活动。而在进行有意识的活动的同时，大脑也会出现许多伴随着诸多因长期的生活训练而出现的无意识的活动。

fMRI实验数据按照实验刺激的不同主要分为静息态数据和功能实验数据。静息态数据无需被扫描者做出任何活动，主要用于观察人处于无意识状态下的大脑活动。而功能数据则是在进行脑科学实验时，实验设计者给予受试者外界刺激，诱导受试者有意识地进行思维活动所产生的fMRI数据，主要用于观察当人处于有意识的行为时，大脑的活动状态。

目前用于分析fMRI行为实验数据的常用软件主要有统计参数映射模型(SPM)、fMRI组独立成分工具箱(GIFT)、静息态数据处理方法(Dparsf)等，这些工具已经用于进行人脑行为认知的研究当中，最为广泛的是SPM(StatisticalParametric Maps)，它可以提取特定刺激因素引起大脑活动的范围，即脑激活区。SPM主要利用刺激时间序列这一先验条件，提出一种线性叠加模型，即广义线性模型(GLM)，来分析不同脑区的fMRI时间序列信号与刺激时间序列的线性叠加关系，从而判定该区域与特定刺激的相关性。SPM框架包括多个组成部分，包括时间校正、头动校正、标准化预处理等数据预处理方法，并集成GLM、pair&no pair T检验、F检验和贝叶斯等多种数据分析方法。Rest、Dparsf、GIFT、MICA等fMRI数据分析工具包直接集成或加载SPM的部分功能，或要求数据由SPM进行预处理后再进行数据分析。

SPM虽然被广泛应用，但也有部分fMRI信号并没有得到很好地转换为脑区激活信息，因为其依赖的广义线性模型(SPM-GLM)对刺激时间序列的要求比较苛刻，这大大限制了fMRI实验设计的灵活性，而且该方法将fMRI时间序列视为一个向量，通过评价拟合参考序列与待测时间序列信号的范数距离作为目标函数，即待测fMRI信号与模型参考时间序列的匹配程度，来评价优化的目的。但是作为时间序列，相邻的时间点间联系是很紧密的，而fMRI信号作为时间序列，由刺激引发的信号发生延迟的个体差异较大，即使同一个体对不同的刺激反应延迟也存在差异，因此当脑神经元活动的时间与刺激呈现的时间存在较大不确定性的差异时，即大脑活动的时间序列本身难以准确描述的时候，传统的SPM-GLM方法很难得到理想的结果。

发明内容

为解决现有技术存在的上述问题，本发明要提出一种简便的、可以克服延迟的功能核磁共振时间序列匹配方法。

为了实现上述目的，本发明的技术方案如下：一种功能核磁共振时间序列匹配方法，包括以下步骤：

A、计算系数矩阵β的初始值

将广义线性模型GLM表述为如下形式：

Y＝Xβ+ε..............(1)

其中，Y为观测信号时间序列，X为设计时间序列，β为系数矩阵，ε为误差，则系数矩阵β的初始值的计算公式如下：

β＝(X^TX)^-1X^Ty………(2)

式中y是Y的一个分量。

B、计算观测信号时间序列Y与设计时间序列X间的范式距离F_dist:

F_dist(Y,Xβ)＝norm(Y-Xβ)*norm(|fft(Y)|-|fft(Xβ)|*w).....(3)

其中，norm()表示二范数，fft()表示快速傅立叶变换，|·|表示取绝对值，w为权重向量，其j分量w(j)由下列公式计算：

w(j)＝1/(e^-ηj-e^θj)...........(4)

w＝w/max(w)............(5)

其中η为正向感兴趣频率常数，θ为负向感兴趣频率常数，它们共同通过公式(5)来确定“最感兴趣频率”的位置。

C、设置初始值

定义初始记录器b_old＝b_min＝β，设置ε＝0.1，Δ＝1，设置迭代计算器t＝0，s_j＝F_dist(y,x*β)，这里j＝1。

D、设置随机方向正确次数记录器

设置随机方向正确次数记录器k＝0。

E、随机获取新的记录点

更新步长Δ＝Δ/2，在当前记录器bold邻域内，从N^m(b_old,Δ)分布中随机获取新的记录点b_new，使得b_new满足(b_min-b_old)*(b_new-b_old)+ε>0。

F、计算新的范式距离F_dist

s_new＝F_dist(y,x*b_new)............(6)

其中x是X的一个分量。

G、如果s_new<s_j，则k增加1位。

H、如果k>v*p，返回步骤D。这里v是采样数，p是0～0.5之间的随机数。

I、更新s_j+1＝s_new；β_j+1＝b_min；j增加1位。

J、转步骤D，迭代直到j到达迭代最大次数，所述的最大次数为90-100次。

K、输出β，完成Y与X的β匹配。

与现有技术相比，本发明具有以下有益效果：

1、频域模变换：如公式(3)所示，本发明在fMRI时间序列匹配中将时域信号经快速离散傅立叶模变换后得到频域序列，从而完全消除相位信息，以达到消除要匹配的fMRI时间序列间相位差的目的。本发明与当前相位校正的方法相比，更为简单，消耗的计算量少。

2、加权约束：考虑离散傅里叶变换，高频部分的离散傅里叶变换相当于高频噪声，本发明对离散傅里叶变换进行加权约束，如公式(4)-(5)所示，以降低高频部分的影响，优先考虑了低频有效信号，更为明确地确定“最感兴趣频率”的位置，从而降低激活检验结果出现“伪阳性”体素的可能性。

3、变步长启发式寻优：在fMRI中，若将GLM看作一个时间序列的匹配问题，欲寻求其最优解则无疑是一种连续空间中的问题。如步骤D到步骤I，本发明通过一个随机给定的初始点，逐步寻求其邻域内的更好的点，并利用近似两个互相穿过对方球心的球的球面所围成的一个空间形成的一个二阶连续空间作为空间禁忌启发式准则，提高了搜索效率，当逐步寻优遇到瓶颈时，逐步减少搜索半径，提高了搜索精度。

附图说明

本发明共有附图4张，其中：

图1是非线性加权向量曲线图。

图2是本发明的寻优流程图。

图3是模拟的时间序列匹配图。

图4是所得的感兴区的激活图。

具体实施方式

下面结合附图对本发明进行进一步地描述。

图1是非线性加权向量曲线图，其中，实线是观测时间序列变换到频域后的位置曲线，虚线是频域匹配非线性加权权重w。在步骤B中，为fMRI时间序列匹配提供频域匹配权重w，也就是图1中虚线所示，提高信噪比，降低激活检验结果出现“伪阳性”体素的可能性。

图2是本发明的寻优流程图，该流程图是在时间序列频域对齐之后通过迭代方式降低匹配误差的过程，如步骤D到J所描述，该寻优方式是在二阶禁忌空间中进行随机邻域优化，使得本发明中的匹配可以稳定收敛到全局最优解，同时避免陷入局部最优解，而且如步骤E所描述，寻优步长Δ可变，提高了本发明的精度，使得最终得到的β解可以稳定在精度范围内，并且根据Abel收敛准则，该收敛于全局最优解。

图3是模拟的时间序列匹配图，其中实线是设计时间序列X，带小方块标记的实线是由本发明进行匹配所得到的fMRI信号时间序列，显示了本发明的fMRI时间序列匹配效果。

图4是所得的感兴区的激活图，该图是经过本发明中步骤A到K完成记忆想象实验观测时间序列与设计时间序列匹配后对脑激活的判别结果，从图4可以看出激活体素准确。

表1是激活体素数目统计结果，从比较结果看，经过本发明的匹配方法检测激活体素，有效降低“伪阳性”。

表1：激活体素数目统计结果表

参考时间序列	GLM	本发明
			序列1	195	89
序列2	206	54

Claims (1)

1.一种功能核磁共振时间序列匹配方法，其特征在于：包括以下步骤：

A、计算系数矩阵β的初始值

将广义线性模型GLM表述为如下形式：

Y＝Xβ+ε……………(1)

其中，Y为观测信号时间序列，X为设计时间序列，β为系数矩阵，ε为误差，则系数矩阵β的初始值的计算公式如下：

β＝(X^TX)^-1X^Ty………(2)

式中y是Y的一个分量；

B、计算观测信号时间序列Y与设计时间序列X间的范式距离F_dist:

F_dist(Y,Xβ)＝norm(Y-Xβ)*norm(|fft(Y)|-|fft(Xβ)|*w)…..(3)

其中，norm()表示二范数，fft()表示快速傅立叶变换，|·|表示取绝对值，w为权重向量，其j分量w(j)由下列公式计算：

w(j)＝1/(e^-ηj-e^θj)………..(4)

w＝w/max(w)…………(5)

其中η为正向感兴趣频率常数，θ为负向感兴趣频率常数，它们共同通过公式(5)来确定“最感兴趣频率”的位置；

C、设置初始值

定义初始记录器b_old＝b_min＝β，设置ε＝0.1，Δ＝1，设置迭代计算器t＝0，s_j＝F_dist(y,x*β)，这里j＝1；

D、设置随机方向正确次数记录器

设置随机方向正确次数记录器k＝0；

E、随机获取新的记录点

更新步长Δ＝Δ/2，在当前记录器bold邻域内，从N^m(b_old,Δ)分布中随机获取新的记录点b_new，使得b_new满足(b_min-b_old)*(b_new-b_old)+ε>0；

F、计算新的范式距离F_dist

s_new＝F_dist(y,x*b_new)…………(6)

其中x是X的一个分量；

G、如果s_new<s_j，则k增加1位；

H、如果k>v*p，返回步骤D；这里v是采样数，p是0～0.5之间的随机数；

I、更新s_j+1＝s_new；β_j+1＝b_min；j增加1位；

J、转步骤D，迭代直到j到达迭代最大次数，所述的最大次数为90-100次；

K、输出β，完成Y与X的β匹配。