脑功能连接频率范围的定量判定方法

摘要

本发明提供了一种自动判别功能连接精确频率范围的方法，包括如下步骤：S1、获取脑功能成像数据；S2、对获取的所述脑功能成像数据进行预处理；S3、根据先验信息生成感兴趣脑功能系统的空间模板；S4、选取所述感兴趣脑功能系统内的某一区域作为种子点，采用信号同步性频谱分析方法计算每个频率点上所有体素和种子点的相关性；S5、逐频率地对所述计算得到的所有体素和种子点的相关性结果进行空间加权求和，得到所述感兴趣脑功能系统的功能连接频率曲线；S6、对步骤S5中功能连接频率曲线逐频率进行显著性判断。根据本发明的方法不仅能够基于数据计算出功能连接-频率分布曲线，并基于该曲线自动判别出精确的功能连接频率范围，还能够克服脑功能成像数据中全局噪声对功能连接频率分析的干扰。

权利要求

1.一种自动判别功能连接精确频率范围的方法，包括如下步骤：S1、获取脑功能成像数据；

S2、对获取的所述脑功能成像数据进行预处理；

S3、根据先验信息生成感兴趣脑功能系统的空间模板；

S4、选取所述感兴趣脑功能系统内的某一区域作为种子点，采用信号同步性频谱分析方法计算每个频率点上所有体素和种子点的相关性；

S5、逐频率地对所述计算得到的所有体素和种子点的相关性结果进行空间加权求和，得到所述感兴趣脑功能系统的功能连接频率曲线；

S6、对步骤S5中功能连接频率曲线逐频率进行显著性判断；

其中，所述根据先验信息生成感兴趣脑功能系统的空间模板为：根据先验信息确定感兴趣功能区域的空间体素位置，生成与脑功能数据匹配的感兴趣功能系统的空间模板。

2.根据权利要求1所述的方法，其特征在于：所述步骤4中计算每个频率点上所有体素和种子点的相关性包括：计算所述区域中所有体素和种子点的同步性频谱分布曲线，组成空间-频率连接矩阵C，逐频率生成同步性空间分布图。

3.根据权利要求1或2所述的方法，其特征在于：所述步骤4中的信号同步性频谱分析方法是coherence相干分析方法。

4.根据权利要求1或2所述的方法，其特征在于：所述步骤4中的信号同步性频谱分析方法是Pearson相关系数谱分解的方法。

5.根据权利要求3所述的方法，其特征在于：所述coherence相干分析方法具体为：其中cxs(λ)为在频率点λ上，时间序列x与s之间的相干系数，其值为0～1，0表示在该频率时间序列之间无线性关系，1表示两个时间序列之间存在理想的线性关系，Gxs表示时间序列x与s之间的互功率谱，Gxx和Gss分别表示时间序列x和s的自功率：式中，Fx和Fs分别表示时间序列x和s傅里叶频谱分析得到的频谱，*表示共轭。

6.根据权利要求4所述的方法，其特征在于：所述Pearson相关系数谱分解的方法包括：将时间序列的点数记为N，Pearson相关系数计算如下：其中，根据离散傅里叶反变换法则：

随后，将cxs分解到频率域上得到：

7.根据权利要求1、2、5或6任意之一所述的方法，其特征在于：所述步骤S5中逐频率地对所述计算得到的所有体素和种子点的相关性结果进行空间加权求和进一步包括：S51、利用所述步骤S3中生成的感兴趣脑功能系统的空间模板m，以及每个频率点λi同步性空间分布图cs(λi)，计算出该频率点上每个体素x对应的空间加权系数wx(λi)其中，M为体素总个数， 为m的空间平均数， 为cs的空间平均数；

S52、利用wx，逐频率地计算同步性空间分布图cs的空间加权和系数，即：SW

所得的c 即为感兴趣功能系统的功能连接频率曲线。

8.根据权利要求2、5或6任意之一所述的方法，其特征在于：所述步骤S6包括：S61、将种子点时间序列的功率谱固定，随机变化它的相位频谱，通过傅里叶反变换得到相位随机化后的种子点时间序列，其他体素的时间序列不做任何变化；

