基于静息态功能连接的杏仁核谱聚类分割方法

摘要

一种基于静息态功能成像的杏仁核谱聚类分割方法，是一种基于谱聚类算法根据杏仁核内部体素功能的相似性对该脑区进行自动化高效率的方法，先对静息态磁共振数据预处理、再杏仁核脑区提取、然后杏仁核内部体素全脑功能连接计算以及最后对功能连接矩阵谱聚类分割，本发明提出的自动分割算法与杏仁核临床解剖结果取得很大程度的一致性，并且在稳定性和抗噪声干扰方面也取得较满意的结果，相对于传统手动分割方法，更加简单，方便，高效、可重复性高等多方面的优点。

说明书

技术领域

本发明属于图像处理领域，具体涉及一种基于静息态功能连接的杏仁核谱聚类分割方法，尤其涉及使用谱聚类算法基于杏仁核各体素的全脑连接模式相似性进行杏仁核子脑区自动化分割的方法。

背景技术

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

杏仁核是人类情绪记忆网络中最为核心的脑区，它在情绪的产生和表达、社会判断和面孔识别等认知活动中扮演了重要的角色。杏仁核结构和功能的异常会引起孤独症、抑郁症、帕金森症等精神疾病，它已经成为神经病理学研究的热点，但目前国内外，杏仁核内部结构和功能的细化研究相对较少，如何快速精确自动化对杏仁核脑区进行划分已经显得尤为重要。人脑是中枢神经系统最高级的部分，它是由140亿个神经细胞构成的具有高度组织性的器官，且具有相似功能的细胞团体活动具有显著的一致性，同样，在杏仁核内部，具有相同功能的细胞核团其功能也会体现出某种程度上的一致性。研究表明静息 态时大脑不同脑区内的低频波动BOLD信号(0.01～0.1Hz)是相关的，运动区、视觉、语言区、听觉区的低频振幅波动的相关性相继被证实，这些相关的低频波动构成了静息态功能连接网络，且不同人群间的静息态功能连接网络是一致的。在杏仁核中，具有相同功能的细胞核团，其静息态功能也具有高度的相似性。基于这样的理论基础本发明通过静息态功能成像根据特定脑区内各体素全脑的功能连接模式的相似性以谱聚类对杏仁核脑区进行分割，尝试深入研究杏仁核内部的功能，全面了解杏仁核的作用机制来为临床上相关疾病的防治提供策略。

目前，杏仁核进行分割主要通过有一定医学解剖基础的临床医生根据杏仁核脑区细胞核团结构的特征手动地对杏仁核脑区进行划分。这样的分割方法一方面对研究者有很高的医学解剖背景要求，另一方面效率低下，耗时长，不利于推广。

发明内容

基于依据杏仁核细胞核团结构特征来手动进行脑区划分的效率低下和耗时长，受主观因素影响较大的现状，本发明的目的在于提出一种基于静息态功能连接的杏仁核谱聚类分割方法，能够稳定、高效、准确地对人脑杏仁核区域进行自动化分割。

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

基于静息态功能连接的杏仁核谱聚类分割方法，具体步骤如下：

(1)、对被试进行静息态磁共振数据采集，并将采集的数据进行预处理，预处理是指提高脑功能图像的信噪比，并在保留图像数据细 节的同时，使被试图像与标准模板进行仿射配准变换；

(2)、提取被试杏仁核脑区，将预处理过的静息态功能磁共振图像配准至MNI标准空间，采用AAL(automated anatomical labeling)模板将皮层和皮层下的灰质划分成90个脑区，左右杏仁核分别对应为41号脑区和42号脑区,同样将杏仁核ROI配准到MNI空间，从而提取出被试杏仁核脑区；

(3)、通过计算杏仁核中每个体素的全脑连接模式，进而得到每两个体素之间的全脑连接模式相似度矩阵，提取杏仁核每个体素的时间序列信号；对其他88个脑区的静息态信号采用脑区内所有体素信号平均的方式提取；其次，通过计算每个杏仁核内部的体素信号与其他88个脑区的信号之间的相关系数，得到杏仁核体素的全脑连接矩阵M，其中元素(i,j)表示第i个杏仁核体素到第j个脑区的功能连接强度，该矩阵的每一行描述了一个杏仁核体素在静息态下和全脑其他脑区之间的连接模式，值越接近1，相关性越大；然后计算每两行连接模式之间的相关系数就能够度量两个体素在全脑功能连接模式上的相似度，记为相似度矩阵N，其中每个元素表示杏仁核第i个体素到第j个体素在全脑功能连接模式上的相似度；

