[0001] 本发明属于大脑成像技术领域，具体涉及一种脑动态功能模式稳定性计算方法。

[0002] 稳定性是大脑功能活动的重要属性之一。大脑通过广泛的神经活动和连接，对信息产生稳定的表征，是其有效运作的关键。虽然近年来有许多研究都考察了大脑功能组织模式的动态特征，但仍缺乏研究描绘其功能组织模式的稳定性，特别是刻画出大脑每个区域的功能稳定性。究其原因，目前没有合适的算法和工具考察体素水平的功能稳定性指标。

[0003] 现有研究采用脑功能成像技术考察了静息态下脑功能组织模式的灵活性和变异度，虽然这两个指标与稳定性有关联，但不能作为稳定性的替代。并且，这些研究依据AAL(Automated Anatomical Labeling)模板划分大脑，计算基于分区的灵活性指标。AAL模板是基于大脑结构信息创建的，无法准确反应大脑的功能组织架构。此外，还有研究关注脑功能不同状态的转换度和稳定性，但是这些属性是对脑整体的和宏观的反映，无法测量单个脑区乃至体素的功能稳定性。

[0004] 中国专利文献CN 107832656中公开了一种基于动态功能脑网络的大脑功能状态信息处理方法，通过采集全脑的各个脑区内各体素的时间信号，求解各脑区的平均时间信号；构建全脑动态功能连接网络，进行动态功能连接网络聚类分析，进行脑网络功能连接微状态特征量化分析。其是通过对脑网络功能连接微状态特征量如微状态停留时间、微状态转换次数和微状态稳定性进行分析实现对大脑功能状态的判别。通过上述专利文献进行分析，其通过求解各脑区的平均时间信号还是基于某个模板的脑分区，如AAL模板，有116个分区，分区数量取决于所选取的模板，通过基于脑分区的分析往往不精确，因为取各脑区的平均时间信号会丢失一些有意义的信号，且所采用的模板划分的分区不一定能够反映大脑的实际功能分区。

[0038] 下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0039] 本发明用于分析脑功能成像数据，测量被试在一个连续状态下(如静息态，或观看一段影片时)的脑功能稳定性。前人研究发现脑区间的连接随着时间发生动态变化，目前以全脑功能连接模式在一段时间内前后一致性作为脑功能稳定性的度量。本发明中所提供的算法是对全脑每个灰质体素进行稳定性逐一计算。对于单个体素，以该体素与其他脑区(体素)动态功能连接作为度量其功能稳定性的特征。

[0040] 如图2所示，具体计算方法包括如下步骤：

[0041] 【S1】.获取被试的功能磁共振成像(fMRI)数据和结构磁共振成像数据，并采用基于MATLAB的工具包SPM12进行脑功能数据预处理；

[0042] 【S2】.通过计算获得被试脑灰质体积图像的平均图像，且提取脑灰质体积为大于0.2的灰质体素作为后续计算的掩膜(mask)；

[0043] 【S3】.采用滑动时间窗方法计算每一灰质体素在每个时间窗下与全脑其余灰质体素动态功能连接的二维矩阵；

[0044] 【S4】.以时间窗为评分者对所得每一个灰质体素的二维矩阵进行计算，得到二维矩阵的肯德尔和谐系数作为该灰质体素的功能稳定性值；

[0045] 【S5】.对被试的全脑灰质体素进行逐一计算后得到所有灰质体素的功能稳定性值；

[0046] 【S6】.对所有灰质体素的功能稳定性值进行Z-score标准化，形成并输出被试的脑动态功能稳定性三维脑图。

[0047] 具体计算步骤细节如下：

