CN106204581A - 基于pca与k均值聚类的动态脑功能连接模式分解方法 - Google Patents

Abstract

基于PCA与K均值聚类的动态脑功能连接模式分解方法，是将动态功能连接分解为一些基本连接模式。先对静息态磁共振数据预处理、再提取脑区的时间序列，然后采用滑动窗方法计算脑区两两之间的相关系数，构造相应的功能连接矩阵，最后对功能连接矩阵进行PCA与K均值聚类分析得到基本连接模式。本发明提出的动态功能连接模式分解算法，能够有效准确地得到分解后的基本功能连接模式，不仅能够捕捉不同基本连接模式之间的转变，并且为临床上精神分裂症、阿尔兹海默症、抑郁症等神经精神疾病的进一步研究和防治提供策略。

Description

基于PCA与K均值聚类的动态脑功能连接模式分解方法

技术领域

本发明属于图像处理领域，具体涉及一种基于PCA与K均值聚类的动态脑功能连接模式分解方法，尤其涉及使用主成分分析(Principal,Component Aanlysis,PCA)和K均值聚类算法对默认模式网络、执行控制网络以及感觉网络的动态脑功能连接模式分解。

背景技术

fMRI是研究脑活动、脑功能的主要的无创方法之一，具有毫米级的空间分辨率，已成为神经科学探索人类大脑神经机制的重要工具。fMRI基于血氧水平依赖对比度(BOLD)间接测量神经元活动，它通过测量由神经活动引起的脑血流和脑血氧等成分变化而造成的磁共振信号变化来反应脑活动。BOLD-fMRI方法的发展对脑认知功能的研究具有突破性的进展。

静息态fMRI刻画了大脑在没有任务时的神经基准活动，具有非常重要的生理意义。研究表明在静息态下，仍存在BOLD信号的波动，在这些自发振荡的信号中，存在某种特定的功能连接，构成了静息态功能连接网络，反映了静息态下人脑的活动。有研究表明，神经精神疾病如阿尔茨海默、轻度认知受损(mild cognitive impairment)、抑郁症(depression)以及精神分裂症(schizophrenia)等的静息态功能 网络连接存在异常，因此研究静息态功能网络连接的异常可以为神经精神疾病的研究提供有效和可靠的生物学指标。静息态功能连接网络主要包含默认模式网络、感觉网络(听觉、视觉和运动)和执行控制网络等，本发明主要针对以上三个网络进行研究。

目前大多研究只关注静态功能连接，对于静息态脑功能成像静止性的假设可能忽略了重要的动态信息。事实上，大脑本身是一个复杂的动态系统，是高度变化的，并且大脑的功能连接具有瞬时变化特性，其包含更丰富的组织信息。研究大脑功能连接的时间动态变化信息，有助于对大脑功能组织结构更综合和全面的认识。

发明内容

为了达到上述目的，本发明的目的在于提供一种基于PCA与K均值聚类的动态脑功能连接模式分解方法，一方面能够捕捉不同基本连接模式之间的转变，其转变可能与神经认知的变化有关；另一方面，对于仅存在于部分连接模式中的，采用静态功能连接难以预测，因此，动态功能连接模式分解也为生理病理学机制研究提供新的思路。

为了达到上述目的，本发明的技术方案为：

基于PCA与K均值聚类的动态功能连接模式分解方法，具体步骤如下：

(1)、对被试进行静息态磁共振数据采集，并进行预处理，磁共振数据的预处理目的是提高脑功能图像的信噪比，并将被试图像与标准模板进行配准变换；

(2)、预处理之后，提取DMN、ECN以及SN网络的时间序列， 采用种子点法提取DMN23个子区的时间序列、ECN12个子区的时间序列、SN13个子区的时间序列，种子点区域选取半径为6mm的球形区域，选取种子点之后，并对其中所有体素的信号进行加权平均，进而得到种子点的时间序列曲线；

(3)、基于PCA进行动态功能连接模式分解，首先通过采用滑动窗方法计算每个被试所有脑区两两之间的皮尔森相关系数，进而得到功能连接矩阵，功能连接矩阵的大小为M*N，其中M为功能连接对的个数，N为滑动窗的个数与被试个数的乘积；其次，采用PCA方法对功能连接矩阵进行模式分解，得到矩阵特征向量，按从大到小顺序保留最上面的k个特征向量，k个特征向量表示动态功能连接的k个主要连接模式；