(4)、运用谱聚类对步骤3中得到的全脑连接模式相似度矩阵N进行聚类分割，实现杏仁核的子区域分割。

步骤(4)具体为：

(1)令样本点x₁,x₂,…,x_n表示任意的被聚类的N每行的功能连接向量，令s_ij表示构造图两个顶点间的连接权值，以高斯相似性函数 来定义顶点之间的权值，其中参数σ称为尺度参数，其中s_ij定义为

s_ij＝exp(-||x_i-x_j||²/2σ²)

则样本的相似性矩阵为S＝(s_ij)(i,j＝1,2,…,n)；

(2)输入：相似性矩阵S∈R^n×n，聚类类别数k；建立样本的相似性连接图，令W为其权值矩阵；

(3)计算非规范化拉普拉斯矩阵L，把W的每一列元素加起来得到N个数，将它们放在对角线上，其他地方都是零，组成一个N*N的矩阵，记为D，并令L＝D-W；计算L的前k个最小的特征值所对应的特征向量v₁,v₂,…v_k；

(4)令V∈R^n×k为v₁,v₂,…v_k由按列排列所组成的矩阵；

(5)对于i＝1,2,…,n，令y_i∈R^k为矩阵V的第i列；

(6)利用k均值算法将属于空间R^k的数据y_i(i＝1,2,…,n)聚成k类，即为C₁,…C_k，便完成了杏仁核各像素的聚类，实现杏仁核的子区域分割。

本发明的创新点在于：与杏仁核解剖图谱对比可以发现本发明提出的基于静息态功能连接的杏仁核谱聚类自动分割算法与杏仁核临床解剖结果取得很大程度的一致性，另一方面定量评价分割结果表明，本发明解剖图谱与生理图谱有70％以上的相似度，并且在稳定性和抗噪声干扰方面也取得较满意的结果，相对于传统手动分割方法，更加简单，方便，高效、可重复性高等多方面的优点。

附图说明

图1是本发明的流程图。

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

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

图4是头动校正结果。

图5是功能像配准前后对比，左侧为配准前功能像，中间为MNI152标准模板，右侧为配准后功能像，经过配准图像与模板之间的差异减小。

图6是杏仁核生理解剖图谱和本发明提出的方法对杏仁核分割结果对比图。

图7是谱聚类分割杏仁核结果与juiche图谱重叠率统计图。

图8是30个被试左侧杏仁核分割结果与脑图谱空间相似性对比。

图9是30个被试右侧杏仁核分割结果与脑图谱空间相似性对比。

图10是不同信噪比下杏仁核分割结果稳定性对比。

具体实施方式

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

本发明基于静息态功能连接进行杏仁核脑区分割方法原理如图1所示。

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

对采集的原始静息态磁共振数据进行预处理,运用fMRI数据预处理，在Linux系统Ubuntu14.04下基于AFNI和FSL软件编写Batch文件(数据批处理)完成。预处理流程见图2，主要包含以下几个方面：

1)结构像头颅分割

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

2)时间对齐

血液动力学函数表明血液对刺激的响应有一定的时间延迟，由于一个TR期间采集到全脑图像，这导致每一层的图像并非在同一时刻采集，而是发生在整个扫描时间段内，时间校正就是通过类似于插值的方法对每一层图像进行层时间处理，使得一个TR周期内各层图像近似在同一时刻获取。

3)头动校正

进行头动校正时一般将被试的大脑看做一个刚体，因此在fMRI实验中被试头部的运动可以近似成一种刚体运动，即只有平移变换与旋转变换的组合。选择单个被试的第一帧图像作为参考图像，通过AFNI的3dvolreg函数使得其余的所有图像与参考图像配准，结果如图4所示。如果头动超过一个体素则去除该被试。

4)空间平滑

空间平滑采用高斯函数进行高斯平滑，能有效地消弱随机噪声对fMRI信号的影响，提高数据的信噪比。三维高斯函数是比较常用的空间平滑方法，其半高全宽决定了空间平滑的力度，本文选用半高全宽(Full Width at Half Maximum,FWHM)为6mm的高斯核函数进行数据平滑。

5)时域带通滤波

静息态fMRI信号是一种低频波动，频率主要集中在0.01～0.1Hz，这种低频波动反映了自发的神经活动。因此采用频率范围为0.01～0.1Hz的带通滤波器去除与呼吸、心跳等有关的生理噪声。