[0048] (1)首先对获取到的功能性磁共振成像(fMRI，functional magnetic resonance imaging)数据(扫描时间不少于5分钟)进行预处理。预处理可使用常见的工具按常规流程进行，常用工具如SPM、FSL和DPARSF等。本发明所采用优选的预处理如下：删除前10秒的功能性磁共振成像数据，采用最小二乘法进行头动校正，再通过线性插值法进行时间层校正；采用线性模型把白质、脑脊液和头动的无关信号回归去除；采用“New Segment”算法把被试结构像划分多种脑组织；采用“DARTEL”算法把被试结构像标准化到MNI(Montreal Neurological Institute)空间，通过配准并结合结构像的变形场把功能像空间标准化到MNI空间，把功能像重采样为3×3×3mm的分辨率空间，滤波保留0.01-0.08频段的信号，以
6mm FWHM(全宽半最大值)的高斯核进行空间平滑。这里的“New Segment”、“DARTEL”等算法为预处理中优选的通用算法，对此本发明不再具体赘述。

[0049] (2)在上述(1)对结构像预处理的基础上进行Jocobian scaled变换，获得每个被试的脑灰质体积图像，并重采样到3×3×3mm的空间，与功能像保持一致。以组间平均灰质体积大于0.2为阈限提取mask掩膜，作为后续分析的掩膜，下面所提及的灰质体素都是该灰质mask掩膜内的体素。灰质体积0.2是通用的用以区分灰质和其他无关脑组织的阈限。mask掩膜只保留了灰质体素(在3×3×3mm空间下约保留44000个体素)，使得稳定性的计算不受到白质和脑脊液等与神经活动无关信号的影响。

[0050] (3)对于预处理后的功能性磁共振成像fMRI数据，采用滑动时间窗方法计算动态功能连接。对于某一体素，计算其在每个时间窗与全脑其他体素(约44000个)的血氧水平依赖(BOLD)信号的皮尔逊相关系数。每个时间窗得到一个功能连接(即皮尔逊相关系数)的一维向量，反映了该体素在当前时间窗的功能连接模式。所有时间窗的功能连接向量在时间维度上连起来构成的二维矩阵，其中一个维度是灰质体素，另一个维度是时间窗，表示该体素的功能连接模式随着时间变化的情况。滑动时间窗方法需要设置三个参数，包括时间窗宽度(如1分钟)、滑动步长(如4秒，2个TR)和时间窗类似(如海明窗，hamming window)。例如8分钟扫描时间、TR为2秒的fMRI数据，按照上述参数设置，时间窗数量为101个。在上述条件下，度量某一体素功能稳定性的特征二维矩阵大小为44000×101。

[0051] (4)根据下列公式，以时间窗为评分者计算功能连接(3)计算得到的二维矩阵的肯德尔和谐系数：

[0052]

[0053]

[0054] W代表某一体素的肯德尔和谐系数；K代表时间窗的数量；N代表该体素的功能连接数量(即全部灰质体素数量减一)；对于每个时间窗，把所有功能连接根据其连接强度值做等级排序，并根据等级进行赋值，连接强度值越高，等级越高(如连接强度值最高的功能连接赋值为44000)；Rn是第n条功能连接的排序等级在所有时间窗的总和。W的取值范围从0到1，用于衡量该体素的全脑功能连接模式在时间上的一致性程度，作为该体素的功能稳定性值。W的值越高，则反映了该体素具有较高的功能稳定性。步骤(3)(4)可参见图1。图1所示大脑某灰质体素在每个时间窗与灰质mask内其他体素的动态功能连接，构成定义其功能稳定性的特征。采用滑动窗方法计算动态功能连接；对于滑动窗方法的参数，图1采用矩形窗，时间窗长度64秒，滑动步长4秒。基于该体素的动态功能连接二维矩阵(矩阵大小是灰质体素数量×时间窗数量)，以时间窗为评分者，计算肯德尔和谐系数，作为量化该体素功能稳定性的数值。

[0055] (5)对单个被试的全脑灰质体素进行逐一计算后，为了降低被试间和条件间的整体差异，对所有灰质体素的稳定性数值(即W值)进行z-score标准化，其是基于原始数据的均值(mean)和标准差(standard deviation)进行数据的标准化，最后输出该被试nii格式的功能稳定性三维脑图，见图3所示多层水平面脑图，其为一名被试的全脑灰质体素的功能稳定性，颜色由浅到深表示稳定性由低到高；图中深色部分主要分布在后部扣带皮层、角回和脑岛前部等脑区，说明它们相对于其他脑区具有更高水平的功能稳定性。