(4)、为研究某些连接模式可能是准稳定的，即它们随时间重复出现以及存在于众多被试中，基于K均值聚类方法对相应的功能连接矩阵进行聚类，每个聚类中心代表存在的一种连接模式，实现基于K均值聚类的动态功能连接模式分解。

本发明的创新点在于：提出基于PCA与K均值聚类的动态功能连接模式分解方法，一方面能够捕捉不同基本连接模式之间的转变；另一方面，基本连接模式能够预测一些神经精神疾病(如阿尔兹海默症、精神分裂症、抑郁症等)，为生理病理学机制研究提供新的思路。

附图说明

图1是静息态功能图像预处理流程图。

图2是头颅分割前后对比图。

图3是功能像配准图，左侧为标准模板，右侧为配准后的图像。

图4是基于PCA的动态功能连接模式分解流程图。

图5是对DMN、ECN以及SN网络的PCA模式分解结果。

图6是对三个网络总体进行PCA模式分解的结果图。

图7是基于K均值聚类的动态功能连接模式分解流程图。

图8是对三个网络总体进行K均值聚类模式分解的结果图。

图9是其中3个被试的状态分配图。

图10是各脑区种子点坐标。

图11是所有被试的平均转移矩阵。

具体实施方式

下面结合附图对本发明做详细叙述。

本发明基于PCA与K均值聚类对静息态动态功能连接进行功能连接模式分解。

1、首先对采集的原始静息态磁共振数据进行预处理，由于磁共振扫描过程中各种各样的噪声的影响，个体自身存在尺度和位置上的差异，非常有必要在分析数据之前对数据做一定的预处理。在整个的实验的数据获取中，主要的噪声信息来源包括：(1)物理头动；(2)图像内层间扫描时间差别；(3)外在磁场的不均匀性。脑功能图像预处理是在保留脑功能图像细节的同时，使用脑功能图像与标准模板进行仿射配准变换方式的预处理，并提高脑功能图像的信噪比。

功能磁共振数据的预处理选取了Linux系统Ubuntu12.04下的AFNI和FSL软件，并通过编写Batch文件实现数据的批处理。预处理 流程见图1，主要包含以下几个方面：

1)结构像头颅分割

采集得到的结构像通常包含头颅信息，需要进行头颅分割，从而消除颅骨部位引入的伪影对后续数据分析的影响。基于FSL软件提供的3drefit、3dresample和fastsegment实现头部颅骨和大脑内部组织结构的分割，分割结果如图2所示。

2)时间对齐

血液动力学函数表明血液对刺激的响应有一定的时间延迟，由于一个TR期间采集到全脑图像，这导致每一层的图像并非在同一时刻采集，而是发生在整个扫描时间段内，时间校正就是通过类似于插值的方法对每一层图像进行层时间处理，使得一个TR周期内各层图像近似在同一时刻获取。输入每帧图像的各扫描参数进行时间校正，包括扫描时间、扫描层数、第一层到最后一层的时间间隔以及扫描顺序，以中间层为参考进行时间校正，将每帧图像校正为几乎是同一时间点采集的各层图像。

3)头动校正

为避免扫描过程中发生微小头动，对扫描质量及结果造成影响，以每个被试的第一帧图像为基准，将其余的所有图像和基准图像对齐。进行头动校正时一般将被试的大脑看作一个刚体，因此在fMRI实验中被试头部的运动可以近似成一种刚体运动，即只有平移变换与旋转变换的组合。通过AFNI的3dvolreg函数使得其余的所有图像与参考图像配准，如果头动超过一个体素则去除该被试。

4)空间平滑

对图像进行高斯平滑，主要目的是提高图像信噪比以及确保图像数据具有随机高斯场性质，这对于需要做平均以及统计结果非常重要。空间平滑采用高斯函数进行高斯平滑，能有效地消弱随机噪声对fMRI信号的影响，提高数据的信噪比。三维高斯函数是比较常用的空间平滑方法，其半高全宽决定了空间平滑的力度，本文选用半高全宽(Full Widthat Half Maximum,FWHM)为6mm的高斯核函数进行数据平滑。

5)时域带通滤波

静息态fMRI信号的低频波动反映了自发的神经活动，因此采用带通为0.01-0.1Hz的带通滤波器去除与呼吸、心跳等有关的生理噪声。

6)去线性

由于机器的长期工作而升温或者被试的不适应，随着时间的积累会存在一个线性趋势。线性趋势也属于噪声引起的，因此需进行去除。

7)图像分割

为了去除脑脊液(Cerebro-Spinal Fluid,CSF)、白质(White Matter,WM)等冗余信号，需要对结构像进行分割，利用分割得到的信息制作脑脊液与白质模板；

