CN103020653A - 基于网络分析的结构和功能磁共振图像联合分类方法 - Google Patents

Abstract

提供了一种基于网络分析的结构和功能磁共振图像联合分类方法。首先建立结构和功能的脑网络模型，计算脑网络的特征路径长度、集群度和网络中心性，用以表征不同的图像模式；然后利用这些网络参数来训练一个自适应提高分类器。本发明能够利用磁共振图像中尽可能多的信息，脑网络参数能够从本质上反应脑的活动，同时采用了多分类器技术，弥补了传统分类方法不能体现脑活动固有属性的不足，能够精确的对脑磁共振图像进行分类。

Description

基于网络分析的结构和功能磁共振图像联合分类方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于网络分析的结构和功能磁共振图像联合分类方法。

背景技术

磁共振成像(Magnetic Resonance Imaging，MRI)以其高时空分辨率，非侵入式等特点在神经疾病诊断治疗方面得到了广泛应用。脑是一个复杂的系统，在受到刺激条件或经历病变时脑的磁共振图像会发生相应的变化。利用图像分类方法，计算磁共振图像具有某种属性的可能性大小，或者自动判别图像的类别属性，是计算机辅助分析的一个重要应用。

传统的磁共振图像分类方法主要有感兴趣区域(ROI)方式和体素(voxel)方式两种分类方法。感兴趣区域方式的分类方法依据目标结构的先验知识，将样本和目标分割成多个目标区域，并据此对目标进行分类；体素方式的分类方法采用复杂的非线性配准，以最大限度地实现个体间的精确对应，然后以图像的每一个空间单位(体素)作为分类依据。这两种方法都假设目标与样本的内部组织结构是一一对应的。前者认为先验的图像区域存在于每一个目标图像当中，并且能够准确分割；后者假定非线性配准后的体素是一一对应的。然而，这样的假设在很多情况下并不合理。人在不同状态下的脑磁共振图像会受到多方面因素的干扰，传统的分类方法都不是根据脑的固有属性对脑磁共振图像进行分类的，因此都会导致分类性能的下降。

发明内容

为了克服已有技术的不足，本发明提供了一种分类准确率高、泛化性能强的脑磁共振图像分类方法。

为实现上述目的，本发明提供了一种基于网络分析的结构和功能磁共振图像联合分类方法，包括：

Sa：对训练样本图像和测试样本图像进行预处理，然后按照标准脑结构模板将全脑分为90个脑区，并提取功能图像各个脑区的平均时间序列和结构图像各个脑区的灰质密度；

Sb：计算各个平均时间序列之间以及各个脑区的灰质密度之间的偏相关系数，分别得到偏相关系数矩阵；

Sc：将所述偏相关系数矩阵二值化，得到复杂网络模型；

Sd：分别计算功能网络和结构网络中各个节点的特征路径长度、集群度和网络中心性；

Se：利用训练样本功能图像和结构图像的特征路径长度、集群度和网络中心性分别作为训练样本图像的特征训练自适应提高分类器，得到多个自适应提高分类器以及对应的分类正确率；

步骤Sf：利用训练好的多个自适应提高分类器以投票的方式进行组合；

Sg：利用训练好的自适应提高分类器对测试样本图像进行分类。

本发明能够利用磁共振图像中尽可能多的信息，脑网络参数能够从本质上反应脑的活动，同时采用了多分类器技术，弥补了传统分类方法不能体现脑活动固有属性的不足，能够精确的对脑磁共振图像进行分类。

附图说明

图1是本发明所述算法流程图；

图2是测试数据上两种方法分类受试者操作特性(ROC)曲线对比。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

基于网络分析的结构和功能磁共振图像联合分类方法是一种全新的磁共振图像分类方法。该方法首先建立结构和功能的脑网络模型，计算脑网络的特征路径长度、集群度和网络中心性，用以表征不同的图像模式；然后利用这些网络参数来训练一个自适应提高(adaboost)分类器。

参照图1，根据本发明所述的一种人脑磁共振成像图像分类方法，能够依据训练样本图像来确定测试样本图像的类别，具体实施步骤如下：

