CN109316188A - 一种偏头痛脑功能连接模式的提取方法 - Google Patents

Abstract

本发明公开了一种偏头痛脑功能连接模式的提取方法，如下：首先采集偏头痛患者和正常人的大脑静息态功能磁共振数据。之后对采集到的数据进行时间层矫正、头动矫正、标准化和平滑四步预处理。然后根据预处理后的数据，使用自动计算方式获取每个被试的全脑稳态功能连接列向量集。之后采用改进的K均值聚类算法对所有被试的全脑稳态功能连接列向量集进行聚类分析。最后根据阈值提取出偏头痛的脑功能连接模式，从而为后续进一步分析提供基础。该发明有助于根据动态性来获取偏头痛的脑功能连接模式。

Description

一种偏头痛脑功能连接模式的提取方法

技术领域

本发明属于医学影像学图像处理技术领域，尤其是涉及一种基于动态功能连接的用于提取偏头痛脑功能连接模式的方法。

背景技术

偏头痛是一种以疼痛反复发作为临床特征的慢性神经综合征，常伴有恶心、呕吐、光敏感、噪声敏感等临床症状。偏头痛发病率比较高，症状严重影响人们的工作效率和生活质量。偏头痛的发病机制至今尚不完全清楚，目前主要有以下几种假说：皮层扩散性抑制学说、三叉神经血管学说、中枢疼痛调节系统功能异常等。并且偏头痛的诊断都基于医生问诊的方式，缺乏明确的生物标志物体征。随着功能影像学技术的发展，有望借助于功能影像学技术进一步了解偏头痛的发病机制，寻找偏头痛患者影像学上的特征性变化，再进一步结合计算机图像处理和数据挖掘技术，发现偏头痛的神经机制，辅助医生进行诊断和治疗。

在各种医学影像学技术中，功能磁共振是一种可以从功能上揭示大脑神经活动的影像学技术，具有无侵入性、无放射性、高空间分辨率和高时间分辨率等优点。特别是基于血氧水平依赖的功能磁共振(BOLD-fMRI)在临床上应用范围最为广泛。而其中的静息态功能磁共振，能够通过神经元的自发活动对静息态下的脑功能连接进行研究，因此对于进行揭示偏头痛的形成机制和病程发展阶段的相关研究更有意义。

目前研究主要集中于静态的功能连接，而人脑是一个复杂的结构。脑区之间的连接和作用是是随着时间动态变化的。当疾病发生时，其动态变化的异常会体现在相应的脑功能连接上。本技术发明旨在基于动态功能连接，通过一定的算法来提取出偏头痛的脑功能连接模式，在此基础上可以进一步研究偏头痛脑功能连接模式特征，从而为探索偏头痛的神经机制异常提供基础。

发明内容

本发明的目的在于提供一种提取偏头痛脑功能连接模式的方法。通过动态功能连接算法提取出全脑稳态功能连接列向量集，使用改进的K均值聚类算法对全脑稳态功能连接列向量集进行聚类。之后根据阈值，提取出偏头痛的脑功能连接模式。为进一步分析提供数据基础。

为了达到上述目的，本发明通过以下技术方案实现：

一种偏头痛脑功能连接模式的提取方法，其特点是，该提取方法包含如下步骤：

步骤1，采集多个偏头痛患者的大脑静息态功能磁共振数据和多个正常人的大脑静息态功能磁共振数据。

步骤2，对采集到的每个被试的静息态功能磁共振数据进行预处理操作，包括时间层矫正、头动矫正、标准化和平滑四个操作。

步骤3，根据上述预处理后的每个被试的静息态功能磁共振数据，使用自动计算方式获得每个被试的全脑稳态功能连接列向量集。

步骤4，使用改进的K均值聚类算法对所有被试的全脑稳态功能连接列向量集进行聚类分析。

步骤5，根据上述聚类分析结果，根据阈值提取偏头痛的脑功能连接模式。

所述的偏头痛脑功能连接模式的提取方法，其特点是，所述的步骤3包含如下步骤：

步骤3.1，对于每一个被试，计算所有脑区的平均时间序列，具体方法如下：根据预处理后的静息态功能磁共振数据，使用现有脑区模板获取全脑的N个脑区B₁、B₂、……B_N所包含的体素的时间信号值，将每个脑区包含的时间信号值进行平均，形成N个脑区的平均时间序列T₁、T₂、……T_N。对于时间长度为M的被试来说，平均时间序列其中 表示第n个脑区在i时刻所有体素的平均时间信号值。