8)冗余去除

去除白质、脑脊液、全脑信号以及头动伪迹等冗余信号。

2、采用种子点方法在标准MNI空间中提取被试脑区时间序列，包含DMN(23个子区)、ECN(12个子区)、SN(13个子区)的时 间序列。为了准确定位各脑区，首先将静息态磁共振图像配准至MNI标准空间；为提高配准精度，采用两步配准：首先将功能像配准到结构像，其次将结构像配准到标准空间，从而利用所得到的变换矩阵将功能像配准到标准空间并将体素重采样为3mm×3mm×3mm；配准模板采用MNI模板，该模板是由加拿大蒙特利尔神经科学研究所(Montreal Neurological Institute,MNI)研发而成；整个配准过程使用FMRIB提供的线性配准工具实现，配准后图像如图3所示。基于fMRI的研究的种子点区域选取以各子脑区的坐标为圆心，半径为6mm的球形区域，种子点区域包含33个体素，各脑区种子点坐标如图10所示。选取种子点之后，分别对每一个被试提取所有种子点区域内全部体素的一维信号，并对其求加权平均作为种子点的时间序列曲线；整个基于种子点的时间序列提取过程基于FSL软件实现。

3、接下来基于PCA方法进行动态功能连接模式分解。PCA的一般过程是：

1)去除平均值：means；

2)计算A的协方差矩阵Σ；

3)计算Σ的特征向量和特征值；

4)将特征值从大到小排序；

5)保留最上面的k个特征向量(这k个特征向量保证了数据映射到特征值最大的特征向量的方向时，数据间的累积方差最大，数据映射到第二大的特征向量时，数据间的累积方差次之)。

整个PCA动态功能连接模式分解流程如图4所示。首先采用滑动 窗方法计算各脑区时间序列两两之间的相关系数，进而构造基于PCA分析的功能连接矩阵。本文采用的滑动窗窗长为90s，滑动步长为1TR＝645ms，相关系数采用皮尔森相关系数。相关系数ρ考察两个变量的相关程度，取值范围在-1～1之间，其中，1表示变量完全正相关，0表示无关，-1表示完全负相关。

${\mathrm{\ρ}}_{X,Y}=\frac{\mathrm{cov}(X,Y)}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}=\frac{E\[(X-{\mathrm{\μ}}_X)(Y-{\mathrm{\μ}}_Y)\]}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}$

式中：ρ(X,Y)——变量X和Y的相关系数；cov(X,Y)——变量X和Y的协方差；σ_X——变量X的标准方差；σ_Y——变量Y的标准方差；μ_X——变量X的均值；μ_Y——变量Y的均值。

功能连接矩阵的具体构造过程如下，首先计算每个被试DMN、ECN和SN网络全部48个子区两两之间的动态功能连接，得到W个48×48的功能连接矩阵，其中W为滑动窗的个数。其次，取每个功能连接矩阵的上三角进行展开，并将展开后的矩阵以时间顺序进行连接，得到每个被试的功能连接矩阵C_S(S＝1,2,...,s)，为减少被试个体之间的差异，对每个被试的C_S减去C_S每行的平均值，即得到最后，将各个被试的去均值功能连接矩阵进行连接，得到用于PCA分析的功能连接矩阵

得到功能连接矩阵后，对功能连接矩阵进行PCA分析，得到分解后的连接模式，这些连接模式尽可能多地功能连接矩阵的信息，为基本连接模式。首先计算矩阵A的协方差矩阵AA^T的特征值和特征向量：AA^T＝U∧U^T，其中U的各列向量为正交特征向量，∧的对角元素为对 应于特征向量的特征值。这些特征向量即为矩阵的主成分，表示功能连接矩阵的主要特征，为基本连接模式，因此称特征向量为“特征连接”。

正交矩阵U的大小为(N²-N)/2×(N²-N)/2，其中N表示脑区的数量，即得到(N²-N)/2个特征连接。然而，对应的特征值大的特征向量表示的动态功能连接矩阵的累积方差最大，因此，少数的特征连接即可有效地代表矩阵的特征。将特征值按从大到小排列，取前10个特征值对应的特征向量来表示功能连接模式分解结果。

整个基于PCA的动态功能连接模式分解流程使用MATLAB编写程序实现。对DMN、ECN以及SN网络的PCA模式分解结果如图5，对三个网络总体进行模式分解的结果如图6。