6)去线性

由于机器的长期工作导致温度升高或者被试的不适应，随着时间的积累会存在线性漂移，需要进行去线性。

7)图像分割

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

8)冗余去除

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

(2)、在标准MNI空间中提取被试杏仁核脑区。为了准确定位杏仁核脑区，首先将静息态磁共振图像配准至MNI标准空间；为提高配准精度，采用两步配准：首先将功能像配准到结构像，其次将结构 像配准到标准空间，从而利用所得到的变换矩阵将功能像配准到标准空间并将体素重采样为3mm×3mm×3mm；配准模板采用MNI模板，该模板是由加拿大蒙特利尔神经科学研究所(Montreal Neurological Institute,MNI)研发而成；整个配准过程使用FMRIB提供的线性配准工具实现，配准前后图像对比如图5所示；AAL(automated anatomical labeling)模板将皮层和皮层下的灰质划分成90个脑区，左右杏仁核分别对应为41号脑区和42号脑区,将41号脑区和42号脑区同样配准至标准MNI空间，作为模板便可以作为我们进一步进行分割的目标脑区。

3、由步骤2在已经获得静息态功能磁共振数据杏仁核的位置基础上，可以获得杏仁核单个体素的时间序列，另一方面，计算AAL模板中除了左右侧杏仁核以外的其它88个脑区平均时间序列，再求得杏仁核每个体素与这88各脑区平均时间序列的相关系数，得到杏仁核体素的全脑连接矩阵M,其中元素(i,j)表示第i个杏仁核体素到第j个脑区的功能连接强度；该矩阵的每一行描述了一个杏仁核体素在静息态下和全脑其他脑区之间的连接模式，值越接近1，相关性越大；整个静息态磁共振图像数据的读取、时间序列的提取和相关性计算均基于python 2.7语言实现，本算法科学计算部分主要使用Numpy(Numeric Python)模块实现，NumPy是一个基础科学的计算包，它提供了许多高级的数值编程工具，如：矩阵数据类型、矢量处理，以及精密的运算库，专为进行严格的数字处理而产生；对磁共振数据的读取和写入则依赖NiBabel模块，NiBabel提供了对常见医学影像数据 格式读取的接口。

4、运用谱聚类对步骤3中得到的全脑连接模式矩阵N进行聚类分割，具体为：

(1)令样本点x₁,x₂,…,x_n表示任意的被聚类的N每行的功能连接向量，令s_ij表示构造图两个顶点间的连接权值，以高斯相似性函数来定义顶点之间的权值，其中参数σ称为尺度参数，其中s_ij定义为

s_ij＝exp(-||x_i-x_j||²/2σ²)

则样本的相似性矩阵为S＝(s_ij)(i,j＝1,2,…,n)；

(2)输入：相似性矩阵S∈R^n×n，聚类类别数k；建立样本的相似性连接图，令W为其权值矩阵；

(3)计算非规范化拉普拉斯矩阵L，把W的每一列元素加起来得到N个数，将它们放在对角线上，其他地方都是零，组成一个N*N的矩阵，记为D，并令L＝D-W；计算L的前k个最小的特征值所对应的特征向量v₁,v₂,…v_k；

(4)令V∈R^n×k为v₁,v₂,…v_k由按列排列所组成的矩阵；

(5)对于i＝1,2,…,n，令y_i∈R^k为矩阵V的第i列；

(6)利用k均值算法将属于空间R^k的数据y_i(i＝1,2,…,n)聚成k类，即为C₁,…C_k，便完成了杏仁核各像素的聚类，实现杏仁核的子区域分割。

下面就本发明提出的基于静息态功能连接的杏仁核分割方法取得的实验结果与杏仁核临床生理解剖图谱做比较，并且对本方法的稳定性、抗噪声干扰等方面做一定的探讨。

根据Amunts等人在05年基于细胞构筑的研究，julich图谱计算并构建了这些杏仁核亚区的概率模板。该图谱将杏仁核分成了三个子区域，分别为基底外侧核(the laterobasal,LB)、浅表核(the superficial,SF)和中央内侧核(the centromedial,CM)，分别对应图6左侧中蓝色、绿色和红色脑区。LB包括副层状核、底内侧核、基底外侧核和外侧核，CM由内侧核和中央核组成，SF包括杏仁核梨状皮质移行区、杏仁前区、杏仁核-海马区以及腹侧和后侧皮质核。图6中julich图谱显示出杏仁核子区域的空间位置关系，图6右侧则为本发明基于杏仁核内部体素与其它各脑区功能连接进行谱聚类对杏仁核划分的结果，二者进行对比可发现本发明的聚类结果和图谱具有高度一致性，表明这一方法的可行性。