[0056] 通过将本发明所提供的算法用于静息态数据，首次描绘出全脑灰质的内在脑功能稳定性分布，发现了相对于单模态的感觉运动皮层，参与高级认知功能加工的联合皮层具有更高的功能稳定性。功能稳定性作为大脑的关键属性，其算法的开发为研究人脑功能正常的运作机制及异常的生物标志提供新的角度。

[0057] 另外，考虑到有些研究人员的计算资源有限，本发明还对算法进行拓展，通过调整具体算法减少稳定性度量的特征数量，大幅减少计算资源消耗。具体方法如下：

[0058] 选取合适的脑功能分区图谱，以单个体素与该图谱的几百个功能分区的动态功能连接，形成调整后的体素-功能分区方法，替代其与几万个体素的动态功能连接(即上述的体素-体素方法)，作为其功能稳定性的度量特征。由于基于结构信息开发的图谱(如Anatomical Automatic Labeling)无法准确反映大脑的功能分区，本发明采用Craddock等(2012)开发的功能分区图谱，该图谱包括200个分区和400个分区两种。对于调整后的算法，上述(3)调整为计算某一灰质体素在每个时间窗与各个分区(200个或400个)的平均BOLD信号的皮尔逊相关系数，所有时间窗的功能连接向量连起来构成的二维数据大小仅为200×101(或400×101)，与44000×101相比，大幅减少运算资源需求。

[0059] 通过分析发现，对于全脑的所有灰质体素，采用体素-功能分区方法计算得到的稳定性数值，与采用体素-体素方法计算得到的稳定性数值，存在显著高相关，相关值均值大于0.92，即保留了大部分功能稳定性信息。调整后的体素-功能分区方法适用于计算资源有限的条件下或探索性分析。

[0060] 由于滑动时间窗的参数设置没有定论和标准，我们还检测本发明算法是否对滑动时间窗参数具有鲁棒性，即不同时间窗参数对功能稳定性计算是否有显著影响。我们逐个改变时间窗参数，包括把时间窗长度减半(30秒)或加倍(2分钟)、把滑动步长减半(2秒)、采用矩形窗等。通过分析真实的静息态fMRI数据，发现时间窗参数的改变对功能稳定性的分布影响很小。

[0061] 实施例：

[0062] 通过磁共振仪器采集了两组被试(如抑郁症患者与正常对照组)的静息态fMRI数据和结构像数据，两组被试各40名，在年龄、性别和教育程度等人口学变量上匹配。研究目的是分析比较抑郁症患者与正常对照组的脑功能稳定性差异。fMRI数据是当个体处于静息态下采集，扫描时间长度6分钟，TR为1秒。

[0063] (1)首先采用SPM12或DPABI等基于MATLAB的工具包对fMRI数据进行常规的静息态脑功能数据预处理。排除功能图像质量不好以及头动过大的被试。计算两组剩余被试的灰质体积图像的平均图像，以0.2为阈限，提取灰质体积大于0.2的体素作为后续计算的mask掩膜。

[0064] (2)设置滑动时间窗参数：时间窗长度1分钟、滑动步长4秒(即4个TR)、海明窗。对于单个被试，采用滑动时间窗方法，计算某一灰质体素与其他灰质体素的动态功能连接，得到m×n矩阵，m为功能连接数量(约44000个)，n为时间窗数量(71个)。并对该矩阵以时间窗为评分者计算肯德尔和谐系数，得到该灰质体素的功能稳定性值。

[0065] (3)对该被试在mask内的灰质体素进行逐一计算后，对所有灰质体素的功能稳定性值进行Z标准化，其他mask外的非灰质体素置0，而后输出一个三维图像的nii文件，代表该被试的功能稳定性脑图。

[0066] (4)对两组的每个被试都进行上述计算，得到各个被试的功能稳定性nii文件。进一步的可将这些获取的nii文件用于统计分析，如抑郁症组与正常对照组的双样本T检验。

[0067] 显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。