S62、按照所述步骤S4的方法计算所有体素与相位随机化后的种子点时间序列的空间-频率同步化矩阵，以及对应的功能连接频率曲线；

S63、将随机化过程重复若干次，利用所有随机得到的功能连接频率曲线上的值，生成频域功能连接强度的零假设分布；

S64、逐频率点地计算感兴趣功能系统的功能连接频率曲线值的显著性水平p，设定阈值p0，其中，0＜p0＜1，在p＜p0时，显著的频率即为功能连接的主要频率。

说明书

脑功能连接频率范围的定量判定方法

技术领域

[0001] 本发明涉及神经影像计算领域，具体涉及一种基于脑功能成像数据的针对脑功能连接的频率范围的定量判定方法。

背景技术

[0002] 功能连接（Functional connectivity）是神经科学中一个越来越受到关注的研究领域。它通过度量空间位置相互分离的大脑活动的同步性，来研究脑区间信息的相互传递及功能的交互联系。现有研究发现，由功能连接构成的脑功能网络与脑解剖网络相似，并且功能连接的具体属性（主要是强度）能够反映疾病引起的脑功能受损，体现发育、学习和疾病康复过程中脑功能整合模式的可塑性变化，预测被试行为的差异。这表明，功能连接是脑部神经元活动的相互作用的体现，能够反映大脑内在的功能构架，具有非常重要的功能意义。

[0003] 除了强度属性，近年来有证据表明，功能连接的频率属性也可能具有一定功能意义。虽然神经元之间的交互活动频率非常宽泛（0.01-500Hz），但是根据神经-血管耦合关系（neurovascular coupling），由血流信号体现出的功能连接往往集中于某一特定的低频范围。并且该频率范围可能与所属脑区的功能有关。例如，初级感觉运动区内部的功能连接集中于0.1Hz以下的频段，而边缘系统（如双侧杏仁核之间）的功能连接分布在相对更宽泛的频率范围（0-0.25Hz）。现在关于功能连接频率属性的研究还处于起步阶段，数量不足十篇，其具体的功能意义还需要进一步探索。

[0004] 现有关于功能连接频率范围的判定方法主要有两种。一是根据傅里叶变换法则，将时域计算的功能连接强度分解到频率域，再通过肉眼主观判断其能量集中的粗略频率范围；二是研究者先将频率段划分为若干小段，再在每一频率小段中计算功能连接强度，最后比较判断其强度相对较高的频率大体范围。这两种方法都依赖于研究者的主观判断，客观性差；并且 判别结果频率分辨率粗糙，可能无法区分出不同功能系统功能连接频率的细致差别。更重要的是，现有两种方法都无法有效避免脑功能成像数据（如功能核磁共振成像fMRI、功能近红外成像fNIRS）中的全局噪声（如呼吸、心跳等生理噪声、头动噪声等空间上大面积存在的噪声）干扰问题。特别是功能近红外成像，它作为一种新兴的非侵入的脑成像技术，近三年间已迅速受到功能连接研究者的关注，被成功应用于健康人、中风患者、及婴儿的功能整合研究。但是由于它是基于从头皮表面出发经过脑组织再回到头皮表面的光子的吸收散射特性成像，数据中包含有大量的头皮部分的血流信号。已有研究发现，该噪声足以扰乱甚至破坏研究者对脑组织活动信号的分析。另外在实际试验过程中，由于光极帽固定不理想导致的成像光极与头皮的偶尔相对滑动也会使数据出现大面积的毛刺或台阶式的头动噪声，从而干扰数据分析的进行。由于头皮血流中的生理噪声和头动噪声的频率分布都很宽泛且复杂，与功能连接很有可能会有频率混叠，所以这些噪声在以上两种频率判别方法中的残留，将很有可能影响功能连接频率范围判断的准确性。现在还没有可靠有效的数据预处理的方法能够解决这个问题。

发明内容