步骤3.2，对于每一个被试，计算全脑动态功能连接矩阵组，具体方法如下：通过采用滑动窗方法，利用一个特定宽度W的滑动窗以步长为1在时间序列T₁、T₂……T_N上滑动，在第j个时间窗下第n个脑区的时间序列表示为 计算该被试所有脑区两两时间序列之间的皮尔森相关系数，得到M-W+1个动态功能连接矩阵dFC，这些动态功能连接矩阵构成了该被试的全脑动态功能连接矩阵组SetdFC＝{dFC₁,dFC₂,…,dFC_b,…,dFC_M-W+1}。

进一步地，脑区两两时间序列之间的皮尔森相关系数具体是指第a个滑动窗下脑区B_x和B_y时间序列dT_xa与dT_ya之间的皮尔森相关系数，其公式如下：

式中Cov(dT_xa,dT_ya)为dT_xa、dT_ya的协方差，D(dT_xa)、D(dT_ya)分别为dT_xa、dT_ya的方差，1≤a≤(M-W+1)，1≤x≤N，1≤y≤N。

进一步地，动态功能连接矩阵dFC_b是指第b个滑动窗下所有脑区两两时间序列之间的皮尔森相关系数组成的动态功能连接矩阵，具体表示为：

其中1≤b≤(M-W+1)，1≤u≤N，1≤v≤N。

步骤3.3，对于每一个被试，计算全脑动态功能连接强度矩阵wdFC，具体方法如下：将全脑动态功能连接矩阵组SetdFC的动态功能连接矩阵dFC_b(1≤b≤(M-W+1))，按照行将每行的数值累加后放到第一列，得到一个全脑动态功能连接强度列向量wdFCV_b(1≤b≤(M-W+1))；每个列向量大小为N×1；将M-W+1个列向量按照时间点从小到大进行拼接，形成全脑动态功能连接强度矩阵wdFC＝[wdFCV₁,wdFCV₂,…,wdFCV_b,…,wdFCV_M-W+1]，大小为N×(M-W+1)。

进一步地，wdFCV_b是指第b个滑动窗下的全脑动态功能连接强度列向量，具体表示为：

其中1≤b≤(M-W+1)，1≤c≤N，1≤d≤N。

步骤3.4，使用自动计算的方式将每一个被试的全脑动态功能连接强度矩阵进行自动分割，具体方法如下：对于每一个被试的全脑动态功能连接强度矩阵wdFC，依次计算wdFCV_k(2≤k≤N-1)与wdFCV_k-1的欧式距离记为Dist1，计算wdFCV_k与wdFCV_k+1记为Dist2；如果Dist2小于Dist1，则将k记录下来，标记为分割点；起始列向量和结尾列向量不做计算。

进一步地，列向量之间欧式距离，其公式如下：

其中E、F为列向量，e_g为E的元素，f_g为F的元素，H为列向量大小。

步骤3.5，对于每一个被试，计算全脑稳态功能连接列向量集，具体方法如下：根据得出的分割点，将动态功能连接强度矩阵wdFC沿着时间轴方向分割成段，对每个段内的值进行平均即可得到全脑稳态功能连接列向量；每个列向量大小为N×1；将所有列向量构成了该被试的全脑稳态功能连接列向量集。

上述的偏头痛脑功能连接模式的提取方法，其特点是，上述的步骤4包含如下步骤：

步骤4.1，将所有被试的全脑稳态功能连接列向量集进行合并构成聚类样本，每一个样本是一个被试的一个全脑稳态功能连接列向量；采用基于自动目标生成过程(ATGP)的方法生成Q个初始类心。

步骤4.2，根据上步的初始类心，采用K均值聚类算法对全脑稳态功能连接列向量样本进行聚类，得到Q个类别。

上述的偏头痛脑功能连接模式的提取方法，其特点是，上述的步骤5包含如下步骤：

步骤5.1，按照每个类别，计算属于偏头痛的样本数目与该类别所有样本数目的比例Propt1和属于正常人的样本数目与该类别所有样本数目的比例Propt2；计算R＝Propt1/(Propt1+Propt2)；如果R大于等于某一阈值，则将该类别标记出来。

