CN103886328A - 基于脑网络模块结构特征的功能磁共振影像数据分类方法 - Google Patents

Abstract

本发明公开了基于脑网络模块结构特征的功能磁共振影像数据分类方法。本发明方法从模块结构角度描述了网络局部聚集特征，体现了网络集团化特性，揭示网络中结构和功能之间的潜在关系，弥补了传统分类方法中对脑局部特性描述的不足，有效的提高了数据分类的准确性。

Description

基于脑网络模块结构特征的功能磁共振影像数据分类方法

技术领域：

本发明涉及的是基于脑网络模块结构特征的功能磁共振影像数据分类方法。

背景技术：

功能磁共振成像（functionalMagneticResonanceImaging,fMRI）是一门图像技术，由于其具有无创伤性、高空间分辩率、使用相对简单等特点，很快被研究人员应用于神经科学，心理学的研究中，并取得了突破性的进展。fMRI主要通过测量血氧水平依赖信号(BloodOxygenationLevelDependent，BOLD)进行脑激活的研究。BLOD主要检测人脑内部血氧的变化，当激发大脑神经系统活动后，会造成脑部部分区域中的血液氧含量发生变化，血氧变化会引起磁场发生局部性改变，进而引起核磁共振信号发生改变。在人脑处于不同状态下，如任务刺激或病变时，脑部功能磁共振图像会出现相应的变化。因此，利用图像分类技术，自动判别图像的类别或状态，是计算机辅助分析的一个重要应用。

人脑是现实世界中最为复杂的网络系统之一。其复杂性不仅体现在数以亿记的神经元及连接的数量，更体现其在不同尺度下的构成以及这些连接在认知功能、思想、感觉及行为时所表现出来的不同模式。近年来，将复杂网络理论应用在神经认知科学中，利用复杂网络基本原理等方法进行属性分析，以期发现网络基本属性及节点间潜在的拓扑关系。复杂网络理论使我们从一个不同的角度来看待人脑这一复杂系统，也为人脑的研究提供了一个新的方向。随着研究的深入，越来越多的研究者将复杂网络应用到脑网络的研究中，从不同角度构建、分析脑网络，并发现无论是结构脑网络，还是功能脑网络，脑网络都具备许多重要的拓扑性质，如“小世界”属性、模块化的组织结构等。模块结构表现了网络的聚集程度或网络集团化程度，是复杂网络的基本属性之一。模块由网络中的部分节点组成，而网络则可划分为若干模块。模块的划分描述了网络局部聚集特性，体现了边的分布的不平衡性，揭示网络了中结构和功能之间的潜在关系。同时，在各类脑疾病的研究中，脑网络方法也得到广泛应用，并取得了重要结论。如：精神分裂症，阿尔茨海默尔症，癫痫，多动症，中风等。结果表明，各种神经疾病状态下，患者的脑网络均存在不同程度的拓扑属性异常。上述结果充分说明，脑网络具备较强的临床价值，其作为重要的分析工具，为脑疾病的诊断提供了新的思路。

传统的功能磁共振影像数据分类方法多为数据驱动方法，主要包括两种。一种为基于体素的分析方法，一种是基于脑网络节点指标的分类方法。前者以体素作为图像分类的基本单位，对该体素在特定时间段内的BOLD信号波动进行分析，并提取特征包括峰值、峰值出现时间点、累积差等，采用复杂的非线性配准，以最大限度地实现个体间的精确对应。这种方法假设目标与样本的内部组织结构是一一对应的。但这种假设并不一定完全合理，人在不同状态下所表现出的BOLD信号波动会受到很多方面因素的影响。后者对影像进行网络建模，并对网络节点指标进行计算及分析，寻找组间显著差异指标并作为分类特征，将其应用到分类模型中。但这种方法只考虑了节点的局部特性，较为局限，并没有考虑到节点间的关联。所以，传统的分类方法不能根据脑的固有属性对磁共振影像进行分类，因此其分类效果并不理想。

发明内容：

本发明的目的在于克服现有技术的不足，提供基于脑网络模块结构特征的功能磁共振影像数据分类方法。

为了解决背景技术所存在的问题，本发明采用以下技术方案：

基于脑网络模块结构特征的功能磁共振影像数据分类方法，它包括以下步骤：

步骤S1：对静息态功能磁共振影像进行预处理，然后根据选定的标准化脑图谱对影像进行区域分割，最后对所分割的各脑区进行平均时间序列的提取；

