基于主成分分析的关键脑区的度量方法

摘要

本发明属于神经影像数据分析技术领域，提出了一种基于主成分分析的关键脑区的度量方法，将采集M个人静息态下的fMRI数据并对其进行预处理；构建Q个稀疏度下的脑网络；计算每个稀疏度下脑网络加权的节点的度、节点效率和介数中心度；利用主成分分析将得到的节点的度、节点效率和介数中心度进行综合降维得到一种新指标，即每个稀疏度下每个人每个节点的中心度得分H；再先后综合Q个稀疏度和M个人的中心度得分，得到每个节点的中心度得分，并根据其大小判定关键脑区。本发明用于度量脑区的重要性，克服了确定功能脑网络的稀疏度主观选择问题，结合多种经典脑区重要性度量方法，更客观、更精确地度量关键脑区。

权利要求

1.一种基于主成分分析的关键脑区的度量方法，其特征在于，包括以下步骤：步骤1：采集M个人静息态下的fMRI数据；

步骤2：对每个fMRI数据进行预处理；

步骤3：基于预处理后的fMRI数据，构建Q个稀疏度下的脑网络；

步骤4：计算每个稀疏度下脑网络的节点的度、节点效率和介数中心度；具体操作包括：步骤4.1：在Q个稀疏度下，计算每个脑网络的节点的度；对于脑网络中任意一个节点i，其节点的度计算公式为：其中，wi是对第i个节点的度的加权值，即第i个节点对应脑区的体积大小占整个大脑的比重；bij是指邻接矩阵B中第(i，j)个元素值；N是指每个脑网络中的节点数目；j是指脑网络中第j个节点；所述邻接矩阵B中的元素bij＝0或1，0是指两个节点之间不存在边，1是指两个节点之间存在边，两个节点之间存在的边即构成每个稀疏度下脑网络的边；公式4.1中，wi的计算公式如下：其中，ni是第i个节点包含的体素数目；

步骤4.2：在Q个稀疏度下，计算每个脑网络的节点效率；

对于脑网络中任意一个节点i，其节点效率计算公式为：其中，fij是节点i到节点j的最短路径的倒数；

步骤4.3：在Q个稀疏度下，计算每个脑网络节点的介数中心度；

对于脑网络中任意一个节点i，其介数中心度计算公式为：其中，σjk表示从节点j到节点k最短路径的条数，σjk(i)表示在这些最短路径中，并经过节点i的路径条数；

步骤5：利用主成分分析将步骤4中得到的节点的度、节点效率和介数中心度进行合并降维得到一种新指标，即每个稀疏度下每个人每个节点的中心度得分H，具体操作如下：步骤5.1：将步骤4中得到的节点的度、节点效率和介数中心度的值重组构成一个二维矩阵R；

R＝(Df,ijk),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3    (5.1)其中，Df,ijk表示第i个人第j个节点第k个稀疏度下的第f个属性值；矩阵R的行数为M×N×Q，M为人的数目，N为节点的数目，Q为稀疏度的数目；矩阵R的列分别为节点的度、节点效率和介数中心度的值；

步骤5.2：基于主成分分析，对矩阵R进行标准化：其中， 是第f个属性值的均值，Sf是第f个属性值的标准差；