步骤5.2，从标记的类别中提取出偏头痛的脑功能连接模式，具体方法如下：对于标记的类别，提取出属于偏头痛的所有样本，也就是属于偏头痛的全脑稳态功能连接列向量；之后逆向定位出步骤3.2的动态功能连接矩阵dFC，将定位出的所有动态功能连接矩阵dFC求平均，得出该类别下偏头痛的脑功能连接模式BFC，具体表示如下：

其中Numb指的是本类别中包含偏头痛的样本数。

本发明与现有技术相比具有以下优点：本发明引入动态功能连接的方法，有助于根据动态性来获取偏头痛的脑功能连接模式；本发明使用了自动分割的方法，提高了整个方法的准确度和效率；本发明提取的偏头痛的脑功能连接模式可以为临床上偏头痛疾病的进一步研究和防治提供分析基础。

附图说明

图1为本发明一种偏头痛脑功能连接模式的提取方法的整体流程图。

具体实施方式

以下结合附图，通过详细说明一个较佳的具体实施例，对本发明做进一步阐述。

如图1所示一种偏头痛脑功能连接模式的提取方法，该提取方法包含如下步骤：

步骤1，采集多个偏头痛患者的大脑静息态功能磁共振数据和多个正常人的大脑静息态功能磁共振数据。

按照临床医学诊断标准，分别采集偏头痛患者组(MP)和正常人对照组(NC)的大脑静息态功能磁共振图像，两组样本的数目各35个，一共70个被试。数据获取过程中被试要求保持大脑清醒，平躺于磁共振仪器内。每个被试的时间点数为160。

步骤2，对采集到的每个被试的静息态功能磁共振数据进行预处理操作，包括时间层矫正、头动矫正、标准化和平滑四个操作。

所有的数据集采用SPM8软件完成时间层矫正、头动矫正、标准化和平滑。

步骤3，根据上述预处理后的每个被试的静息态功能磁共振数据，使用自动计算方式获得每个被试的全脑稳态功能连接列向量集。

步骤3.1，对于每一个被试，计算所有脑区的平均时间序列，具体方法如下：根据预处理后的静息态功能磁共振数据，使用Brainnetome Atlas获取全脑的N＝246个脑区B₁、B₂、……B₂₄₆所包含的体素的时间信号值，将每个脑区包含的时间信号值进行平均，形成246个脑区的平均时间序列T₁、T₂、……T₂₄₆。对于时间长度为M＝160的被试来说，平均时间序列 其中表示第n个脑区在i时刻所有体素的平均时间信号值。

步骤3.2，对于每一个被试，计算全脑动态功能连接矩阵组，具体方法如下：通过采用滑动窗方法，利用一个特定宽度W＝12的滑动窗以步长为1在时间序列T₁、T₂……T₂₄₆上滑动，在第j个时间窗下第n个脑区的时间序列表示为

计算该被试所有脑区两两时间序列之间的皮尔森相关系数，得到149个动态功能连接矩阵dFC，这些动态功能连接矩阵构成了该被试的全脑动态功能连接矩阵组SetdFC＝{dFC₁,dFC₂,…,dFC_b,…,dFC₁₄₉}。

进一步地，脑区两两时间序列之间的皮尔森相关系数具体是指第a个滑动窗下脑区B_x和B_y时间序列dT_xa与dT_ya之间的皮尔森相关系数，其公式如下：

式中Cov(dT_xa,dT_ya)为dT_xa、dT_ya的协方差，D(dT_xa)、D(dT_ya)分别为dT_xa、dT_ya的方差，1≤a≤149，1≤x≤246，1≤y≤246。

进一步地，动态功能连接矩阵dFC_b是指第b个滑动窗下所有脑区两两时间序列之间的皮尔森相关系数组成的动态功能连接矩阵，具体表示为：

其中1≤b≤149，1≤u≤246，1≤v≤246，N＝246。

步骤3.3，对于每一个被试，计算全脑动态功能连接强度矩阵wdFC，具体方法如下：将全脑动态功能连接矩阵组SetdFC的动态功能连接矩阵dFC_b(1≤b≤149)，按照行将每行的数值累加后放到第一列，得到一个全脑动态功能连接强度列向量wdFCV_b(1≤b≤149)；每个列向量大小为246×1；将149个列向量按照时间点从小到大进行拼接，形成全脑动态功能连接强度矩阵wdFC＝[wdFCV₁,wdFCV₂,…,wdFCV_b,…,wdFCV₁₄₉]，大小为246×149。