步骤Sa，对训练样本图像和测试样本图像进行预处理，然后按照标准脑结构模板将全脑分为90个脑区，并提取功能图像各个脑区的平均时间序列和结构图像各个脑区的灰质密度；

1.脑磁共振图像的预处理

由于磁共振扫描过程中各种各样的噪声的影响，被试个体自身存在尺度和位置上的差异，非常有必要在分析数据之前对数据做一定的预处理。在整个的实验的数据获取中，主要的噪声信息来源有：(1)物理头动；(2)图像内层间扫描时间差别；(3)外在磁场的不均匀性等。脑磁共振图像预处理的常见步骤有：图像序列对齐，联合配准，标准化(或称均一化)，空间平滑滤波等。

2.脑磁共振图像的分割

采用国际通用的结构标记模板(AAL)，将全脑分为90个脑区。结构标记模板是磁共振图像研究领域使用最为广泛的脑结构模板。

3.提取各脑区的平均时间序列

依据预处理后的脑功能磁共振图像的数据，提取包含于相应脑区内部的各个体素在不同时间点上激活值的时间序列Y(矩阵维数D×N)，其中D为包含于球体内部的体素数目，N为时间点数。所述激活值是指各个体素在不同时间点上的血氧水平依赖(BOLD)强度。然后，将各个体素的激活值进行平均，得到脑区平均时间序列。

4.提取各脑区的灰质密度

依据预处理后的脑结构磁共振图像的数据，制作半径为4mm的球形感兴趣区(ROI)，提取包含于感兴趣区内的33个2mm×2mm×2mm的体素的灰质密度Y，其中，33个2mm×2mm×2mm的体素的灰质密度Y距离每个脑区内部的中心点最近。

步骤Sb：计算各个平均时间序列之间以及各脑区灰质密度之间的偏相关系数。该步骤Sb具体包括如下步骤：

1.计算协方差系数

依据步骤Sa提取的各个脑区的时间序列或灰质密度，计算各个平均时间序列之间以及各脑区灰质密度之间的协方差矩阵S，S的每个元素s_i，j为第i个和第j个时间序列之间或灰质密度之间的协方差系数，

$s_{i,j}=\frac{1}{M}\underset{t=1}{\overset{M}{\mathrm{\Σ}}}(x_i\left(t\right)-\stackrel{\‾}{x_i})(x_j\left(t\right)-\stackrel{\‾}{x_j})$

其中，对于功能图像，M为时间点数目，x_i(t)(i＝1，…，M)为第i个时间序列，

为第i个时间序列的平均值，

为第j个时间序列的平均值；对于结构图像，M为脑区内体素数目，x_i(t)(i＝1，…，M)为第i个脑区的灰质密度，为第i个脑区灰质密度的平均值，

为第j个脑区灰质密度的平均值。

2.计算偏相关系数

依据协方差系数矩阵S(矩阵维度为90×90)，计算偏相关系数矩阵R(矩阵维度为90×90)，R的每个元素r_i，j为：

$r_{i,j}=-\frac{s_{i,j}^{-1}}{\sqrt{s_{i,i}^{-1}{s}_{j,j}^{-1}}}$

其中，

为协方差矩阵S(矩阵维度为90×90)的逆矩阵的第{i，j}个元素。

3.对偏相关系数进行Fisher变换

依据偏相关系数矩阵R(矩阵维度为90×90)，计算经过Fisher变换的偏相关系数矩阵F(矩阵维度为90×90)，F的每个元素f_ij为：

$f_{i,j}=\frac{1}{2}\left(\frac{1+r_{i,j}}{1-r_{i,j}}\right),$

其中，f_ij为经过Fisher变换后的偏相关系数矩阵F(矩阵维度为90×90)的第{i，j}个元素，r_ij为偏相关系数矩阵R(矩阵维度为90×90)的第{i，j}个元素。

步骤Sc：将偏相关系数矩阵二值化，得到复杂网络模型；

设定阈值T’，令经过Fisher变换后的偏相关系数矩阵F(矩阵维度为90×90)中大于等于T’的值为1，小于T’的值为0，得到复杂网络模型。二值化后的矩阵中1表示两个脑区之间有连接，即网络中两个节点之间的边存在，0则表示两个脑区之间没有连接，即网络中的两个节点之间没有边。阈值选取的方法为：使网络中实际存在的边的数量是网络中可能存在的边的数量