[0005] 为了克服上述已有技术均依赖于研究者的主观判断，客观性差；并且判别结果频率分辨率粗糙的缺陷；同时，更重要的是，为了克服现有两种方法都无法有效避免脑功能成像数据（如功能核磁共振成像fMRI、功能近红外成像fNIRS）中的全局噪声（如呼吸、心跳等生理噪声、头动噪声等空间上大面积存在的噪声）干扰问题。本发明提出了一种自动判别功能连接精确频率范围的方法，包括如下步骤：S1、获取脑功能成像数据；S2、对获取的所述脑功能成像数据进行预处理；S3、根据先验信息生成感兴趣脑功能系统的空间模板；S4、选取所述感兴趣脑功能系统内的某一区域作为种子点，采用信号同步性频谱分析方法计算每个频率点上所有体素和种子点的相关性；S5、逐频率地对所述计算得到的所有体素和种子点的相关性结果进行空间加权求和，得到所述感兴趣脑功能系统的功能连接频率曲线；S6、对步骤S5中功能连接频率曲线逐频率进行显著性判断。

[0006] 本发明利用相干分析等信号同步性频谱分析方法，分析每个频率点上所有体素和种子点的相关性；并将空间信息引入频率分析中，通过空间加权求和的方法消除全局噪声SW对同步性频率分析的影响；基于逐频率点对c 进行统计检验，自动判别出功能连接的精确频率范围。该方法计算过程快速，在普通微机上即可完成。这种方法主要用于各种相对稳定状态（如静息状态）的脑功能连接频率范围探索性研究，也适用于其他非线性状态（如药物应用对脑功能的影响）的脑功能连接频率范围变化的分析。

[0007] 根据本发明的自动判别功能连接精确频率范围的方法，它不仅能够基于数据计算出功能连接-频率分布曲线，并基于该曲线自动判别出精确的功能连接频率范围，还能够克服脑功能成像数据中全局噪声对功能连接频率分析的干扰。

附图说明

[0008] 图1是本发明实现过程的流程图；

[0009] 图2是空间-频率连接矩阵C的示意图；

[0010] 图3是具体实例中近红外光谱成像具体光极摆放位置示意图；

[0011] 图4是具体实例中不同频率的同步性空间分布图；

[0012] 图5a是根据本发明的方法所得的感觉运动系统功能连接频率曲线图；

[0013] 图5b是传统方法所得的感觉运动系统功能连接频率曲线图；

[0014] 图6a是模拟脑自发波动的频谱分布（左）和时间曲线示例（右）；

[0015] 图6b是模拟脑自发波动的空间分布图示例；

[0016] 图6c是模拟系统生理噪声的频谱分布（左）和时间曲线示例（右）；

[0017] 图6d是模拟系统生理噪声的振幅（上）和相位（下）的空间分布图示例；

[0018] 图6e是模拟头动噪声的时间曲线示例；

[0019] 图6f是模拟头动噪声的空间分布图示例；

[0020] 图7a是在系统生理噪声干扰下，感兴趣区1和感兴趣区2的平均时间序列示例；

[0021] 图7b是在系统生理噪声干扰下，传统的相干分析方法得到的结果；

[0022] 图7c是在系统生理噪声干扰下，空间-频率连接矩阵示例；

[0023] 图7d是在系统生理噪声干扰下，本发明的方法得到的结果；

[0024] 图7e是在系统生理噪声干扰下，传统相干分析方法和本发明方法的ROC曲线对比图；

[0025] 图8a是在头动噪声干扰下，传统相干方法得到的结果；

[0026] 图8b是在头动噪声干扰下，本发明的方法得到的结果；

[0027] 图8c是在头动噪声干扰下，传统相干分析方法和本发明方法的ROC曲线对比图；

[0028] 图9a是本发明方法对空间模板范围估计过大的鲁棒性分析结果；

[0029] 图9b是本发明方法对空间模板范围估计过小的鲁棒性分析结果；

[0030] 图9c是本发明方法对空间模板位置估计偏移的鲁棒性分析结果。

具体实施方式

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

[0032] 根据本发明的方法的理论前提基于以下本领域的已有研究结论：（1）脑功能系统内部的体素之间活动相似性高，而功能系统内部体素与外部体素之间活动相似性低；（2）全局噪声会使所有体素间的脑功能连接整体升高；（3）在功能连接的主要频率上，实际试验中的全局噪声强度不足以完全淹没功能连接在功能系统内外部的差异。