进一步地，wdFCV_b是指第b个滑动窗下的全脑动态功能连接强度列向量，具体表示为：

其中1≤b≤149，1≤c≤246，1≤d≤246，N＝246。

步骤3.4，使用自动计算的方式将每一个被试的全脑动态功能连接强度矩阵进行自动分割，具体方法如下：对于每一个被试的全脑动态功能连接强度矩阵wdFC，依次计算wdFCV_k(2≤k≤(N-1＝245))与wdFCV_k-1的欧式距离记为Dist1，计算wdFCV_k与wdFCV_k+1记为Dist2；如果Dist2小于Dist1，则将k记录下来，标记为分割点；起始列向量和结尾列向量不做计算。

进一步地，列向量之间欧式距离，其公式如下：

其中E、F为列向量，e_g为E的元素，f_g为F的元素，H为列向量大小。在此实施例中H为246。

步骤3.5，对于每一个被试，计算全脑稳态功能连接列向量集，具体方法如下：根据得出的分割点，将动态功能连接强度矩阵wdFC沿着时间轴方向分割成段，对每个段内的值进行平均即可得到全脑稳态功能连接列向量；每个列向量大小为246×1；将所有列向量构成了该被试的全脑稳态功能连接列向量集。

步骤4，使用改进的K均值聚类算法对所有被试的全脑稳态功能连接列向量集进行聚类分析。

步骤4.1，将所有被试的全脑稳态功能连接列向量集进行合并构成聚类样本，每一个样本是一个被试的一个全脑稳态功能连接列向量；采用基于自动目标生成过程(ATGP)的方法生成20个初始类心。

步骤4.2，根据上步的初始类心，采用K均值聚类算法对全脑稳态功能连接列向量样本进行聚类，得到20个类别。

步骤5，根据上述聚类分析结果，根据阈值提取偏头痛的脑功能连接模式。

步骤5.1，按照每个类别，计算属于偏头痛的样本数目与该类别所有样本数目的比例Propt1和属于正常人的样本数目与该类别所有样本数目的比例Propt2；计算R＝Propt1/(Propt1+Propt2)；如果R大于等于某一阈值，则将该类别标记出来；本实施例中阈值定为0.8，将符合判定规则的类别标记出来。

步骤5.2，从标记的类别中提取出偏头痛的脑功能连接模式，具体方法如下：对于标记的类别，提取出属于偏头痛的所有样本，也就是属于偏头痛的全脑稳态功能连接列向量；之后逆向定位出步骤3.2的动态功能连接矩阵dFC，将定位出的所有动态功能连接矩阵dFC求平均，得出该类别下偏头痛的脑功能连接模式BFC，具体表示如下：

其中Numb指的是本类别中包含偏头痛的样本数。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (4)

1.一种偏头痛脑功能连接模式的提取方法，其特征在于，该提取方法包含如下步骤：

步骤1，采集多个偏头痛患者的大脑静息态功能磁共振数据和多个正常人的大脑静息态功能磁共振数据；

步骤2，对采集到的每个被试的静息态功能磁共振数据进行预处理操作，包括时间层矫正、头动矫正、标准化和平滑四个操作；

步骤3，根据上述预处理后的每个被试的静息态功能磁共振数据，使用自动计算方式获得每个被试的全脑稳态功能连接列向量集；

步骤4，使用改进的K均值聚类算法对所有被试的全脑稳态功能连接列向量集进行聚类分析；

步骤5，根据上述聚类分析结果，根据阈值提取偏头痛的脑功能连接模式。

2.如权利要求1所述的一种偏头痛脑功能连接模式的提取方法，其特征在于，所述的步骤3包含如下步骤：

步骤3.1，对于每一个被试，计算所有脑区的平均时间序列，具体方法如下：根据预处理后的静息态功能磁共振数据，使用现有脑区模板获取全脑的N个脑区B₁、B₂、……B_N所包含的体素的时间信号值，将每个脑区包含的时间信号值进行平均，形成N个脑区的平均时间序列T₁、T₂、……T_N；对于时间长度为M的被试来说，平均时间序列其中 表示第n个脑区在i时刻所有体素的平均时间信号值；