其中N为网络中节点的数目)的十分之一。二值化的过程可描述为令

$w_{i,j}=\left\{\begin{array}{c}1,\left|f_{i,j}\right|\&GreaterEqual;T^{\&prime;}\\ 0,\left|f_{i,j}\right|<T^{\&prime;}\end{array}\right.,$

其中，w_ij为二值化后的网络的第{i，j}个元素，f_ij为经过Fisher变换的偏相关系数矩阵F(矩阵维度为90×90)的第{i，j}个元素，T’为选取的阈值，|·|为绝对值计算符号。

步骤Sd：计算网络中各个节点的特征路径长度、集群度和网络中心性；

依据复杂网络模型，计算网络的特征路径长度、集群度和网络中心性，作为磁共振图像的特征。

特征路径长度提供了网络中某一节点的信息到达另一节点的最优路径。我们可以用特征路径长度矩阵描述网络中任意两个节点i，j的特征路径长度lij。网络平均特征路径长度L描述了网络中任意两个节点的特征路径长度的平均值，即

$L=\frac{1}{N(N-1)}\underset{i,j\&Element;V,i\&NotEqual;j}{\mathrm{\&Sigma;}}{l}_{\mathrm{ij}}$

其中，N为网络中节点的个数，即分割的脑区数90；l_ij为节点i，j之间的特征路径长度，V为网络中所有节点的集合。

集群度是度量网络性质的另一个重要特征，用来量度某一节点的相邻节点互为邻居的可能性。某一节点i的集群度C_i的值等于它的相邻节点之间存在的边的数目与它们之间所有可能的边数的比值，即

$C_i=\frac{e_i}{\frac{k_i(k_i-1)}{2}}=\frac{{2e}_i}{k_i(k_i-1)}$

其中，e_i表示节点i的邻点之间存在的边数，k_i表示节点i的邻点的数目，

就表示节点i的邻点之间可能存在的边数。

网络中各节点的网络中心性描述了网络中各个节点的重要程度；节点i的网络中心性定义为网络中除i以外的任意两个节点之间通过了节点i的特征路径的数目与该两个节点之间所有特征路径数目的比值，即

$B_i=\underset{i\&NotEqual;j\&NotEqual;k}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&delta;}}_{\mathrm{jk}}\left(i\right)}{{\mathrm{\&delta;}}_{\mathrm{jk}}}$

其中，B_i为节点i的网络中心性，δ_jk为节点j和k之间的最短路径的数目，δ_jk(i)为节点j和k之间的最短路径中经过了节点i的数目。

步骤Se：利用训练样本功能图像和结构图像的特征路径长度、集群度和网络中心性分别作为训练样本图像的特征各训练一自适应提高(adaboost)分类器，得到多个自适应提高(adaboost)分类器以及对应的分类正确率，自适应提高(adaboost)分类器可以为6个，或根据需要可以增加或减少自适应提高(adaboost)分类器的数量；

得到训练样本图像的特征后，首先将每类特征中每个节点的信息作为一个线性分类器，用线性分类器的加权和组成一个新的自适应提高(adaboost)分类器，最初每个分类器的权重设为

(m为样本图像的数目)，自适应提高(adaboost)分类器在训练过程中逐渐调整每个线性分类器的权重，最后得到一个最优的自适应提高(adaboost)分类器。具体实施步骤如下：

对给定的样本(x₁，y₁)，…，(x_m，y_m)，其中x_i∈X，y_i∈Y＝(-1，1)，X为训练样本图像的网络特征，Y为图像类别，首先设定初始化分类器的权重为

之后进行T次迭代，迭代过程如下：

变量t从1开始增加到T，每次迭代首先计算每个特征h_t对训练样本图像进行分类得到的分类误差ε_t，然后计算新的样本权重，

${\mathrm{\&alpha;}}_t=\frac{1}{2}\ln \left(\frac{{1-\mathrm{\&epsiv;}}_t}{{\mathrm{\&epsiv;}}_t}\right),$

最后，更新各线性分类器的权重，