对矩阵R标准化后，得到一个二维矩阵Y，其矩阵形式如下：Y＝(Df,ijk*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3   (5.3)计算矩阵Y中两两列向量间的相关系数：其中， 分别是对矩阵R标准化后的第m个属性值和第n个属性值的均值；

将得到的矩阵Y中两两列向量间的相关系数组成矩阵X，其矩阵形式如下：X＝(xmn),m＝1,2,3；n＝1,2,3     (5.5)利用齐次线性方程组求非零解的方法，得出矩阵X的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X的主特征向量P：其中，λmax和eλmax分别是矩阵X的最大特征值及其对应的单位特征向量；

利用得到的矩阵X的主特征向量P，计算每个稀疏度下每个人每个节点的中心度得分H：H＝R×P   (5.7)

步骤6：利用主成分分析，综合Q个稀疏度下每个人每个节点的中心度得分H，对其降维得到每个人每个节点的中心度得分H1；具体操作如下：步骤6.1：在Q个稀疏度下，将每个人每个节点的中心度得分H重组构成一个二维矩阵R1；

R1＝(Dk,ij),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q   (6.1)其中，Dk,ij表示第k个稀疏度下的第i个人第j个节点的中心度得分，矩阵R1的行数为M×N，列数为Q；

步骤6.2：基于主成分分析，对矩阵R1进行标准化：其中， 是第k个稀疏度下的中心度得分的均值，Sk是第k个稀疏度下的中心度得分的标准差；

对矩阵R1标准化后，得到一个二维矩阵Y1，其矩阵形式如下：Y1＝(Dk,ij*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q     (6.3)计算矩阵Y1中两两列向量间的相关系数：

其中， 分别是对矩阵R1标准化后的第m个稀疏度下和第n个稀疏度下的中心度得分的均值；

将得到的矩阵Y1两两列向量间的相关系数组成一个矩阵X1，其矩阵形式如下：X1＝(xmn),m＝1,2,…,Q；n＝1,2,…,Q   (6.5)利用齐次线性方程组求非零解的方法，得出矩阵X1的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X1的主特征向量P1：其中，λmax和eλmax分别为矩阵X1的最大特征值及其对应的单位特征向量；

利用得到的矩阵X1的主特征向量P1，计算每个人每个节点的中心度得分H1：H1＝R1×P1    (6.7)

步骤7：利用主成分分析，综合M个人每个节点的中心度得分H1，并对其降维得到每个节点的中心度得分H2；具体操作如下：步骤7.1：将M个人的每个节点的中心度得分H1重组构成一个二维矩阵R2；

R2＝(Di,j),i＝1,2,…,M；j＝1,2,…,N    (7.1)其中，Di,j表示第i个人的第j个节点的中心度得分，矩阵R2的行数为N，列数为M；

步骤7.2：基于主成分分析，对矩阵R2进行标准化：其中， 是第i个人的中心度得分的均值，Si是第i个人的中心度得分的标准差；

对矩阵R2标准化后，得到一个二维矩阵Y2，其矩阵形式如下：Y2＝(Di,j*),i＝1,2,…,M；j＝1,2,…,N    (7.3)计算矩阵Y2中两两列向量间的相关系数：

其中， 分别是对矩阵R2标准化后的第m个人和第n个人的中心度得分的均值；

将得到的矩阵Y2中两两列向量间的相关系数组成矩阵X2，其矩阵形式如下：X2＝(xmn),m＝1,2,…,M；n＝1,2,…,M    (7.5)利用齐次线性方程组求非零解的方法，得出矩阵X2的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，计算矩阵X2的主特征向量P2：其中，λmax和eλmax分别是矩阵X2的最大特征值及其对应的单位特征向量；

利用得到的矩阵X2的主特征向量P2，计算每个节点的中心度得分H2：H2＝R2×P2    (7.7)

步骤8：将步骤7.2得到的每个节点的中心度得分H2进行归一化：选取1个最大的中心度得分H2max，将所得的中心度得分H2分别除以H2max，得到H2＇，0＜H2＇≤1；对所得H2＇进行降序排列，前C个H2＇对应的脑区判定为关键脑区。

2.根据权利要求1所述的基于主成分分析的关键脑区的度量方法，其特征在于，所述的对每个fMRI数据进行预处理，具体包括以下步骤：对采集到的fMRI数据去除时间点；

对去除时间点后的fMRI数据进行时间层校正；

对时间层校正后的fMRI数据进行头动校正；

对头动校正后的fMRI数据进行空间标准化；

对空间标准化后的fMRI数据进行去线性漂移；

对去线性漂移后的fMRI数据进行带通滤波。

3.根据权利要求2所述的基于主成分分析的关键脑区的度量方法，其特征在于，所述的对采集到的fMRI数据去除时间点，包括：去除时间点的个数为10；

所述的对去线性漂移后的fMRI数据进行带通滤波，包括：带通滤波的频率区间为0.01～0.08Hz。

4.根据权利要求1所述的基于主成分分析的关键脑区的度量方法，其特征在于，所述的步骤3，具体包括以下步骤：步骤3.1：

选定一个脑模板和Q个稀疏度；

根据选定的脑模板将每个fMRI数据划分为N个脑区；

定义每个脑区作为每个稀疏度下用于构建脑网络的节点，每个脑网络共有N个节点；

步骤3.2：提取每个节点的平均时间序列信号τi；

其中，τi是指第i个节点的平均时间序列，Gj是指第i个节点中第j个体素的时间序列，ni是第i个节点包含的体素数目；

步骤3.3：将步骤3.2中每个节点的平均时间序列信号τi组成矩阵T，T＝(τij)，并计算出两两节点之间的相关系数；

节点i和节点j之间的相关系数aij的计算公式：其中， 和 分别是i脑区和j脑区时间序列的均值；

步骤3.4：将步骤3.3中计算得到的相关系数组成一个N×N的相关系数矩阵A，A＝(aij)；

对相关系数矩阵A二值化：对矩阵A中的元素，利用Q个稀疏度确定值为1的元素个数；对一个稀疏度q％，取矩阵A中元素总个数的前q％个数值大的元素，将其值转化为1，其余的元素值转化为0，从而得到一个邻接矩阵B，矩阵B＝(bij)；

步骤3.5：根据步骤3.1和3.4中的节点和边，在Q个稀疏度下构建Q个脑网络。

5.根据权利要求4所述的基于主成分分析的关键脑区的度量方法，其特征在于，所述的步骤3.1中，选定的脑模板为AAL模板。

说明书

基于主成分分析的关键脑区的度量方法

技术领域

[0001] 本发明属于神经影像数据分析技术领域，具体涉及静息态下fMRI人脑网络中关键脑区的度量方法，尤其涉及一种基于主成分分析的关键脑区的度量方法。

背景技术

[0002] 功能磁共振成像的原理是基于血氧依赖水平的，当大脑某区域的神经元兴奋，则该区域需要的氧气量会显著增加，从而引起该区域的血流量增加，而氧气需要红细胞蛋白运送，当氧合血红蛋白来到大脑的某个组织，该组织的磁敏感度将会降低，此时磁共振成像显示较高的局部信号，当该组织取走氧气后，氧合血红蛋白会变成脱氧血红蛋白，导致该组织的磁敏感性增强，在磁共振成像中显示较低的局部信号，由于其具有的较高的时间和空间分辨率、无创性和可重复性等特点，因而被广泛的用于人脑功能的研究。

[0003] 图论作为描述网络特征的重要工具，目前已被广泛的用于脑网络的研究。基于图论，可以将人脑看成是一个脑网络，对人脑功能活动的研究转化为对脑网络的定量分析，对脑网络属性的分析在研究大脑功能活动具有非常重要的意义。通过比较脑疾病患者和健康对照人群的脑网络属性，有可能找到一些脑疾病的病因及其治疗方法。

[0004] 衡量人类大脑关键脑区通常有节点的度、节点效率和介数中心度等测量方法。然而这些传统的单一测量指标并未考虑节点本身在大脑中所占的比重，且它们是从不同的角度出发来进行测量，时常会得到不一致的结果。即使是相同的测量指标，在不同阈值下得到的结果也不尽相同。而关键脑区的确定是脑网络分析中的一个重要方面，其关系着脑网络架构的整体变化，与脑疾病的发生发展有重要的关联，是目前全局脑网络研究中亟待解决的关键技术问题。

发明内容

[0005] 本发明的目的在于克服上述现有技术中存在的利用不同指标度量关键脑区，其结果不一致的现象，及稀疏度阈值的主观选择问题，从而提出一种基于主成分分析的关键脑区的度量方法。

[0006] 为了便于理解，对本发明中出现的部分名词作以下解释说明：

[0007] 关键脑区：其改变被认为与某一行为或某一疾病发生发展有关，是度量大脑网络改变的重要指标。

[0008] 静息态：指人类在清醒、不受任何刺激、不集中做任何事情的放松状态。

[0009] fMRI：英文全称为functional magnetic resonance imaging，功能磁共振成像，是一种新兴的神经影像学技术，其原理是利用磁振造影来测量神经元活动所引发之血液动力的改变。

[0010] 稀疏度：是指脑网络中节点之间的连接线数与最大可能存在的边数之比。

[0011] 体素：大脑的每个脑区中含有多个体素，每个体素是指一小块立方体区域，本发明中的每个体素的大小为3×3×3mm3。

[0012] AAL模板：英文全称为Anatomical Automatic Labeling，将人脑分为116个区域，但只有90个区域属于大脑，剩余26个区域属于小脑结构，在此发明中只用到90个大脑区域。

[0013] 为了实现上述目的，本发明采用以下技术方案：

[0014] 一种基于主成分分析的关键脑区的度量方法，包括以下步骤：

[0015] 步骤1：采集M个人静息态下的fMRI数据；

[0016] 步骤2：对每个fMRI数据进行预处理；

[0017] 步骤3：基于预处理后的fMRI数据，构建Q个稀疏度下的脑网络；

[0018] 步骤4：计算每个稀疏度下脑网络的节点的度、节点效率和介数中心度；具体操作包括：

[0019] 步骤4.1：在Q个稀疏度下，计算每个脑网络的节点的度；对于脑网络中任意一个节点i，其节点的度计算公式为：

[0020]

[0021] 其中，wi是对第i个节点的度的加权值，即第i个节点对应脑区的体积大小占整个大脑的比重；bij是指邻接矩阵B中第(i，j)个元素值；N是指每个脑网络中的节点数目；j是指脑网络中第j个节点；所述邻接矩阵B中的元素bij＝0或1，0是指两个节点之间不存在边，1是指两个节点之间存在边，两个节点之间存在的边即构成每个稀疏度下脑网络的边；公式4.1中，wi的计算公式如下：

[0022]

[0023] 其中，ni是第i个节点包含的体素数目；

[0024] 步骤4.2：在Q个稀疏度下，计算每个脑网络的节点效率；

[0025] 对于脑网络中任意一个节点i，其节点效率计算公式为：

[0026]

[0027] 其中，fij是节点i到节点j的最短路径的倒数；

[0028] 步骤4.3：在Q个稀疏度下，计算每个脑网络节点的介数中心度；

[0029] 对于脑网络中任意一个节点i，其介数中心度计算公式为：

[0030]

[0031] 其中，σjk表示从节点j到节点k最短路径的条数，σjk(i)表示在这些最短路径中，并经过节点i的路径条数；

[0032] 步骤5：利用主成分分析将步骤4中得到的节点的度、节点效率和介数中心度进行合并降维得到一种新指标，即每个稀疏度下每个人每个节点的中心度得分H，具体操作如下：

[0033] 步骤5.1：将步骤4中得到的节点的度、节点效率和介数中心度的值重组构成一个二维矩阵R；

[0034] R＝(Df,ijk),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.1)[0035] 其中，Df,ijk表示第i个人第j个节点第k个稀疏度下的第f个属性值；矩阵R的行数为M×N×Q，M为人的数目，N为节点的数目，Q为稀疏度的数目；矩阵R的列分别为节点的度、节点效率和介数中心度的值；

[0036] 步骤5.2：基于主成分分析，对矩阵R进行标准化：

[0037]

[0038] 其中， 是第f个属性值的均值，Sf是第f个属性值的标准差；

[0039] 对矩阵R标准化后，得到一个二维矩阵Y，其矩阵形式如下：

[0040] Y＝(Df,ijk*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.3)[0041] 计算矩阵Y中两两列向量间的相关系数：

[0042]

[0043] 其中， 分别是对矩阵R标准化后的第m个属性值和第n个属性值的均值；

[0044] 将得到的矩阵Y中两两列向量间的相关系数组成矩阵X，其矩阵形式如下：

[0045] X＝(xmn),m＝1,2,3；n＝1,2,3  (5.5)

[0046] 利用齐次线性方程组求非零解的方法，得出矩阵X的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X的主特征向量P：

[0047]

[0048] 其中，λmax和eλmax分别是矩阵X的最大特征值及其对应的单位特征向量；

[0049] 利用得到的矩阵X的主特征向量P，计算每个稀疏度下每个人每个节点的中心度得分H：

[0050] H＝R×P  (5.7)

[0051] 步骤6：利用主成分分析，综合Q个稀疏度下每个人每个节点的中心度得分H，对其降维得到每个人每个节点的中心度得分H1；具体操作如下：

[0052] 步骤6.1：在Q个稀疏度下，将每个人每个节点的中心度得分H重组构成一个二维矩阵R1；

[0053] R1＝(Dk,ij),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.1)