步骤3.2，对于每一个被试，计算全脑动态功能连接矩阵组，具体方法如下：通过采用滑动窗方法，利用一个特定宽度W的滑动窗以步长为1在时间序列T₁、T₂……T_N上滑动，在第j个时间窗下第n个脑区的时间序列表示为 计算该被试所有脑区两两时间序列之间的皮尔森相关系数，得到M-W+1个动态功能连接矩阵dFC，这些动态功能连接矩阵构成了该被试的全脑动态功能连接矩阵组SetdFC＝{dFC₁，dFC₂，...，dFC_b，...，dFC_M-W+1}；

进一步地，脑区两两时间序列之间的皮尔森相关系数具体是指第a个滑动窗下脑区B_x和B_y时间序列dT_xa与dT_ya之间的皮尔森相关系数，其公式如下：

式中Cov(dT_xa，dT_ya)为dT_xa、dT_ya的协方差，D(dT_xa)、D(dT_ya)分别为dT_xa、dT_ya的方差，1≤a≤(M-W+1)，1≤x≤N，1≤y≤N；

进一步地，动态功能连接矩阵dFC_b是指第b个滑动窗下所有脑区两两时间序列之间的皮尔森相关系数组成的动态功能连接矩阵，具体表示为：

其中1≤b≤(M-W+1)，1≤u≤N，1≤v≤N；

步骤3.3，对于每一个被试，计算全脑动态功能连接强度矩阵wdFC，具体方法如下：将全脑动态功能连接矩阵组SetdFC的动态功能连接矩阵dFC_b(1≤b≤(M-W+1))，按照行将每行的数值累加后放到第一列，得到一个全脑动态功能连接强度列向量wdFCV_b(1≤b≤(M-W+1))；每个列向量大小为N×1；将M-W+1个列向量按照时间点从小到大进行拼接，形成全脑动态功能连接强度矩阵wdFC＝[wdFCV₁，wdFCV₂，…，wdFCV_b，…，wdFCV_M-W+1]，大小为N×(M-W+1)；

进一步地，wdFCV_b是指第b个滑动窗下的全脑动态功能连接强度列向量，具体表示为：

其中1≤b≤(M-W+1)，1≤c≤N，1≤d≤N；

步骤3.4，使用自动计算的方式将每一个被试的全脑动态功能连接强度矩阵进行自动分割，具体方法如下：对于每一个被试的全脑动态功能连接强度矩阵wdFC，依次计算wdFCV_k(2≤k≤N-1)与wdFCV_k-1的欧式距离记为Dist1，计算wdFCV_k与wdFCV_k+1记为Dist2；如果Dist2小于Dist1，则将k记录下来，标记为分割点；起始列向量和结尾列向量不做计算；

进一步地，列向量之间欧式距离，其公式如下：

其中E、F为列向量，e_g为E的元素，f_g为F的元素，H为列向量大小；

步骤3.5，对于每一个被试，计算全脑稳态功能连接列向量集，具体方法如下：根据得出的分割点，将动态功能连接强度矩阵wdFC沿着时间轴方向分割成段，对每个段内的值进行平均即可得到全脑稳态功能连接列向量；每个列向量大小为N×1；将所有列向量构成了该被试的全脑稳态功能连接列向量集。

3.如权利要求2所述的偏头痛脑功能连接模式的提取方法，其特征在于，所述的步骤4包含如下步骤：

步骤4.1，将所有被试的全脑稳态功能连接列向量集进行合并构成聚类样本，每一个样本是一个被试的一个全脑稳态功能连接列向量；采用基于自动目标生成过程(ATGP)的方法生成Q个初始类心；

步骤4.2，根据上步的初始类心，采用K均值聚类算法对全脑稳态功能连接列向量样本进行聚类，得到Q个类别。

4.如权利要求3所述的偏头痛脑功能连接模式的提取方法，其特征在于，所述的步骤5包含如下步骤：

步骤5.1，按照每个类别，计算属于偏头痛的样本数目与该类别所有样本数目的比例Propt1和属于正常人的样本数目与该类别所有样本数目的比例Propt2；计算R＝＝Propt1/(Propt1+Propt2)；如果R大于等于某一阈值，则将该类别标记出来；

步骤5.2，从标记的类别中提取出偏头痛的脑功能连接模式，具体方法如下：对于标记的类别，提取出属于偏头痛的所有样本，也就是属于偏头痛的全脑稳态功能连接列向量；之后逆向定位出步骤3.2的动态功能连接矩阵dFC，将定位出的所有动态功能连接矩阵dFC求平均，得出该类别下偏头痛的脑功能连接模式BFC，具体表示如下：

其中Numb指的是本类别中包含偏头痛的样本数。