$D_{t+1}\left(i\right)=\frac{D_t\left(i\right)}{Z_t}\left\{\begin{array}{c}{e}^{-{\mathrm{\&alpha;}}_t},h_t\left(x_i\right)=y_i\\ e^{{\mathrm{\&alpha;}}_t},h_t\left(x_i\right)\&NotEqual;y_i\end{array}\right.,$

其中Z_t为归一化因子。

循环结束后得到最优自适应提高(adaboost)分类器：

$H\left(x\right)=\mathrm{sign}\left(\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t{h}_t\left(x\right)\right).$

步骤Sf：利用训练好的多个自适应提高(adaboost)分类器以投票的方式进行组合得到：

$F\left(x_i\right)=\mathrm{sign}\left(\underset{t=1}{\overset{6}{\mathrm{\&Sigma;}}}(w_t\&CenterDot;F_t\left(x_i^t\right))\right)$

其中，x_i是第i个样本的特征，是第i个样本的第t类特征，w_t是利用第t类特征得到的分类器的分类正确率，F_t是第t类特征的分类判别率，F是组合分类器的输出。

步骤Sg：利用训练好的该自适应提高(adaboost)分类器对测试样本图像进行分类。

将测试样本输入上述步骤得到的最优自适应提高(adaboost)分类器，对测试样本图像进行分类，分类结果通过分类正确率、真阳性率和假阳性率输出。

本发明所述的基于网络分析的结构和功能磁共振图像联合分类方法的效果，可通过真实的磁共振成像数据得以说明：

(1)真实数据实验过程

为展示本发明的效果，在实施方案中采用真实数据集作测试，共41个被试参与了实验，21男、20女。被试年龄段及临床痴呆分级信息见表格1。实验采用T2*加权梯度回波平面成像(Echo-Planar Imaging，EPI)序列获取针刺刺激后BOLD fMRI静息数据。

采用统计参数图(SPM)软件(http://www.fil.ion.ucl.ac.uk/spm/)对数据进行预处理，包括图像序列对齐，联合配准，标准化(或称均一化)、空间平滑滤波。使用本发明所述方法(方法A)对比现有的基于局部特征的分类方法(方法B)，获取分类方法的受试者操作特性(ROC)曲线及其曲线下面积(AUC)，并将ROC曲线和AUC作为分类器性能的度量。

表格1被试信息

(2)实验结果

在真实实验数据集上两种方法的分类ROC曲线分别在图2中显示，其中，图2中的真阳性率是指实际为阳性而按该筛检试验的标准被正确地判为阳性的百分比，假阳性率是指实际为阴性而按该筛检试验的标准被错误地判为阳性的百分比。如图2所示，方法A的ROC曲线在大部分阈值范围内高于方法B；AUC值对比情况：方法A的AUC值为0.9396，方法B的AUC值为0.8578。曲线下面积(AUC)能度量总体分类性能、后验概率和排序性能，AUC值越大，则该分类方法的总体性能越好。由此，方法A效果好于方法B。

实验结果说明，本发明所述的基于网络分析的结构和功能磁共振图像联合分类方法，有效地提高了脑磁共振图像的分类性能。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的权利要求书的保护范围之内。

Claims (10)

1.一种基于网络分析的结构和功能磁共振图像联合分类方法，包括：

Sa：对训练样本图像和测试样本图像进行预处理，然后按照标准脑结构模板将全脑分为90个脑区，并提取功能图像各个脑区的平均时间序列和结构图像各个脑区的灰质密度；

Sb：计算各个平均时间序列之间以及各个脑区的灰质密度之间的偏相关系数，分别得到偏相关系数矩阵；

Sc：将所述偏相关系数矩阵二值化，得到复杂网络模型；

Sd：分别计算功能网络和结构网络中各个节点的特征路径长度、集群度和网络中心性；

Se：利用训练样本功能图像和结构图像的特征路径长度、集群度和网络中心性分别作为训练样本图像的特征(训练)自适应提高分类器，得到多个自适应提高分类器以及对应的分类正确率；

步骤Sf：利用训练好的多个自适应提高分类器以投票的方式进行组合；

Sg：利用训练好的自适应提高分类器对测试样本图像进行分类。

2.根据权利要求1所述的方法，其特征在于所述提取功能图像各个脑区的平均时间序列包括：

分别提取每个脑区内部各个体素在不同时间点上的激活值，再将各个体素的激活值进行平均，得到脑区平均时间序列。

3.根根据权利要求1所述的方法，其特征在于所述提取结构图像各个脑区的灰质密度包括：

提取每个脑区内部的距离中心点最近的33个体素的灰质密度作为脑区的灰质密度。

4.根据权利要求1所述的方法，其特征在于所述计算各个平均时间序列之间以及各个脑区的灰质密度之间的偏相关系数包括：

计算协方差矩阵S，该协方差矩阵维度为90×90，S的每个元素s_i，j为第i个和第j个时间序列或灰质密度之间的协方差系数，

$s_{i,j}=\frac{1}{M}\underset{t=1}{\overset{M}{\mathrm{\&Sigma;}}}(x_i\left(t\right)-\stackrel{\&OverBar;}{x_i})(x_j\left(t\right)-\stackrel{\&OverBar;}{x_j}),$

其中，对于功能图像，M为时间点数目，x_i(t)(i＝1，…，M)为第i个时间序列，

为第i个时间序列的平均值，

为第j个时间序列的平均值；对于结构图像，M为脑区内体素数目，x_i(t)(i＝1，…，M)为第i个脑区的灰质密度，

为第i个脑区灰质密度的平均值，

为第j个脑区灰质密度的平均值。

然后，计算偏相关系数矩阵R，该偏相关系数矩阵R的维度为90×90，R的每个元素r_i，j为：

$r_{i,j}=-\frac{s_{i,j}^{-1}}{\sqrt{s_{i,i}^{-1}{s}_{j,j}^{-1}}}.$

其中，为协方差矩阵S的逆矩阵的第{i，j}个元素。

最后，对偏相关系数进行Fisher变换，得到经过Fisher变换后的偏相关系数矩阵F该变换后的偏相关系数矩阵维度为90×90。

5.根据权利要求4所述的方法，其特征在于：

选取阈值将经过Fisher变换的偏相关系数矩阵F二值化，该变换后的偏相关系数矩阵维度为90×90，二值化后1表示网络中两个节点之间的边存在，0则表示网络中的两个节点之间没有边。

6.根据权利要求5所述的方法，其特征在于所述选取阈值包括：

使网络中实际存在的边的数量是网络中可能存在的边的数量的十分之一。

7.根据权利要求1所述的方法，其特征在于所述计算特征路径长度包括：

用特征路径长度矩阵描述网络中任意两个节点i，j的特征路径长度l_ij，网络平均特征路径长度L描述了网络中任意两个节点的特征路径长度的平均值，即

$L=\frac{1}{N(N-1)}\underset{i,j\&Element;V,i\&NotEqual;j}{\mathrm{\&Sigma;}}{l}_{\mathrm{ij}}$

其中，N为网络中节点的个数，即分割的脑区数90；l_ij为节点i，j之间的特征路径长度，V为网络中所有节点的集合。

8.根据权利要求1所述的方法，其特征在于所述计算集群度的包括：

某一节点i的集群度C_i的值等于它的相邻节点之间存在的边的数目与它们之间所有可能的边数的比值，即

$C_i=\frac{e_i}{\frac{k_i(k_i-1)}{2}}=\frac{2e_i}{k_i(k_i-1)}$

其中，e_i表示节点i的邻点之间存在的边数，k_i表示节点i的邻点的数目，

就表示节点i的邻点之间可能存在的边数。

9.根据权利要求1所述的方法，其特征在于所述计算网络中心性包括：

节点i的网络中心性定义为网络中除i以外的任意两个节点之间通过了节点i的特征路径的数目与该两个节点之间所有特征路径数目的比值，即

$B_i=\underset{i\&NotEqual;j\&NotEqual;k}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&delta;}}_{\mathrm{jk}}\left(i\right)}{{\mathrm{\&delta;}}_{\mathrm{jk}}}$

其中，B_i为节点i的网络中心性，δ_jk为节点j和k之间的最短路径的数目，δ_jk(i)为节点j和k之间的最短路径中经过了节点i的数目。

10.根据权利要求1所述的方法，其特征在于所述自适应提高分类器为6个。

Images