另一方面，为了进一步定量评价分割结果，本发明重复对30个被试的静息态功能图像数据进行杏仁核脑区分割，从算法的精确性和稳定性两个角度对其进行探讨。

卡尔·皮尔逊基于弗朗西斯·高尔顿的算法提出了皮尔逊线性相关系数，通常该相关系数也被称作“皮尔森相关系数”。相关系数ρ考察两个变量的相关程度，取值范围在-1～1之间，其中，1表示变量完全正相关，0表示无关，-1表示完全负相关。

${\rho }_{X,Y}=\frac{\mathrm{cov}(X,Y)}{{\sigma }_X{\sigma }_Y}=\frac{E\left[(X-{\mu }_X)(Y-{\mu }_Y)\right]}{{\sigma }_X{\sigma }_Y}$

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

我们将大小m×n×o的待测三维图像和参考图像B看作为随机变量，则两个图像之间的相关系数为

$CC(A,B)=\frac{\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}\underset{k=1}{\overset{o}{\Sigma }}(A_{i,j,k}-\bar{A})(B_{i,j,k}-\bar{B})}{\sqrt{\left(\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}\underset{k=1}{\overset{o}{\Sigma }}{\left(A_{i,j,k}-\bar{A}\right)}^2\right)\left(\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}\underset{k=1}{\overset{o}{\Sigma }}{\left(B_{i,j,k}-\bar{B}\right)}^2\right)}}$

式中A表示待评价三维图像，B表示参考图像，A_i,j,k表示待测图像中的像素点，B_i,j,k表示参考图像中的像素点，表示图像A的均值，表示图像B的均值，m表示图像行数，n表示图像的列数，o表示图像高度。图像A和图像B的相关系数越小(越接近0)，说明图像A和图像B的相似度越小；相关系数越大(越接近1)，说明图像A和B的相似度越大。计算本发明的聚类组分析结果分别与julich图谱对应的杏仁核子区域的皮尔森相关系数如图7所示，结果越接近1表示二者空间分布越相似，由图可知右侧CM亚区分割结果与图谱最接近(CC＝0.78)，右侧SF亚区相对谱图相似度较低(CC＝0.65)。

本发明将每个被试杏仁核亚区的分割结果与该数据集的组分析结果进行空间相关，以各子区域与组分析结果的空间重叠率的变化剧烈程度来反映本文分割算法的稳定性，图8和图9描述了左侧和右侧杏仁核子区单个被试与组分析结果的重叠情况，通过计算右侧CM空间相关系数的方差为0.0303最大，表明该亚区分割的结果稳定性比较弱，而右侧LB脑区空间相关系数的方差为0.0017最小，所以这个脑区聚类算法的分割结果可重复性最高，稳定性也最好(其它脑区的方差分别为var(CML)＝0.0041,var(LBL)＝0.0185,var(SFR)＝ 0.0098,var(SFL)＝0.0194)

在讨论了分割算法的精确性和稳定性的基础上，叠加不同强度的噪声到杏仁核体素的时间序列中再进行脑区划分来考察本发明算法的健壮性。首先将杏仁核左右脑区各体素的时间序列提取出来，通过计算将信噪比分别为10db,30db,50db,70db,90db的高斯白噪声叠加到各时间序列中再进行图1所示的基于谱聚类基础上的杏仁核亚区分割算法，整个实验模拟过程在python 2.7平台实现。

叠加噪声后杏仁核6个亚区的聚类分割如图10所示，其中(A)图反映了不同噪声下基于ROI体素时间序列与ALL模板其余88个脑区平均序列计算的功能连接矩阵进行分割结果与没有加入噪声的情况下各子区空间相似性变化情况，由图可知不同信噪比的图像数据均取得了比较满意的分割结果，在信噪比大于70db的情况下，各子区都保持了近乎与没有加入噪声的结果100％的相似性，在10db的低信噪比下各区皮尔森相关系数也能保持0.7左右。

综上所述，本发明目的在于提出一种基于静息态功能磁共振成像的能够稳定、高效、准确的自动化人脑杏仁核区域分割方法，为临床上抑郁症、孤独症等情绪相关精神疾病的进一步研究和防治提供策略。