步骤S2：利用偏相关方法，计算各脑区平均时间序列两两间的关联程度，得到偏相关关联矩阵；

步骤S3：设定特定阈值，将关联矩阵转化为二值矩阵，完成静息态功能脑网络的构建；

步骤S4：对网络进行模块结构分析，划分对应模块并计算相关指标，包括模块度、模块内度及参与系数等；

步骤S5：利用支持向量机SVM分类算法，选择网络模块化指标为分类特征，对图象进行分类，并采用交叉验证方法对分构建的分类器进行检验；

步骤S6：利用敏感性分析方法，量化所选特征在分类模型中的重要程度，对所选特征进行二次筛选，优化模型。

进一步的，步骤S1中，

所述的对于静息态功能磁共振影像数据进行预处理，要最大程度的去除采集过程中由于设备、被试头动、生物噪音等带来的噪音信号，保留功能影像的细节，提高了信噪比；并通过仿射变换和局部非线性变换方法将将影响同一化到所选取的标准空间；预处理步骤至少包括时间层校正、头动校正、联合配准、空间标准化及低频滤波；

所述的选择的标准化脑图谱，是国际通用的解剖标记模板，模板将全脑共分割90个区域，其中左右半脑各45个区域；

所述提取各脑区的平均时间序列，步骤包括：提取每个脑区所包含的所有体素在不同时间点上的激活信号，再将各体素在不同时间点上的激活信号进行算术平均，得到脑区的平均时间序列；激活信号指的是不同时间点上的血氧水平依赖强度。

进一步的，所述步骤S2中包括步骤：

S21：计算平均时间序列两两之间的协方差系数，得到协方差矩阵：

根据步骤S1所提取的平均时间序列，计算各平均时间序列两两之间的协方差矩阵S，元素Si_j为脑区i与脑区j的时间序列间的协方差系数:

$S_{\mathrm{ij}}=\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为时间点个数，xi(t)(i=1,2,……,M)为第i个脑区的时间序列，为第i个脑区的时间序列的平均值，协方差矩阵S的维度为90×90；

S22：计算平均时间序列两两之间的偏相关系数，得到偏相关关联矩阵：

根据S21所计算的协方差系数矩阵S，计算偏相关系数矩阵R，元素r_ij为脑区i与脑区j的偏相关系:

$r_{\mathrm{ij}}=-\frac{s_{\mathrm{ij}}^{-1}}{\sqrt{s_{\mathrm{ii}}^{-1}{s}_{\mathrm{jj}}^{-1}}}$

其中，为协方差矩阵S的逆矩阵的第i行第j列的元素。偏相关系数矩阵R的维度为90×90。

进一步的，所述步骤S3中，设定阈值τ，若偏相关关联矩阵元素r_ij的值大于τ，则二值矩阵B对应元素b_ij的值设为1，表示脑区i和脑区j之间存在连接，否则设为0，表示脑区i和脑区j之间不存在连接；二值化过程可描述为：