4、为研究某些连接模式可能是准稳定的，即它们随时间重复出现以及存在于众多被试中，基于K均值聚类方法对相应的功能连接矩阵进行聚类分析。

整个K均值聚类的动态功能连接模式分解流程如图7所示。首先采用滑动窗方法计算各脑区时间序列两两之间的相关系数，进而构造基于K均值聚类分析的功能连接矩阵。功能连接矩阵的具体构造过程如下，首先计算每个被试DMN、ECN和SN网络全部48个子区两两之间的动态功能连接，得到W个48×48的功能连接矩阵，其中W为滑动窗的个数。其次，取每个功能连接矩阵的上三角进行展开，并将展开后的矩阵以窗的时间顺序进行连接。最后，将各个被试的功能连接矩阵进行连接，得到用于K均值聚类分析的功能连接矩阵B，B的 大小为P×Q，其中P为窗的个数与被试个数的乘积，Q为功能连接对的个数。

得到功能连接矩阵后，采用K均值聚类方法对功能连接矩阵进行行聚类分析，具体为：

1)从n个数据对象中任意选择k个对象作为初始聚类中心M_i(i＝1,2,...,k)。

2)根据每个聚类对象的初值(中心对象)，计算每个对象X_k(k＝1,2,...,n)与这些中心对象的距离D_k＝|M_i-X_k|，并根据最小距离重新对相应的对象进行划分，形成类簇C_i(i＝1,2,...N_i)，N_i表示类簇C_i中包含的对象个数。

3)计算每个聚类的均值作为更新后的聚类中心：

$M_i=\underset{X_k\∈C_i}{\Σ}{X}_k/N_i$

4)重复步骤(2)和步骤(3)，直到每个聚类不再发生变化为止。

整个基于K均值聚类法的动态功能连接模式分解流程使用MATLAB编写程序实现。每个聚类中心代表存在的一种连接模式，即一种“状态”，本文选取的聚类个数为7，对DMN、ECN以及SN三个网络总体进行K均值聚类模式分解的结果如图8。

接下来根据聚类结果可以得到每个被试以时间为函数的状态分配图。图9所示为其中3个被试的状态分配图，可以看出，功能连接在长时间内趋于单一状态，尽管在状态转换过渡时的短时间内也通常会出现其它状态。

除了状态转移图，也可采用马尔可夫链描述其转换行为，将状态转换过程看作马尔科夫链。马尔可夫链具体为：设{X(t),t∈T}是随机过程，其参数集T＝{0,1,2,...}，状态空间S为可数集，设S＝{1,2,...}，若对任意的m≥1及任意的i₀,i₁,...,i_m,j∈S，当时，有

P{X(m+1)＝j|X(0)＝i₀,X(1)＝i₁,...,X(m)＝i_m}

＝P{X(m+1)＝j|X(m)＝i_m}

则称{X(t),t∈T}为离散时间马尔可夫链。通过马尔可夫链得到每个被试的转移矩阵，转移矩阵中的元素表示由一种状态转移到另一种状态的概率。然后将所有被试的转移矩阵进行平均，得到平均转移矩阵，其矩阵的每一行元素乘以被试个数然后相加，结果为整数。

综上所述，本发明目的在于提出基于PCA与K均值聚类的动态功能连接模式方法，不仅能够捕捉不同基本连接模式之间的转变，并且为临床上精神分裂症、阿尔兹海默症、抑郁症等神经精神疾病的进一步研究和防治提供策略。

Claims (3)

1.基于PCA与K均值聚类的动态功能连接模式分解方法，其特征在于，具体步骤如下：

(1)、对被试进行静息态磁共振数据采集，并进行预处理，磁共振数据的预处理目的是提高脑功能图像的信噪比，并将被试图像与标准模板进行配准变换；

(2)、预处理之后，提取DMN、ECN以及SN网络的时间序列，采用种子点法提取DMN23个子区的时间序列、ECN12个子区的时间序列、SN13个子区的时间序列，种子点区域选取半径为6mm的球形区域，选取种子点之后，并对其中所有的体素信号加权平均进而得到种子点的时间序列曲线；

(3)、基于PCA进行动态功能连接模式分解，首先通过采用滑动窗方法计算每个被试所有脑区两两之间的皮尔森相关系数，进而得到功能连接矩阵，功能连接矩阵的大小为M*N，其中M为功能连接对的个数，N为滑动窗的个数与被试个数的乘积；其次，采用PCA方法对功能连接矩阵进行模式分解，得到矩阵特征向量，按从大到小顺序保留最上面的k个特征向量，k个特征向量表示动态功能连接的k个主要连接模式；