[0054] 其中，Dk,ij表示第k个稀疏度下的第i个人第j个节点的中心度得分，矩阵R1的行数为M×N，列数为Q；

[0055] 步骤6.2：基于主成分分析，对矩阵R1进行标准化：

[0056]

[0057] 其中， 是第k个稀疏度下的中心度得分的均值，Sk是第k个稀疏度下的中心度得分的标准差；

[0058] 对矩阵R1标准化后，得到一个二维矩阵Y1，其矩阵形式如下：

[0059] Y1＝(Dk,ij*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.3)

[0060] 计算矩阵Y1中两两列向量间的相关系数：

[0061]

[0062] 其中， 分别是对矩阵R1标准化后的第m个稀疏度下和第n个稀疏度下的中心度得分的均值；

[0063] 将得到的矩阵Y1两两列向量间的相关系数组成一个矩阵X1，其矩阵形式如下：

[0064] X1＝(xmn),m＝1,2,…,Q；n＝1,2,…,Q  (6.5)

[0065] 利用齐次线性方程组求非零解的方法，得出矩阵X1的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X1的主特征向量P1：

[0066]

[0067] 其中，λmax和eλmax分别为矩阵X1的最大特征值及其对应的单位特征向量；

[0068] 利用得到的矩阵X1的主特征向量P1，计算每个人每个节点的中心度得分H1：