$b_{\mathrm{ij}}=\left\{\begin{array}{cc}1,& r_{\mathrm{ij}}\≥\mathrm{\τ}\\ 0,& r_{\mathrm{ij}}\≤\mathrm{\τ}\end{array}\right.$

其中b_ij为二值矩阵中第i行第j列元素，r_ij为偏相关关联矩阵中第i行第j列元素。二值矩阵B为静息态功能脑网络模型，维度为90×90。

进一步的，所述步骤S4中包括：

步骤S41：利用基于贪婪思想的CNM算法，对步骤S3所得到的静息态功能脑网络模型进行块划分：所述的基于贪婪思想的CNM算法，利用堆结构对模块度进行计算、存储及更新；构造稀疏矩阵△Q，通过对该矩阵的更新，以获得模块度最优，并记录此时的模块结构；步骤包括：

（1）首先初始化：

将网络中每个节点均视为一个独立的模块。模块度增量矩阵△Q初始化定义为：

其中，m为网络中的总边数，k_i为节点i的度。由上式可完成模块度增量矩阵△Q的初始化，并得到对应的最大堆H；

（2）然后计算模块度Q：

扫描当前最大堆H，选择其中最大△Qi_j并合并对应的模块i和j。计算新的模块度Q=Q+△Qi_j。更新矩阵△Q、最大堆H和辅助向量a。△Q更新规则可描述为：

更新过程中，当模块度增量矩阵△Q中的最大元素小于零时，则意味着此时网络划分的模块度Q已达到峰值，停止迭代，记录此时的模块划分结果；

（3）遍历下一节点：重复上述步骤，直至网络中所有节点均完成划分；

步骤S42：计算静息态功能脑网络模型的模块度

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设将网络划分为k个模块，定义一个k×k维的对称矩阵E＝e_ij，其中元素e_ij表示网络中连接两个不同模块的节点的边在所有边中所占的比例，这两个节点分别位于第i个模块和第j个模块；设矩阵中对角线上各元素之和为：

它给出了网络中连接某个模块内部各节点的边占所有边的数目中所占比例；定义每行（或者列）中各元素之和为它表示与第i个模块中的节点相连的边在所有边中所占的比例。在此基础上，模块度可被描述为：

$Q=\frac{1}{2m}\underset{s}{\mathrm{\Σ}}(e_{\mathrm{ss}}-a_s^2)$

其中，e_ss表示第s个模块内部边的权重之和占网络中所有边权重总和的比例；

表示网络中全部边的总权重；es_s-a_s ²对应着第s个社区对Q的贡献，Q值越大模块结构越明显；上式的物理意义是，网络中某模块结构中连接两个节点的边的比例，减去在同样的模块结构下任意连接这两个节点的边的比例的期望值；Q的上限为1，而Q越接近这个值，就说明模块结构越明显。

步骤S43：计算静息态功能脑网络模型中各脑区的模块内度。

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的模块内度为

$Z_i=\frac{k_{n_i}-\stackrel{\‾}{k_n}}{{\mathrm{\σ}}_{k_n}}$

其中，

是模块n中的节点i与在本模块中其他节点的连接边的数量；

为模块n中

(i=1,2,…,N)的均值；

为模块n中

(i=1,2,…,N)的标准差；模块内度衡量了节点与其他同属于同一模块的节点的连接数量；当某个节点在模块内相比较其他节点具有更高的模块内连接数量时，该节点的模块内度就相对较高；

步骤S44：计算静息态功能脑网络模型中各脑区的参与系数。

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的参与系数为：

$P_i=1-\underset{n=1}{\overset{N_m}{\mathrm{\Σ}}}{\left(\frac{k_{n_i}}{k_i}\right)}^2$

其中，Kn_i是模块n中节点i的模块内度。K_i是节点i的总度数。N_m是模块总数。参与系数衡量了模块之间的连接情况。若模块n中的节点i与其他模块有大量的连接，则参与系数接近1，否则参与系数的值接近0。

进一步的，步骤S5中所述采用交叉验证方法对分构建的分类器进行检验：

每次验证过程将样本集中随机选择70%的样本作为训练样本，剩余30%为测试样本，进行分类测试，得到测试结果；重复进行100次测试后的分类准确率的算术平均值，作为分类器的分类正确率；

进一步的，步骤S6中所述的利用敏感性分析方法，计算每个特征在目标类别中的方差变化，以判断该特征在分类过程中的重要性并标准化及进行排名，量化所选特征在分类模型中的重要程度；步骤包括：X_j为特征,j=1,……k，k为特征总数，Y为基于特征X₁到X_k的模型。特征i的敏感性为：

$S_i=\frac{V_i}{V\left(Y\right)}=\frac{V\left(E\left(Y\right|X_i)\right)}{V\left(Y\right)}$

其中，V(Y)为绝对输出变量。数学期望E要求X-i的积分，即除Xi外的所有元素。然后，计算标准化特征敏感性，用以表征特征在分类中的重要程度。标准化的特征敏感性可描述为：

$I_i=\frac{S_i}{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{S}_j}$

其中，S_i为特征i的敏感性，k为特征总数。

本发明对比现有技术，有如下的有益效果：本发明提供的基于脑网络模块结构特征的功能磁共振影像数据分类方法，从模块结构角度描述了网络局部聚集特征，体现了网络集团化特性，揭示网络中结构和功能之间的潜在关系，弥补了传统分类方法中对脑局部特性描述的不足，有效的提高了数据分类的准确性。

附图说明：

图1是本发明方法流程图。

图2是本发明的方法与基于体素的分类方法及基于网络节点特征的分类方法性能比较结果图。

具体实施方式：

下面结合附图和具体实施方式对本发明作进一步描述：

本发明提出的一种基于脑网络模块结构特征的功能磁共振影像数据分类方法，利用脑网络模块结构，分析网络局部聚集特性，揭示结构和功能之间的潜在关系，有效的提高了数据分类的准确性。

本发明的基于脑网络模块结构特征的功能磁共振影像数据分类方法的具体实施过程如图1所示，包括如下步骤：

步骤S1：对静息态功能磁共振影像进行预处理，然后根据选定的标准化脑图谱对影像进行区域分割，最后对所分割的各脑区进行平均时间序列的提取。具体包括如下步骤：

步骤S11：对静息态功能磁共振影像数据进行预处理；

图像采集过程中，不可避免会产生各种类型的噪音，包括被试头动、呼吸、心跳等。这些噪音的存在对图像的质量产生重要的影响，所以对于原始数据需进行去噪处理，提高信噪比。图像预处理采用SPM（StatisticalParametricMapping）软件进行，主要步骤包括时间层校正、头动校正、联合配准、空间标准化及低频滤波等。

步骤S12：利用所选定的标准化脑图谱，对预处理后的影像数据进行区域分割；

采用国际通用的标准脑解剖图谱AAL图谱（AutomatedAnatomicalLabeling）。该图谱对全脑进行了区域级的大尺度节点分割，共定义90个脑区，左右半脑各45个。AAL图谱是相关研究领域广泛使用的脑图谱。

步骤S13：对所划分的各脑区，进行平均时间序列的提取。

计算每个脑区中所包含的所有体素的BOLD信号的算术平均值来表征该节点的信号值。提取每个脑区所包含的所有体素在不同时间点上的激活信号，再将各体素在不同时间点上的激活信号进行算术平均，得到脑区的平均时间序列。激活信号指的是不同时间点上的BOLD信号强度。

步骤S2：利用偏相关方法，计算各脑区平均时间序列两两间的关联程度，得到偏相关关联矩阵。具体包括如下步骤：

步骤S21：计算平均时间序列两两之间的协方差系数，得到协方差矩阵；

根据步骤S1所提取的平均时间序列，计算各平均时间序列两两之间的协方差矩阵S，元素Si_j为脑区i与脑区j的时间序列间的协方差系数，

$S_{\mathrm{ij}}=\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,2,……,M)为第i个脑区的时间序列，