[0033] 图1、本发明实现过程的流程图。如图1所示，本发明提出了一种自动判别功能连接精确频率范围的方法，它不仅能够基于数据计算出功能连接-频率分布曲线，并基于该曲线自动判别出精确的功能连接频率范围，还能够克服脑功能成像数据中全局噪声对功能连接频率分析的干扰。

[0034] 参见图1，在采集到脑功能数据后，选取感兴趣功能系统内的某一体素/区域作为种子点，计算得到每个频率点上各体素与种子点的同步性空间分布图，再利用先验的感兴趣功能系统的空间模板，对每个频率点上的同步性空间分布图进行空间加权求和，得到空SW SW间加权的同步性系数，构成功能连接频率曲线c ，最后逐频率点对c 进行统计检验，强度显著的频率即为感兴趣脑功能系统的功能连接频率范围。

[0035] 下面参照附图详细描述本发明提出的自动判别功能连接精确频率范围的方法。图1为本发明实现过程的流程图。如图1所示，根据本发明的方法，包括如下步骤：

[0036] S1：获取脑功能成像数据

[0037] 根据本发明脑功能成像数据可以是脑功能核磁共振成像数据，也可以是脑功能近红外光谱成像数据。

[0038] 脑功能核磁共振成像数据的采集在具备平面回波成像（EPI）序列的磁共振扫描仪2
上完成。成像参数设置如下，空间分辨率一般为数毫米，如3*3mm，重复时间（repetition time，TR）一般为数秒，如2s，采样时间最好多于5分钟，成像范围足够覆盖感兴趣功能区和部分功能区外部区域。

[0039] 脑功能近红外光谱成像数据的采集在功能近红外光谱成像仪上完成。成像参数设置如下：，时间分辨率一般为数十或数百毫秒，如100ms，光极间距一般为3~5cm，采样时间最好多于5分钟，成像范围足够覆盖感兴趣功能区的外部脑皮层部分和部分功能区外部区域。

[0040] 以上采集过程，如果为静息状态，要求受试者安静，并尽可能保持头部不动。

[0041] S2：预处理

[0042] 数据采集完毕后需要进行常规预处理。功能核磁共振数据的预处理包括：头动矫正、去线性漂移、空间滤波（也称为空间平滑）、空间标准化等。功能近红外数据的预处理包括：由光强数据到血氧浓度数据的解算、去线性漂移等。但上述这些过程要根据具体使用目的决定，除近红外光强数据的解算外，并非必需的过程。

[0043] S3：根据先验信息生成感兴趣脑功能系统的空间模板：

[0044] 根据先验信息确定感兴趣功能区域的空间体素位置，生成与脑功能数据匹配的感兴趣功能系统的空间模板。先验信息可以来自已有技术的研究结果，也可以来自受试者附加的结构像扫描提供的解剖信息，或功能定位任务扫描的分析结果。例如，感兴趣功能系统为感觉运动系统，对受试者进行手指运动任务的脑功能扫描，将该任务显著激活的脑区标记为1，没有显著激活的区域标记为0，即可生成感觉运动系统空间模板m。根据具体使用情况，可以将m进行适当尺度（如8mm）的空间平滑。

[0045] S4：选取种子点，计算空间-频率连接矩阵C，逐频率生成同步性空间分布图。

[0046] 选取感兴趣功能系统内的某一体素或区域作为种子点，将该体素的时间序列或者该区域的平均时间序列定义为种子点时间序列s，计算成像单元中每个体素的时间序列x与s的同步性频谱分布cxs(λ)，将cxs(λ)作为列向量，按照体素的空间位置顺序排列在一起，即可生成空间-频率连接矩阵C。图2是空间-频率连接矩阵C的示意图，如图2所示，C的第i行由频率点λi上每个体素与种子点的同步性大小构成。按照空间位置重构该行向量，即可得到频率点λi的同步性空间分布图cs(λi)，如图2c所示。

[0047] 计算同步性频谱分布可以采用相干（coherence）分析方法，也可以采用Pearson相关系数谱分解的方法。

[0048] 其中，相干（coherence）分析方法具体如下：

[0049]