[0069] H1＝R1×P1  (6.7)

[0070] 步骤7：利用主成分分析，综合M个人每个节点的中心度得分H1，并对其降维得到每个节点的中心度得分H2；具体操作如下：

[0071] 步骤7.1：将M个人的每个节点的中心度得分H1重组构成一个二维矩阵R2；

[0072] R2＝(Di,j),i＝1,2,…,M；j＝1,2,…,N  (7.1)

[0073] 其中，Di,j表示第i个人的第j个节点的中心度得分，矩阵R2的行数为N，列数为M；

[0074] 步骤7.2：基于主成分分析，对矩阵R2进行标准化：

[0075]

[0076] 其中， 是第i个人的中心度得分的均值，Si是第i个人的中心度得分的标准差；

[0077] 对矩阵R2标准化后，得到一个二维矩阵Y2，其矩阵形式如下：

[0078] Y2＝(Di,j*),i＝1,2,…,M；j＝1,2,…,N  (7.3)

[0079] 计算矩阵Y2中两两列向量间的相关系数：

[0080]

[0081] 其中， 分别是对矩阵R2标准化后的第m个人和第n个人的中心度得分的均值；

[0082] 将得到的矩阵Y2中两两列向量间的相关系数组成矩阵X2，其矩阵形式如下：

[0083] X2＝(xmn),m＝1,2,…,M；n＝1,2,…,M  (7.5)

[0084] 利用齐次线性方程组求非零解的方法，得出矩阵X2的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，计算矩阵X2的主特征向量P2：

[0085]

[0086] 其中，λmax和eλmax分别是矩阵X2的最大特征值及其对应的单位特征向量；

[0087] 利用得到的矩阵X2的主特征向量P2，计算每个节点的中心度得分H2：

[0088] H2＝R2×P2  (7.7)

[0089] 步骤8：将步骤7.2得到的每个节点的中心度得分H2进行归一化：选取1个最大的中心度得分H2max，将所得的中心度得分H2分别除以H2max，得到H2＇，0＜H2＇≤1；对所得H2＇进行降序排列，前C个H2＇对应的脑区判定为关键脑区。

[0090] 优选地，所述的对每个fMRI数据进行预处理，具体包括以下步骤：

[0091] 对采集到的fMRI数据去除时间点；

[0092] 对去除时间点后的fMRI数据进行时间层校正；

[0093] 对时间层校正后的fMRI数据进行头动校正；

[0094] 对头动校正后的fMRI数据进行空间标准化；

[0095] 对空间标准化后的fMRI数据进行去线性漂移；

[0096] 对去线性漂移后的fMRI数据进行带通滤波。

[0097] 优选地，所述的对采集到的fMRI数据去除时间点，包括：去除时间点的个数为10；所述的对去线性漂移后的fMRI数据进行带通滤波，包括：带通滤波的频率区间为0.01～
0.08Hz。

[0098] 优选地，所述的步骤3，具体包括以下步骤：

[0099] 步骤3.1：

[0100] 选定一个脑模板和Q个稀疏度；

[0101] 根据选定的脑模板将每个fMRI数据划分为N个脑区；

[0102] 定义每个脑区作为每个稀疏度下用于构建脑网络的节点，每个脑网络共有N个节点；

[0103] 步骤3.2：提取每个节点的平均时间序列信号τi；

[0104]

[0105] 其中，τi是指第i个节点的平均时间序列，Gj是指第i个节点中第j个体素的时间序列，ni是第i个节点包含的体素数目；

[0106] 步骤3.3：将步骤3.2中每个节点的平均时间序列信号τi组成矩阵T，T＝(τij)，并计算出两两节点之间的相关系数；

[0107] 节点i和节点j之间的相关系数aij的计算公式：

[0108]

[0109] 其中， 和 分别是i脑区和j脑区时间序列的均值；

[0110] 步骤3.4：将步骤3.3中计算得到的相关系数组成一个N×N的相关系数矩阵A，A＝(aij)；对相关系数矩阵A二值化：对矩阵A中的元素，利用Q个稀疏度确定值为1的元素个数；对一个稀疏度q％，取矩阵A中元素总个数的前q％个数值大的元素，将其值转化为1，其余的元素值转化为0，从而得到一个邻接矩阵B，矩阵B＝(bij)；

[0111] 步骤3.5：根据步骤3.1和3.4中的节点和边，在Q个稀疏度下构建Q个脑网络。

[0112] 优选地，所述的步骤3.1中，选定的脑模板为AAL模板。

[0113] 与现有技术相比，本发明具有的有益效果：

[0114] 1.度量脑区的传统方法是从一个单独的属性来度量脑区的重要性，节点的度和节点效率分别是从局部度量关键脑区，介数中心度是从全局度量关键脑区，这样得到的结果可能会因从不同属性来度量而有所不同，本发明结合了多个属性综合考虑，比如节点的度、节点效率和介数中心度，避免了结果不一致的问题，得到的度量关键脑区的结果更加精准，不但解决了上述问题，而且得到的结果更具代表性。

[0115] 2.本发明综合了多个不同稀疏度，解决了功能脑网络稀疏度阈值的主观选择问题，得到的关键脑区结果更加全面、客观。