为第i个脑区的时间序列的平均值，协方差矩阵S的维度为90×90。

步骤S22：计算平均时间序列两两之间的偏相关系数，得到偏相关关联矩阵。

根据S21所计算的协方差系数矩阵S，计算偏相关系数矩阵R，元素ri_j为脑区i与脑区j的偏相关系数，

$r_{\mathrm{ij}}=-\frac{s_{\mathrm{ij}}^{-1}}{\sqrt{s_{\mathrm{ii}}^{-1}{s}_{\mathrm{jj}}^{-1}}}$

其中，

为协方差矩阵S的逆矩阵的第i行第j列的元素。偏相关系数矩阵R的维度为90×90。

步骤S3：特定阈值，将关联矩阵转化为二值矩阵，完成静息态功能脑网络的构建。

设定阈值τ，若偏相关关联矩阵元素r_ij的值大于τ，则二值矩阵B对应元素b_ij的值设为1，表示脑区i和脑区j之间存在连接，否则设为0，表示脑区i和脑区j之间不存在连接。二值化过程可描述为：

$b_{\mathrm{ij}}=\left\{\begin{array}{cc}1,& r_{\mathrm{ij}}\≥\mathrm{\τ}\\ 0,& r_{\mathrm{ij}}\≤\mathrm{\τ}\end{array}\right.$

其中b_ij为二值矩阵中第i行第j列元素，r_ij为偏相关关联矩阵中第i行第j列元素。二值矩阵B为静息态功能脑网络模型，维度为90×90。

步骤S4：对网络进行模块结构分析，划分对应模块并计算相关指标，包括模块度、模块内度及参与系数等。

步骤S41：利用基于贪婪思想的CNM算法，对步骤S3所得到的静息态功能脑网络模型进行模块划分。

步骤S41所述的基于贪婪思想的CNM算法，利用堆结构对模块度进行计算、存储及更新。构造稀疏矩阵△Q，通过对该矩阵的更新，以获得模块度最优，并记录此时的模块结构。

初始化。将网络中每个节点均视为一个独立的模块。模块度增量矩阵△Q初始化定义为：

其中，m为网络中的总边数，k_i为节点i的度。由上式可完成模块度增量矩阵△Q的初始化，并得到对应的最大堆H。

计算模块度Q。扫描当前最大堆H，选择其中最大△Q_ij并合并对应的模块i和j。计算新的模块度Q=Q+△Q_ij。更新矩阵△Q、最大堆H和辅助向量a。△Q更新规则可描述为：

更新过程中，当模块度增量矩阵△Q中的最大元素小于零时，则意味着此时网络划分的模块度Q已达到峰值，停止迭代，记录此时的模块划分结果。

遍历下一节点。重复上述步骤，直至网络中所有节点均完成划分。

步骤S42：计算静息态功能脑网络模型的模块度。

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设将网络划分为k个模块。定义一个k×k维的对称矩阵E＝e_ij，其中元素e_ij表示网络中连接两个不同模块的节点的边在所有边中所占的比例，这两个节点分别位于第i个模块和第j个模块。设矩阵中对角线上各元素之和为：

它给出了网络中连接某个模块内部各节点的边占所有边的数目中所占比例；定义每行（或者列）中各元素之和为

它表示与第i个模块中的节点相连的边在所有边中所占的比例。在此基础上，模块度可被描述为：

$\frac{1}{2m}\underset{s}{\mathrm{\Σ}}(e_{\mathrm{ss}}-a_s^2)$

其中，e_ss表示第s个模块内部边的权重之和占网络中所有边权重总和的比例；

表示网络中全部边的总权重。e_ss-a_s ²对应着第s个社区对Q的贡献，Q值越大模块结构越明显。上式的物理意义是，网络中某模块结构中连接两个节点的边的比例，减去在同样的模块结构下任意连接这两个节点的边的比例的期望值。Q的上限为1，而Q越接近这个值，就说明模块结构越明显。

步骤S43：计算静息态功能脑网络模型中各脑区的模块内度。

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的模块内度为

$Z_i=\frac{k_{n_i}-\stackrel{\‾}{k_n}}{{\mathrm{\σ}}_{k_n}}$

其中，k_ni是模块n中的节点i与在本模块中其他节点的连接边的数量。

为模块n中

(i=1,2,…,N)的均值。为模块n中

(i=1,2,…,N)的标准差。模块内度衡量了节点与其他同属于同一模块的节点的连接数量。当某个节点在模块内相比较其他节点具有更高的模块内连接数量时，该节点的模块内度就相对较高。

步骤S44：计算静息态功能脑网络模型中各脑区的参与系数。

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的参与系数为

$P_i=1-\underset{n=1}{\overset{N_m}{\mathrm{\Σ}}}{\left(\frac{k_{n_i}}{k_i}\right)}^2$

其中，是模块n中节点i的模块内度。k_i是节点i的总度数。N_m是模块总数。参与系数衡量了模块之间的连接情况。若模块n中的节点i与其他模块有大量的连接，则参与系数接近1，否则参与系数的值接近0。

步骤S5：利用支持向量机分类算法，选择网络模块化指标为分类特征，对图象进行分类，并采用交叉验证方法对分构建的分类器进行检验。

每次验证过程将样本集中随机选择70%的样本作为训练样本，剩余30%为测试样本，进行分类测试，得到测试结果；重复进行100次测试后的分类准确率的算术平均值，作为分类器的分类正确率。

步骤S6：利用敏感性分析方法，量化所选特征在分类模型中的重要程度。

步骤S6所述的利用敏感性分析方法，计算每个特征在目标类别中的方差变化，以判断该特征在分类过程中的重要性并标准化及进行排名，量化所选特征在分类模型中的重要程度。步骤包括：X_j为特征,j=1,……k，k为特征总数，Y为基于特征X₁到X_k的模型。特征i的敏感性为：

$S_i=\frac{V_i}{V\left(Y\right)}=\frac{V\left(E\left(Y\right|X_i)\right)}{V\left(Y\right)}$

其中，V(Y)为绝对输出变量。数学期望E要求X-i的积分，即除X_i外的所有元素。然后，计算标准化特征敏感性，用以表征特征在分类中的重要程度。标准化的特征敏感性可描述为：

$I_i=\frac{S_i}{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{S}_j}$

其中，S_i为特征i的敏感性，k为特征总数。

将本发明的方法与基于体素的分类方法及基于网络节点特征的分类方法，在不同的特征数目选择下，进行性能比较，结果如图2所示。结果表明，在选择10个特征数条件下，三种方法正确率接近。随着特征数目的增多，本发明方法及基于网络节点特征方法，较基于体素方法，正确率呈显著增加。在50个特征数目之前，本发明方法较其他两种方法，均表现出最优。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，它包括以下步骤：

步骤S1：对静息态功能磁共振影像进行预处理，然后根据选定的标准化脑图谱对影像进行区域分割，最后对所分割的各脑区进行平均时间序列的提取；

步骤S2：利用偏相关方法，计算各脑区平均时间序列两两间的关联程度，得到偏相关关联矩阵；

步骤S3：设定特定阈值，将关联矩阵转化为二值矩阵，完成静息态功能脑网络的构建；

步骤S4：对网络进行模块结构分析，划分对应模块并计算相关指标，包括模块度、模块内度及参与系数等；

步骤S5：利用支持向量机SVM分类算法，选择网络模块化指标为分类特征，对图象进行分类，并采用交叉验证方法对分构建的分类器进行检验；步骤S6：利用敏感性分析方法，量化所选特征在分类模型中的重要程度，对所选特征进行二次筛选，优化模型。

2.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，步骤S1中，

所述的对于静息态功能磁共振影像数据进行预处理，要最大程度的去除采集过程中由于设备、被试头动、生物噪音等带来的噪音信号，保留功能影像的细节，提高了信噪比；并通过仿射变换和局部非线性变换方法将将影响同一化到所选取的标准空间；预处理步骤至少包括时间层校正、头动校正、联合配准、空间标准化及低频滤波；

所述的选择的标准化脑图谱，是国际通用的解剖标记模板，模板将全脑共分割90个区域，其中左右半脑各45个区域；

所述提取各脑区的平均时间序列，步骤包括：提取每个脑区所包含的所有体素在不同时间点上的激活信号，再将各体素在不同时间点上的激活信号进行算术平均，得到脑区的平均时间序列；激活信号指的是不同时间点上的血氧水平依赖强度。

3.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，所述步骤S2中包括步骤：

S21：计算平均时间序列两两之间的协方差系数，得到协方差矩阵：

根据步骤S1所提取的平均时间序列，计算各平均时间序列两两之间的协方差矩阵S，元素S_ij为脑区i与脑区j的时间序列间的协方差系数:

$S_{\mathrm{ij}}=\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为时间点个数，xi(t)(i=1,2,……,M)为第i个脑区的时间序列，为第i个脑区的时间序列的平均值，协方差矩阵S的维度为90×90；

S22：计算平均时间序列两两之间的偏相关系数，得到偏相关关联矩阵：

根据S21所计算的协方差系数矩阵S，计算偏相关系数矩阵R，元素r_ij为脑区i与脑区j的偏相关系数:

$r_{\mathrm{ij}}=-\frac{s_{\mathrm{ij}}^{-1}}{\sqrt{s_{\mathrm{ii}}^{-1}{s}_{\mathrm{jj}}^{-1}}}$

其中，为协方差矩阵S的逆矩阵的第i行第j列的元素；偏相关系数矩阵R的维度为90×90。

4.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，所述步骤S3中，

设定阈值τ，若偏相关关联矩阵元素r_ij的值大于τ，则二值矩阵B对应元素b_ij的值设为1，表示脑区i和脑区j之间存在连接，否则设为0，表示脑区i和脑区j之间不存在连接；二值化过程可描述为：

$b_{\mathrm{ij}}=\left\{\begin{array}{cc}1,& r_{\mathrm{ij}}\≥\mathrm{\τ}\\ 0,& r_{\mathrm{ij}}\≤\mathrm{\τ}\end{array}\right.$

其中b_ij为二值矩阵中第i行第j列元素，r_ij为偏相关关联矩阵中第i行第j列元素；二值矩阵B为静息态功能脑网络模型，维度为90×90。

5.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，所述步骤S4中包括：

步骤S41：利用基于贪婪思想的CNM算法，对步骤S3所得到的静息态功能脑网络模型进行模块划分：所述的基于贪婪思想的CNM算法，利用堆结构对模块度进行计算、存储及更新；构造稀疏矩阵△Q，通过对该矩阵的更新，以获得模块度最优，并记录此时的模块结构；步骤包括：

（1）首先初始化：

将网络中每个节点均视为一个独立的模块；模块度增量矩阵△Q初始化定义为：

其中，m为网络中的总边数，k_i为节点i的度；由上式可完成模块度增量矩阵△Q的初始化，并得到对应的最大堆H；

（2）然后计算模块度Q：

扫描当前最大堆H，选择其中最大△Q_ij并合并对应的模块i和j；计算新的模块度Q=Q+△Q_ij；更新矩阵△Q、最大堆H和辅助向量a；△Q更新规则可描述为：

更新过程中，当模块度增量矩阵△Q中的最大元素小于零时，则意味着此时网络划分的模块度Q已达到峰值，停止迭代，记录此时的模块划分结果；

（3）遍历下一节点：重复上述步骤，直至网络中所有节点均完成划分；

步骤S42：计算静息态功能脑网络模型的模块度

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设将网络划分为k个模块，定义一个k×k维的对称矩阵E＝e_ij，其中元素e_ij表示网络中连接两个不同模块的节点的边在所有边中所占的比例，这两个节点分别位于第i个模块和第j个模块；设矩阵中对角线上各元素之和为：

它给出了网络中连接某个模块内部各节点的边占所有边的数目中所占比例；定义每行（或者列）中各元素之和为

它表示与第i个模块中的节点相连的边在所有边中所占的比例；在此基础上，模块度可被描述为：

$\frac{1}{2m}\underset{s}{\mathrm{\Σ}}(e_{\mathrm{ss}}-a_s^2)$

其中，e_ss表示第s个模块内部边的权重之和占网络中所有边权重总和的比例；

表示网络中全部边的总权重；e_ss-a_s ²对应着第s个社区对Q的贡献，Q值越大模块结构越明显；上式的物理意义是，网络中某模块结构中连接两个节点的边的比例，减去在同样的模块结构下任意连接这两个节点的边的比例的期望值；Q的上限为1，而Q越接近这个值，就说明模块结构越明显；

步骤S43：计算静息态功能脑网络模型中各脑区的模块内度；

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的模块内度为

$Z_i=\frac{k_{n_i}-\stackrel{\‾}{k_n}}{{\mathrm{\σ}}_{k_n}}$

其中，Kn_i是模块n中的节点i与在本模块中其他节点的连接边的数量；

为模块n中Kn_i(i=1,2,…,N)的均值；

为模块n中Kni(i=1,2,…,N)的标准差；模块内度衡量了节点与其他同属于同一模块的节点的连接数量；当某个节点在模块内相比较其他节点具有更高的模块内连接数量时，该节点的模块内度就相对较高；

步骤S44：计算静息态功能脑网络模型中各脑区的参与系数；

根据步骤S41所得到的静息态功能脑网络模块划分结果，假设节点i属于模块n，节点i的参与系数为：

$P_i=1-\underset{n=1}{\overset{N_m}{\mathrm{\Σ}}}{\left(\frac{k_{n_i}}{k_i}\right)}^2$

其中，Kn_i是模块n中节点i的模块内度；K_i是节点i的总度数；N_m是模块总数；参与系数衡量了模块之间的连接情况；若模块n中的节点i与其他模块有大量的连接，则参与系数接近1，否则参与系数的值接近0。

6.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，步骤S5中所述采用交叉验证方法对分构建的分类器进行检验：每次验证过程将样本集中随机选择70%的样本作为训练样本，剩余30%为测试样本，进行分类测试，得到测试结果；重复进行100次测试后的分类准确率的算术平均值，作为分类器的分类正确率。

7.根据权利要求1所述的基于脑网络模块结构特征的功能磁共振影像数据分类方法，其特征在于，步骤S6中所述的利用敏感性分析方法，计算每个特征在目标类别中的方差变化，以判断该特征在分类过程中的重要性并标准化及进行排名，量化所选特征在分类模型中的重要程度；步骤包括：X_j为特征,j=1,……k，k为特征总数，Y为基于特征X₁到X_k的模型；特征i的敏感性为：

$S_i=\frac{V_i}{V\left(Y\right)}=\frac{V\left(E\left(Y\right|X_i)\right)}{V\left(Y\right)}$

其中，V(Y)为绝对输出变量；数学期望E要求X-i的积分，即除Xi外的所有元素；然后，计算标准化特征敏感性，用以表征特征在分类中的重要程度；标准化的特征敏感性可描述为：

$I_i=\frac{S_i}{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{S}_j}$

其中，S_i为特征i的敏感性，k为特征总数。

Images