[0050] 式中cxs(λ)为在频率点λ上，时间序列x与s之间的相干系数，其值为0~1，0表示在该频率时间序列之间无线性关系，1表示两个时间序列之间存在理想的线性关系。Gxs表示时间序列x与s之间的互功率谱，Gxx和Gss分别表示时间序列x和s的自功率谱：

[0051]

[0052]

[0053]

[0054] 式中，Fx和Fs分别表示时间序列x和s傅里叶频谱分析得到的频谱，*表示共轭。

[0055] Pearson相关系数谱分解的方法具体如下：

[0056] 将时间序列的点数记为N（例如N=300），Pearson相关系数计算如下：

[0057]

[0058] 根据离散傅里叶反变换法则：

[0059]

[0060]

[0061] 可将cxs分解到频率域上：

[0062]

[0063] S5：逐频率地对同步性空间分布图进行空间加权求和：

[0064] 利用步骤S3中生成的感兴趣脑功能系统的空间模板p，结合每个频率点λi上同步性空间分布图cs(λi)的特点，计算出该频率点上每个体素x对应的空间加权系数wx(λi)

[0065]

[0066] 式中M为体素总个数（例如52个）， 为m的空间平均数， 为cs的空间平均数。

[0067] 利用wx，逐频率地计算同步性空间分布图cs的空间加权和系数，即：

[0068]SW

[0069] 所得的c 即为感兴趣功能系统的功能连接频率曲线。

[0070] S6：对功能连接频率曲线逐频率进行显著性判断：

[0071] 为了避免统计检验中的第一类错误（“弃真”），这里采用非参数的多重比较方法。为了构建频域功能连接强度的零假设分布，将种子点时间序列的功率谱固定，随机变化它的相位频谱，通过傅里叶反变换法则得到相位随机化后的种子点时间序列，其他体素的时间序列不做任何变化，按照上面步骤S4的方法计算所有体素与相位随机化后的种子点时间序列的空间-频率同步化矩阵，以及对应的功能连接频率曲线。将随机化过程重复若干次（如100次），利用所有随机得到的功能连接频率曲线上的值，生成频域功能连接强度的零假设分布。逐频率点地计算感兴趣功能系统的功能连接频率曲线值的显著性水平p，设定阈值p0,其中，0

[0072] 下文提供了采用根据本发明的方法对功能连接精确频率范围进行自动判别实例的具体过程。

[0073] 本实例计划研究静息状态下感觉运动系统的功能连接频率范围。选取16名受试者参与实验，19~26岁，男8人，女8人。所有受试者都没有任何运动或者神经系统疾病。经爱丁堡利手量表测验，所有受试者均为右利手（标准为大于40）。对所有受试者进行预定的两段脑功能近红外光谱扫描。前一段为11分钟的静息状态，要求受试者安静放松闭眼，尽量保持头部不动，不进行系统地思考。后一段为双手手指序列运动任务。在进行脑功能近红外光谱扫描过程中，首先是63秒的休息（看“+”），然后是任务block，接着是基线block。共有7个任务block和8个基线block，每一个block的长度可以从20秒到28秒随机选取，但是每一个任务block和与其对应的基线block的长度是相同的。在任务block，屏幕上首先出现0.5s的“~***~”序列，然后该序列消失，出现有“F”和“J”随机顺序组成的四个字母序列，例如“FJFJ”。要求受试者看到这个字母序列之后用左手和右手的食指又快又准的按下键盘上的FJFJ，顺序和次数都要求与屏幕上的字母序列相同。1.5s后字母序列会消失，然后出现下一组刺激。设定扫描参数为：时间分辨率为10Hz，17个发射光极，16个接收光极，共组成52导，覆盖双侧感觉运动皮层及部分颞上回等其他区域。具体光极摆放位置如附图3所示。

[0074] 对采集到的功能近红外光谱数据进行如下预处理，包括：光强数据向脑血氧浓度数据的转换，去线性漂移。然后对双手手指序列运动任务数据进行去除前30s的不稳定数据，高通滤波（0.017Hz），逐导进行广义线性模型（GLM）分析，最后在受试者群体水平上进行单样本t检验，显著性水平p<0.05（经过false discovery rate多重比较校正）的导为任务显著激活的导。