[0116] 3.利用主成分分析法，将多种不相关指标利用主成分分析综合在一起来度量关键脑区，同时对多种指标值组成的高维数据逐步降维转变成低维数据，使计算得到简化；并通过综合多种指标、多个稀疏度、多个人得到中心度得分，并根据中心度得分判定关键脑区，该计算方法涵盖面广且精确度高。

附图说明

[0117] 图1为本发明基于主成分分析的关键脑区的度量方法的基本流程示意图。

[0118] 图2为本发明基于主成分分析的关键脑区的度量方法得到的功能性腹泻患者组关键脑区的示意图。

[0119] 图3为本发明由节点的度得到功能性腹泻患者组关键脑区的示意图。

[0120] 图4为本发明由节点效率得到功能性腹泻患者组关键脑区的示意图。

[0121] 图5为本发明由介数中心度得到功能性腹泻患者组关键脑区的示意图。

具体实施方式

[0122] 为了便于理解，对本发明中出现的部分名词作以下解释说明：

[0123] 功能性腹泻患者：英文全称为Functional Diarrhea patients，以下简称为FDI。

[0124] 健康对照组：英文全称为Healthy control group，以下简称为HEA。

[0125] 下面结合附图和具体的实施例对本发明做进一步的解释说明。

[0126] 所采集的fMRI数据是采用德国西门子公司的3T磁共振扫描仪，利用梯度回波平面成像序列对被试者进行采集，其设置的参数如下：TR＝2s；TE＝30ms；翻转角度＝90°；层内分辨率＝3.75×3.75mm2；层厚度＝5mm；矩阵大小＝64；轴状位上扫描层数＝30；总共扫描时间为6分钟，每个被试者包含了180个三维功能像。在扫描过程中，使用乌龙头线圈和泡沫垫子固定被试者头部避免头动；在扫描完成后，被试者会被询问其在扫描期间是否保持放松并清醒的状态，以避免不合格的扫描序列纳入。

[0127] 参考图1～图5所示，实施例一：本发明的一种基于主成分分析的关键脑区的度量方法，包括以下步骤：

[0128] 步骤1：采集M个人静息态下的fMRI数据；

[0129] 步骤2：对每个fMRI数据进行预处理；

[0130] 步骤3：基于预处理后的fMRI数据，构建Q个稀疏度下的脑网络；

[0131] 步骤4：计算每个稀疏度下脑网络的节点的度、节点效率和介数中心度；具体操作包括：

[0132] 步骤4.1：在Q个稀疏度下，计算每个脑网络的节点的度；对于脑网络中任意一个节点i，其节点的度计算公式为：

[0133]

[0134] 其中，wi是对第i个节点的度的加权值，即第i个节点对应脑区的体积大小占整个大脑的比重；bij是指邻接矩阵B中第(i，j)个元素值；N是指每个脑网络中的节点数目；j是指脑网络中第j个节点；所述邻接矩阵B中的元素bij＝0或1，0是指两个节点之间不存在边，1是指两个节点之间存在边，两个节点之间存在的边即构成每个稀疏度下脑网络的边；

[0135] 公式4.1中，wi的计算公式如下：

[0136]

[0137] 其中，ni是第i个节点包含的体素数目；

[0138] 步骤4.2：在Q个稀疏度下，计算每个脑网络的节点效率；

[0139] 对于脑网络中任意一个节点i，其节点效率计算公式为：

[0140]

[0141] 其中，fij是节点i到节点j的最短路径的倒数；

[0142] 步骤4.3：在Q个稀疏度下，计算每个脑网络节点的介数中心度；

[0143] 对于脑网络中任意一个节点i，其介数中心度计算公式为：

[0144]

[0145] 其中，σjk表示从节点j到节点k最短路径的条数，σjk(i)表示在这些最短路径中，并经过节点i的路径条数；

[0146] 步骤5：利用主成分分析将步骤4中得到的节点的度、节点效率和介数中心度进行合并降维得到一种新指标，即每个稀疏度下每个人每个节点的中心度得分H，具体操作如下：

[0147] 步骤5.1：将步骤4中得到的节点的度、节点效率和介数中心度的值重组构成一个二维矩阵R；

[0148] R＝(Df,ijk),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.1)[0149] 其中，Df,ijk表示第i个人第j个节点第k个稀疏度下的第f个属性值；矩阵R的行数为M×N×Q，M为人的数目，N为节点的数目，Q为稀疏度的数目；矩阵R的列分别为节点的度、节点效率和介数中心度的值；

[0150] 步骤5.2：基于主成分分析，对矩阵R进行标准化：

[0151]

[0152] 其中， 是第f个属性值的均值，Sf是第f个属性值的标准差；

[0153] 对矩阵R标准化后，得到一个二维矩阵Y，其矩阵形式如下：

[0154] Y＝(Df,ijk*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.3)[0155] 计算矩阵Y中两两列向量间的相关系数：

[0156]

[0157] 其中， 分别是对矩阵R标准化后的第m个属性值和第n个属性值的均值；

[0158] 将得到的矩阵Y中两两列向量间的相关系数组成矩阵X，其矩阵形式如下：

[0159] X＝(xmn),m＝1,2,3；n＝1,2,3  (5.5)

[0160] 利用齐次线性方程组求非零解的方法，得出矩阵X的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X的主特征向量P：

[0161]

[0162] 其中，λmax和eλmax分别是矩阵X的最大特征值及其对应的单位特征向量；

[0163] 利用得到的矩阵X的主特征向量P，计算每个稀疏度下每个人每个节点的中心度得分H：

[0164] H＝R×P  (5.7)

[0165] 步骤6：利用主成分分析，综合Q个稀疏度下每个人每个节点的中心度得分H，对其降维得到每个人每个节点的中心度得分H1；具体操作如下：