(4)、为研究某些连接模式可能是准稳定的，即它们随时间重复出现以及存在于众多被试中，基于K均值聚类方法对相应的功能连接矩阵进行聚类，每个聚类中心代表存在的一种连接模式，实现基于K均值聚类的动态功能连接模式分解。

2.根据权利要求1所述的基于PCA与K均值聚类的动态功能连接模式分解，其特征在于，步骤(3)具体为：

整个PCA动态功能连接模式分解首先采用滑动窗方法计算各脑区时间序列两两之间的相关系数，进而构造基于PCA分析的功能连接矩阵，采用的滑动窗窗长为90s，滑动步长为1TR＝645ms，相关系数采用皮尔森相关系数，相关系数ρ考察两个变量的相关程度，取值范围在-1～1之间，其中，1表示变量完全正相关，0表示无关，-1表示完全负相关：

${\mathrm{\ρ}}_{X,Y}=\frac{\mathrm{cov}(X,Y)}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}=\frac{E\[(X-{\mathrm{\μ}}_X)(Y-{\mathrm{\μ}}_Y)\]}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}$

式中：ρ(X,Y)——变量X和Y的相关系数；cov(X,Y)——变量X和Y的协方差；σ_X——变量X的标准方差；σ_Y——变量Y的标准方差；μ_X——变量X的均值；μ_Y——变量Y的均值；

功能连接矩阵的具体构造过程如下，首先计算每个被试DMN、ECN和SN网络全部48个子区两两之间的动态功能连接，得到W个48×48的功能连接矩阵，其中W为滑动窗的个数；其次，取每个功能连接矩阵的上三角进行展开，并将展开后的矩阵以时间顺序进行连接，得到每个被试的功能连接矩阵C_S(S＝1,2,...,s)，为减少被试个体之间的差异，对每个被试的C_S减去C_S每行的平均值，即得到最后，将各个被试的去均值功能连接矩阵进行连接，得到用于PCA分析的功能连接矩阵

得到功能连接矩阵后，对功能连接矩阵进行PCA分析，得到分解后的连接模式，这些连接模式尽可能多地功能连接矩阵的信息，为基本连接模式，首先计算矩阵A的协方差矩阵AA^T的特征值和特征向量：AA^T＝U∧U^T，其中U的各列向量为正交特征向量，∧的对角元素为对应于特征向量的特征值。这些特征向量即为矩阵的主成分，表示功能连接矩阵的主要特征，为基本连接模式，因此称特征向量为“特征连接”；

正交矩阵U的大小为(N²-N)/2×(N²-N)/2，其中N表示脑区的数量，即得到(N²-N)/2个特征连接，然而，对应的特征值大的特征向量表示的动态功能连接矩阵的累积方差最大，因此，少数的特征连接即可有效地代表矩阵的特征，将特征值按从大到小排列，取前10个特征值对应的特征向量来表示功能连接模式分解结果。

3.根据权利要求1所述的基于PCA与K均值聚类的动态功能连接模式分解，其特征在于，步骤(4)具体为：

首先采用滑动窗方法计算各脑区时间序列两两之间的相关系数，进而构造基于K均值聚类分析的功能连接矩阵，功能连接矩阵的具体构造过程如下，首先计算每个被试DMN、ECN和SN网络全部48个子区两两之间的动态功能连接，得到W个48×48的功能连接矩阵，其中W为滑动窗的个数；其次，取每个功能连接矩阵的上三角进行展开，并将展开后的矩阵以窗的时间顺序进行连接；最后，将各个被试的功能连接矩阵进行连接，得到用于K均值聚类分析的功能连接矩阵B，B的大小为P×Q，其中P为窗的个数与被试个数的乘积，Q为功能连接对的个数；

得到功能连接矩阵后，采用K均值聚类方法对功能连接矩阵进行行聚类分析，具体为：

1)从n个数据对象中任意选择k个对象作为初始聚类中心M_i(i＝1,2,...,k)；

2)根据每个聚类对象的初值，计算每个对象X_k(k＝1,2,...,n)与这些中心对象的距离D_k＝|M_i-X_k|，并根据最小距离重新对相应的对象进行划分，形成类簇C_i(i＝1,2,...N_i)，N_i表示类簇C_i中包含的对象个数；

3)计算每个聚类的均值作为更新后的聚类中心：

$M_i=\underset{X_k\∈C_i}{\Σ}{X}_k/N_i$

4)重复步骤(2)和步骤(3)，直到每个聚类不再发生变化为止。