[0075] 将所有任务显著激活的导标记为1，未显著激活的导标记为0，并利用3*3的矩阵：[0.25，0.5，0.25；0.5，1，0.5；0.25，0.5，0.25]为平滑核进行空间平滑，生成感觉运动系统的空间模板m。

[0076] 将定位任务激活最显著的第24导定义为种子点。提取第24导时间序列作为种子点时间序列，利用相干分析方法计算其与其他51个导静息状 态下的同步性频谱分布。将每个导的同步性频谱分布曲线作为列向量，按照导的编号顺序排列在一起形成空间-频率连接矩阵C。提取空间-频率连接矩阵的一行，按照空间位置重构，即得到该频率的同步性空间分布图。图4显示了4个频率点上的同步性空间分布图。

[0077] 利用空间模板m和空间-频率连接矩阵C，根据公式（7）计算出每个频率上每个导的空间加权系数，逐频率计算同步性空间分布图（即空间-频率连接矩阵的行向量）的空间加权和，生成功能连接频率曲线。最后对功能连接频率曲线逐频率进行显著性判断（p<0.05）。如图5a所示，感觉运动系统的静息状态下功能连接频率范围为0.01-0.0732Hz。

[0078] 为了对比，这里也呈现了传统的判断功能连接频率范围的方法分析的结果。基于相同的静息态脑功能近红外光谱成像数据，分别平均左右两侧感觉运动系统区域（如图3所示）所有导的时间序列作为左右两侧感觉运动区的时间序列，并利用相干分析计算两条时间序列之间的同步性频谱分布，将同步性显著高的频率作为功能连接频率范围。如图5b所示，基于传统方法，功能连接的频率范围被判定到了0-0.1025Hz和1.3Hz左右。

[0079] 从上文将传统方法与本发明方法结果的比较发现，本发明抑制了传统方法中存在的多处由噪声引起的虚假同步，例如<0.01Hz的超低频同步（可能与数据漂移有关，也有可能与生理的超低频噪声有关），0.1Hz左右的低频同步（可能与0.1Hz左右的血压和心率波动的变化，及Mayer’swave有关），以及1.3Hz左右的高频同步（可能与心跳波动有关）。这说明本发明既可自动精确判别功能连接频率范围，又可克服数据中的全局性噪声对功能连接频率分析的影响。

[0080] 下面基于本发明的方法，进行了一系列的模拟实验。其中分方法部分和结果部分来阐述。

[0081] 一、方法部分：

[0082] 模拟数据D由三个主要部分构成：脑自发波动DSF，全局噪声DGN，以及机器噪声DIN，SF GN IN即D＝D +D +D 。时间序列长度为15分钟，采样率为10Hz。观测矩阵为10×25。

[0083] 脑自发波动DSF根据以下模型生成：

[0084] DSF＝TSF×SSF （9）

[0085] 其中，如图6a中所示，TSF是低频自发波动的基本时间序列，它由0.01-0.1Hz的相SF位随机的正弦波线性叠加构成，在0.01-0.1Hz段，频谱分布服从0.69*1/f。S 表示自发波动的空间分布。如图6b所示，它包含两个方形区，表示两个感兴趣区，共同构成感兴趣功能SF
系统。S 为二值化矩阵，即感兴趣区内为1，感兴趣区外为0。式（6）中符号“×”表示向量外基运算。

[0086] 机器噪声DIN由信噪比为1的随机高斯白噪声生成，每个观测导之间相互独立。

[0087] 全局噪声DGN包含两类，即系统生理噪声DGN-1和头动噪声DGN-2。他们分别根据各自噪声的时-空间特点模拟产生。具体方法如下：

[0088] 1.第一类全局噪声——系统生理噪声DGN-1

[0089] 系统生理噪声主要包括由系统性生理活动，如心跳、呼吸以及一些低频的血管运动等引起的血氧活动。为了评估这些复杂的系统生理噪声可能对检测功能连接频率范围产GN-1生的影响，我们在每个空间位置k按照如下模型模拟生成D ：