[0166] 步骤6.1：在Q个稀疏度下，将每个人每个节点的中心度得分重组构成一个二维矩阵R1；

[0167] R1＝(Dk,ij),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.1)

[0168] 其中，Dk,ij表示第k个稀疏度下的第i个人第j个节点的中心度得分，矩阵R1的行数为M×N，列数为Q；

[0169] 步骤6.2：基于主成分分析，对矩阵R1进行标准化：

[0170]

[0171] 其中， 是第k个稀疏度下的中心度得分的均值，Sk是第k个稀疏度下的中心度得分的标准差；

[0172] 对矩阵R1标准化后，得到一个二维矩阵Y1，其矩阵形式如下：

[0173] Y1＝(Dk,ij*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.3)

[0174] 计算矩阵Y1中两两列向量间的相关系数：

[0175]

[0176] 其中， 分别是对矩阵R1标准化后的第m个稀疏度下和第n个稀疏度下的中心度得分的均值；

[0177] 将得到的矩阵Y1两两列向量间的相关系数组成一个矩阵X1，其矩阵形式如下：

[0178] X1＝(xmn),m＝1,2,…,Q；n＝1,2,…,Q  (6.5)

[0179] 利用齐次线性方程组求非零解的方法，得出矩阵X1的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X1的主特征向量P1：

[0180]

[0181] 其中，λmax和eλmax分别为矩阵X1的最大特征值及其对应的单位特征向量；

[0182] 利用得到的矩阵X1的主特征向量P1，计算每个人每个节点的中心度得分H1：

[0183] H1＝R1×P1  (6.7)

[0184] 步骤7：利用主成分分析，综合M个人每个节点的中心度得分H1，并对其降维得到每个节点的中心度得分H2；具体操作如下：

[0185] 步骤7.1：将M个人的每个节点的中心度得分H1重组构成一个二维矩阵R2；

[0186] R2＝(Di,j),i＝1,2,…,M；j＝1,2,…,N  (7.1)

[0187] 其中，Di,j表示第i个人的第j个节点的中心度得分，矩阵R2的行数为N，列数为M；

[0188] 步骤7.2：基于主成分分析，对矩阵R2进行标准化：

[0189]

[0190] 其中，是第i个人的中心度得分的均值，Si是第i个人的中心度得分的标准差；

[0191] 对矩阵R2标准化后，得到一个二维矩阵Y2，其矩阵形式如下：

[0192] Y2＝(Di,j*),i＝1,2,…,M；j＝1,2,…,N  (7.3)

[0193] 计算矩阵Y2中两两列向量间的相关系数：

[0194]

[0195] 其中， 分别是对矩阵R2标准化后的第m个人和第n个人的中心度得分的均值；将得到的矩阵Y2中两两列向量间的相关系数组成矩阵X2，其矩阵形式如下：

[0196] X2＝(xmn),m＝1,2,…,M；n＝1,2,…,M  (7.5)

[0197] 利用齐次线性方程组求非零解的方法，得出矩阵X2的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，计算矩阵X2的主特征向量P2：

[0198]

[0199] 其中，λmax和eλmax分别是矩阵X2的最大特征值及其对应的单位特征向量；

[0200] 利用得到的矩阵X2的主特征向量P2，计算每个节点的中心度得分H2：

[0201] H2＝R2×P2  (7.7)

[0202] 步骤8：将步骤7.2得到的每个节点的中心度得分H2进行归一化：选取1个最大的中心度得分H2max，将所得的中心度得分H2分别除以H2max，得到H2＇，0＜H2＇≤1；对所得H2＇进行降序排列，前C个H2＇对应的脑区判定为关键脑区。

[0203] 作为一种可实施的方式，本实施例中的人数M设为20，稀疏度个数Q设为30，脑区个数或脑网络中的节点数N设为90；前C个H2＇对应的脑区判定为关键脑区，其中，C设为10。

[0204] 实施例二：本发明的一种基于主成分分析的关键脑区的度量方法，包括以下步骤：

[0205] 步骤1：采集M名功能性腹泻患者静息态下的fMRI数据；

[0206] 步骤2：对每个fMRI数据进行预处理；具体操作如下：

[0207] 步骤2.1：对采集到的fMRI数据去除前10个时间点；

[0208] 步骤2.2：对去除时间点后的fMRI数据进行时间层校正；

[0209] 步骤2.3：对时间层校正后的fMRI数据进行头动校正；

[0210] 步骤2.4：对头动校正后的fMRI数据进行空间标准化；

[0211] 步骤2.5：对空间标准化后的fMRI数据进行去线性漂移；

[0212] 步骤2.6：对去线性漂移后的fMRI数据进行带通滤波，带通滤波器的频率区间为0.01～0.08Hz；

[0213] 步骤3：基于预处理后的fMRI数据，构建Q个稀疏度下的脑网络；具体操作如下：

[0214] 步骤3.1：

[0215] 选定一个脑模板和Q个稀疏度；

[0216] 根据选定的脑模板将每个fMRI数据划分为N个脑区；

[0217] 定义每个脑区作为每个稀疏度下用于构建脑网络的节点，每个脑网络共有N个节点；

[0218] 步骤3.2：提取每个节点的平均时间序列信号τi；

[0219]

[0220] 其中，τi是指第i个节点的平均时间序列，Gj是指第i个节点中第j个体素的时间序列，ni是第i个节点包含的体素数目；

[0221] 步骤3.3：将步骤3.2中每个节点的平均时间序列信号τi组成矩阵T，T＝(τij)，并计算出两两节点之间的相关系数；

[0222] 节点i和节点j之间的相关系数aij的计算公式：

[0223]

[0224] 其中， 和 分别是i脑区和j脑区时间序列的均值；