[0090] DkGN-1=TGN-1(t+PkGN-1)×SkGN-1 (10)

[0091] 其中，TGN-1表示系统生理噪声的基本时间序列，如图6c所示。已有研究表明，系统生理噪声的频率范围跨度很广，其中极低频波动（~0.04Hz），动脉血压及心率的变化波动SF（~0.1Hz），以及呼吸相关的波动（~0.25Hz）可能会和脑自发波动D 的频率有重叠。而心跳SF
相关的波动（~1Hz）相对高频，基本不会与D 有频率重叠。因此，为了全面考量信号频率重叠和频率不重叠两种情况，我们模拟生理噪声的频率范围为0-0.15Hz。这个频率范围不仅GN-1
覆盖脑自发波动的频率范围（0.01-0.1Hz），还比它宽，有频率不重叠部分。T 由幅度、相位均随机的0-0.15Hz的正弦波线性叠加生成。对于系统生理噪声的空间分布，已有研究显示，尽管观测导之间存在很强的一致性，不同观测导之间在相位和振幅上还是存在一定的差异。 因此，我们模拟了相位空间分布图和振幅空间分布图。他们在每个观测导上都有不同的值。具体地，他们通过被二维高斯核（半高全宽为10）空间平滑后的高斯白噪声矩阵生GN-1 GN-1
成。P 矩阵的幅值被限定在0到2之间。S 矩阵的幅值被限定在1/3到1/2之间。为了尽量全面地模拟出实际实验中系统生理噪声的复杂的时空间特性，随机模拟实验重复了
100次。

[0092] 2.第二类全局噪声——头动噪声DGN-2

[0093] 为了模拟头动噪声可能对功能连接频率范围判别的影响，我们按照以下模型模拟GN-2生成头动噪声D ：

[0094] DGN-2=TGN-2×SGN-2 (11)

[0095] 其中，TGN-2是头动噪声的基本时间序列，如图6e所示。它由时间上随机出现的GN-21-12个脉冲信号构成。这个头动噪声的出现频率与实际实验中接近。S 是头动噪声振幅的空间分布图。如图6f所示，在不同观测导之间，头动噪声幅值略有不同，这是由于实际实验中头动噪声是由于光极帽与头皮的随机相对位移造成的，而不同位置光极帽的曲率不一样，由此可能造成头动噪声的幅度在不同位置不一样。该空间分布图通过被二维高斯核（半高全宽为10）空间平滑后的高斯白噪声矩阵生成。矩阵的幅值被限定在3到10之间。为了尽可能全面的模拟出实际实验中头动噪声时间上随机产生的特性，以及不同被试间的差异，头动噪声被随机模拟100次。

[0096] 3.空间模版

[0097] 在本发明中，空间模版m对去除全局噪声干扰有很关键的作用。因此，本发明对m的鲁棒性非常重要。为了模拟m对功能系统大小的估计误差，我们将模拟的功能系统范围在原有基础上分别扩大或缩小0%，10%，以及20%。另外，为了模拟m中对功能系统位置的估计误差，我们将模拟的功能系统位置在原有位置上分别水平移动0%，10%，以及20%。对于每种情况，都估计相应的功能连接频率范围，并度量m的误差对结果的影响。

[0098] 为了定量估计所得功能连接频率范围的检测结果的准确性，我们在模拟实验中按照如下公式计算了检测结果的特异性（specificity）和敏感性（sensitivity）：

[0099]

[0100] (12)

[0101]

[0102] 其中 表示检测得到的功能连接频率范围， 表示真实的功能连接频率范围，在这里我们模拟为0.01-0.1Hz。I表示特异性和敏感性计算的频率段，这里我们设定为0-0.25Hz。另外，考虑到阈值选择可能对功能连接频率范围检测结果的影响，我们还采用了受试者工作特征曲线（receiver operator characteristic(ROC)）进一步分析该方法的特异性和敏感性。

[0103] 二、结果部分：

[0104] 1、第一类全局噪声——系统生理噪声DGN-1

[0105] 在第一部分的模拟实验中，数据按照D＝DSF+DGN-1+DIN生成，用于度量系统生理噪声对功能连接频率范围检测的干扰。