[0225] 步骤3.4：将步骤3.3中计算得到的相关系数组成一个N×N的相关系数矩阵A，A＝(aij)；对相关系数矩阵A二值化：对矩阵A中的元素，利用Q个稀疏度确定值为1的元素个数；例如：稀疏度为1％时，取，即将矩阵A中前81个数值大的元素的值转化为1，其余的元素转化为0，从而得到一个邻接矩阵B，矩阵B＝(bij)；

[0226] 步骤3.5：根据步骤3.1和3.4中的节点和边，在Q个稀疏度下构建Q个脑网络；

[0227] 步骤4：在Q个稀疏度下，计算每个脑网络的节点的度、节点效率和介数中心度；具体操作包括：

[0228] 步骤4.1：在Q个稀疏度下，计算每个脑网络的节点的度；对于脑网络中任意一个节点i，其节点的度计算公式为：

[0229]

[0230] 其中，wi是对第i个节点的度的加权值，即第i个节点对应脑区的体积大小占整个大脑的比重；bij是指邻接矩阵B中第(i，j)个元素值；N是指每个脑网络中的节点数目；j是指脑网络中第j个节点；所述邻接矩阵B中的元素bij＝0或1，0是指两个节点之间不存在边，1是指两个节点之间存在边，两个节点之间存在的边即构成每个稀疏度下脑网络的边；

[0231] 公式4.1中，wi的计算公式如下：

[0232]

[0233] 其中，ni是第i个节点包含的体素数目；

[0234] 步骤4.2：在Q个稀疏度下，计算每个脑网络的节点效率；

[0235] 对于脑网络中任意一个节点i，其节点效率计算公式为：

[0236]

[0237] 其中，fij是节点i到节点j的最短路径的倒数；

[0238] 步骤4.3：在Q个稀疏度下，计算每个脑网络节点的介数中心度；

[0239] 对于脑网络中任意一个节点i，其介数中心度计算公式为：

[0240]

[0241] 其中，σjk表示从节点j到节点k最短路径的条数，σjk(i)表示在这些最短路径中，并经过节点i的路径条数；

[0242] 步骤5：利用主成分分析将步骤4中得到的节点的度、节点效率和介数中心度进行合并降维得到一种新指标，即每个稀疏度下每个人每个节点的中心度得分H，具体操作如下：

[0243] 步骤5.1：将步骤4中得到的节点的度、节点效率和介数中心度的值重组构成一个二维矩阵R；

[0244] R＝(Df,ijk),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.1)[0245] 其中，Df,ijk表示第i个人第j个节点第k个稀疏度下的第f个属性值；矩阵R的行数为M×N×Q，M为人的数目，N为节点的数目，Q为稀疏度的数目；矩阵R的列分别为节点的度、节点效率和介数中心度三个指标的值；

[0246] 步骤5.2：基于主成分分析，对矩阵R进行标准化：

[0247]

[0248] 其中， 是第f个属性值的均值，Sf是第f个属性值的标准差；

[0249] 对矩阵R标准化后，得到一个二维矩阵Y，其矩阵形式如下：

[0250] Y＝(Df,ijk*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q；f＝1,2,3  (5.3)[0251] 计算矩阵Y中两两列向量间的相关系数：

[0252]

[0253] 其中， 分别是对矩阵R标准化后的第m个属性值和第n个属性值的均值；

[0254] 将得到的矩阵Y中两两列向量间的相关系数组成矩阵X，其矩阵形式如下：

[0255] X＝(xmn),m＝1,2,3；n＝1,2,3  (5.5)

[0256] 利用齐次线性方程组求非零解的方法，得出矩阵的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应的单位特征向量eλmax，并计算矩阵的主特征向量P：

[0257]

[0258] 其中，λmax和eλmax分别是矩阵X的最大特征值及其对应的单位特征向量；

[0259] 利用得到的矩阵X的主特征向量P，计算每个稀疏度下每个人每个节点的中心度得分H：

[0260] H＝R×P  (5.7)

[0261] 步骤6：利用主成分分析，综合Q个稀疏度下每个人每个节点的中心度得分H，对其降维得到每个人每个节点的中心度得分H1；具体操作如下：

[0262] 步骤6.1：在Q个稀疏度下，将每个人每个节点的中心度得分重组构成一个二维矩阵R1；

[0263] R1＝(Dk,ij),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.1)

[0264] 其中，Dk,ij表示第k个稀疏度下的第i个人第j个节点的中心度得分，矩阵R1的行数为M×N，列数为Q；

[0265] 步骤6.2：基于主成分分析，对矩阵R1进行标准化：

[0266]

[0267] 其中， 是第k个稀疏度下的中心度得分的均值，Sk是第k个稀疏度下的中心度得分的标准差；

[0268] 对矩阵R1标准化后，得到一个二维矩阵Y1，其矩阵形式如下：

[0269] Y1＝(Dk,ij*),i＝1,2,…,M；j＝1,2,…,N；k＝1,2,…,Q  (6.3)

[0270] 计算矩阵Y1中两两列向量间的相关系数：

[0271]

[0272] 其中， 分别是对矩阵R1标准化后的第m个稀疏度下和第n个稀疏度下的中心度得分的均值；

[0273] 将得到的矩阵Y1两两列向量间的相关系数组成一个矩阵X1，其矩阵形式如下：

[0274] X1＝(xmn),m＝1,2,…,Q；n＝1,2,…,Q  (6.5)

[0275] 利用齐次线性方程组求非零解的方法，得出矩阵X1的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，并计算矩阵X1的主特征向量P1：

[0276]

[0277] 其中，λmax和eλmax分别为矩阵X1的最大特征值及其对应的单位特征向量；

[0278] 利用得到的矩阵X1的主特征向量P1，计算每个人每个节点的中心度得分H1：

[0279] H1＝R1×P1  (6.7)

[0280] 步骤7：利用主成分分析，综合M个人每个节点的中心度得分H1，并对其降维得到每个节点的中心度得分H2；具体操作如下：

[0281] 步骤7.1：将M个人的每个节点的中心度得分H1重组构成一个二维矩阵R2；

[0282] R2＝(Di,j),i＝1,2,…,M；j＝1,2,…,N  (7.1)

[0283] 其中，Di,j表示第i个人下的第j个节点的中心度得分，矩阵R2的行数为Q，列数为M；

[0284] 步骤7.2：基于主成分分析，对矩阵R2进行标准化：

[0285]

[0286] 其中， 是第i个人的中心度得分的均值，Si是第i个人的中心度得分的标准差；

[0287] 对矩阵R2标准化后，得到一个二维矩阵Y2，其矩阵形式如下：

[0288] Y2＝(Di,j*),i＝1,2,…,M；j＝1,2,…,N  (7.3)

[0289] 计算矩阵Y2中两两列向量间的相关系数：

[0290]

[0291] 其中， 分别是对矩阵R2标准化后的第m个人和第n个人的中心度得分的均值；

[0292] 将得到的矩阵Y2中两两列向量间的相关系数组成矩阵X2，其矩阵形式如下：

[0293] X2＝(xmn),m＝1,2,…,M；n＝1,2,…,M  (7.5)

[0294] 利用齐次线性方程组求非零解的方法，得出矩阵的所有特征值λi及其对应的单位特征向量ei，在得到的特征值λi及其对应的单位特征向量ei中，选取最大特征值λmax及其对应单位特征向量eλmax，计算矩阵X2的主特征向量P2：

[0295]

[0296] 其中，λmax和eλmax分别是矩阵X2的最大特征值及其对应的单位特征向量；

[0297] 利用得到的矩阵X2的主特征向量P2，计算每个节点的中心度得分H2：

[0298] H2＝R2×P2  (7.7)

[0299] 步骤8：将步骤7.2得到的每个节点的中心度得分X2进行归一化：选取1个最大的中心度得分H2max，将所得的中心度得分H2分别除以H2max，得到H2，0＜H2＇≤1；对所得H2＇进行降序排列，前C个H2＇对应的脑区判定为关键脑区。

[0300] 作为一种可实施的方式，本实施例选取了一组功能性腹泻患者共20个人作为被试对象，即M＝20；在允许的稀疏度0～50％范围内，从1％～30％平均选定30个稀疏度，即Q＝30；脑区个数或脑网络中的节点数N设为90；前C个H2＇对应的脑区判定为关键脑区，C设为
10。

[0301] 作为一种可实施的方式，步骤6和步骤7的顺序不受限制，除上述提出的步骤之外，也可以先进行步骤7，即先综合20个人的每个稀疏度下每个节点的中心度得分H，并对其降维得到每个稀疏度下每个节点的中心度得分；再进行步骤6，综合30个稀疏度下的每个节点的中心度得分，对其降维得到每个节点的中心度得分，所得结果并不会影响对关键脑区的判定。

[0302] 如表1所示，本发明对20名功能性腹泻患者不同脑区关键性的度量结果，中心度得分1为从步骤6到步骤7计算的中心度得分，中心度得分2为先进行步骤7再进行步骤6得到的中心度得分；其中节点的度、节点效率和介数中心度，是由每项指标分别计算出的结果；表1中第一列为90个脑区的中文名称，其中L代表大脑中的左脑，R为大脑中的右脑，例如：中央前回L，代表左脑中的中央前回脑区；中央前回R代表右脑中的中央前回脑区。

[0303] 表1不同度量方法得到的功能性腹泻患者90个脑区的关键性得分

[0304]

[0305]

[0306] 实施例三：本实施例选取了一组健康人作为被试对象，人数为20个，作为功能性腹泻患者的健康对照组，并得出脑区关键性的度量结果；本发明在实施例三中，只是换取一组被试对象，操作步骤与实施例二完全相同，此处不再赘述；最终得出的90个脑区的中心度得分值如表2所示。

[0307] 表2不同度量方法得到的健康对照人群90个脑区的中心度得分

[0308]

[0309]

[0310]

[0311] 表3不同度量方法得到的前10个关键脑区

[0312]

[0313]

[0314] 根据实施例二和实施例三中的表1和表2得到的数据，将功能性腹泻患者与健康对照组的不同得分结果进行统计，将得分值最大的前10个脑区判定为关键脑区，最后得到几种方法的关键脑区列表，其中打勾的脑区为关键脑区，结果如表3所示。由表3可以看出，无论是健康对照组还是功能性腹泻患者组，通过中心度得分1确定的关键脑区与中心度得分2确定的关键脑区都是一致的，这说明本发明在计算健康组和患者组的关键脑区的过程中，无论是先综合稀疏度再综合所有人，还是先综合所有人再综合稀疏度，对其得到的结果都是不影响的。传统上是通过三种基本方法：节点的度、节点效率或介数中心度，从三个角度来度量关键脑区的，由表3可以看出，无论是健康对照组还是功能性腹泻患者组，以上三种基本方法所得到的结果并非完全一致，无法确定哪种方法的结果更具有合理性。本发明对三种基本方法进行了加权处理，即在三种基本方法基础上考虑了各个脑区在大脑中所占的比重，且综合三者特征，得到的结果更为全面、精确、合理。在确定患者组的关键脑区时，本发明得到的关键脑区与介数中心度方法得到的关键脑区完全一致，但这只是一种巧合，并不能说明三种基本方法中介数中心度方法更好。通过健康对照组与功能性腹泻患者组的关键脑区的对比，发现两组人之间的关键脑区还是存在较大的差异，其差异可能与功能性腹泻疾病的病理生理学有关，甚至在该疾病的不同阶段，患者的关键脑区也可能会发生显著变化，因此，检测患者的关键脑区具有重要意义。本发明基于主成分分析的关键脑区的度量方法，对于关键脑区的判定是合理且可行的，对其进行的研究将对脑疾病相关的临床医学具有重要的指导意义。

[0315] 以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。