[0106] 为了与本发明的方法比较，也给出了利用传统的相干方法检测功能连接频率范围的相关结果。传统的相干方法直接将功能系统内两个代表性区域（即模拟的两个感兴趣区）中的平均时间序列之间的相干程度cROI1-ROI2作为该系统的功能连接频率特性。100次随机模拟实验的平均cROI1-ROI2曲线在图7b中呈现（方差为图中细黑线）。如图7b中下方黑色横条所示，传统方法得到的平均静息功能连接频率范围是从0.0001Hz到0.1510Hz，上下限标准差分别为0.0008Hz和0.0032Hz。该范围不仅在超低频（图中的左侧浅灰色部分）处而且在高频（图中的右侧浅灰色部分）处都大大超过了实验设定的功能连接频率范围的金标准（0.01-0.1Hz，图中的中部深灰色部分）。检测得到的功能连接频率范围基本与实验模拟的系统生理噪 声的频率范围（0-0.15Hz）吻合。检测结果的特异性为0.591±0.020（均值±标准差），敏感性为1。这说明，利用传统的相干方法，系统生理噪声会明显干扰功能连接频率范围检测。

[0107] 图7c为利用本发明，种子点选择左侧感兴趣区的正中间的体素，得到的空间-频率连接矩阵及其低频部分放大显示。基于该连接矩阵，以及如图6b中显示的空间模版m，得到功能连接频率曲线cSW。如图7d中所示，得到的功能连接频率范围（图中下方黑色横条）平均从0.0090Hz到0.1033Hz，上下限标准差为0.0026Hz及0.0041Hz。该结果与实验设定的功能连接频率范围的金标准（0.01-0.1Hz，图中的中部深灰色部分）非常接近。检测结果的特异性为0.931±0.026（均值±标准差），远远高于传统的相干方法。检测结果的敏感性与传统的相干方法类似，为0.999±0.011。这说明本发明可以成功抑制全局生理噪声的不良影响。另外，如图7e所示，本发明的ROC曲线下面积为0.975，显著高于传统相干方法的ROC曲线下面积0.862。该结果进一步证明了本发明与传统方法比，对于功能连接频率范围检测具有更高的特异性和敏感性。

[0108] 2.第二类全局噪声——头动噪声DGN-2

[0109] 在第二部分的模拟实验中，数据按照D＝DSF+DGN-2+DIN生成，用于研究头动噪声的干扰。

[0110] 传统的相干方法得到的cROI1-ROI2曲线在图8a中呈现（指数坐标显示）。它明显覆盖了很广的频率范围，得到的功能连接频率范围（图中下方黑色横条）平均位于0.0002Hz和0.754Hz之间，上下限标准差为0.0014Hz和0.173Hz。该范围远远广于实验设定的功能连接频率范围的金标准（0.01-0.1Hz，图中的中部深灰色部分）。检测结果的特异性为0.037±0.178（均值±标准差），敏感性为1。对头动噪声的频率谱分析显示，传统相干方法得到的功能连接的频率范围与头动噪声的频率范围基本一致（结果没有显示在图上）。这说明当利用传统相干方法时，头动噪声会严重干扰功能连接频率范围检测。

[0111] 本发明得到的cSW曲线在图8b中呈现。它集中在低频部分，得到的功能连接频率范围位于0.0077±0.0039Hz和0.113±0.032Hz之间。特异性为0.866±0.197，显著高于传统方法。敏感性为1，与传统方法一样。另外，如图8c所示，本发明ROC曲线下面积为0.976，高于传统方法的0.850。进一步证明本发明相对于传统方法对克服全局头动噪声的有效性。

[0112] 3.空间模版

[0113] 本发明对空间模版m的鲁棒性分析结果显示在图9中。如图9a和9b所示，当对SW功能系统大小的过大或过小估计时，尽管c 的强度会略有降低，但是最终的功能连接频率范围结果基本还是与理想情况一致。类似的，如图9c所示，当对功能系统位置偏差估计时，最终的功能连接频率范围结果基本还是与理想情况一致。因此，实验结果表明，本发明对空间模版对功能系统大小及空间位置的估计具有一定鲁棒性。

[0114] 以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。