用于构建脑磁图网络的信号提取方法

杨春兰, 吴文晓, 吴水才, 任洁钏

(1.北京工业大学生命科学与生物工程学院，智能化生理测量与临床转化北京市国际科研合作基地，北京 100124；2.首都医科大学附属北京天坛医院，北京 100050)

作为人体最重要的器官之一，大脑的功能多样，结构复杂，人体的多种生理功能都与大脑的神经活动紧密相关. 成年人的大脑由大约1011个神经元细胞组成，众多神经元细胞之间通过约1015个神经突触相互连接进而形成一个高度复杂的脑网络. 美国著名的脑网络分析专家Sporns[1]教授于2005年提出脑连接组的理论并指出：根据元素和连接不同，脑连接组可划分为微观尺度、介观尺度和宏观尺度3种不同级别. 微观尺度即神经元细胞和神经突触之间的连接；介观尺度即神经元集群之间的连接；宏观尺度即各个大脑分区之间的连接. 大脑的功能网络可以反映大脑在处理不同任务时的协同关系，这为理解大脑的工作模式提供了重要的支撑. 由于当前科技水平限制，人类对大脑的研究主要集中在宏观尺度,大脑功能连接网络的构建主要通过脑电图(electroencephalogram,EEG)、脑磁图(magnetoencephalography,MEG)及功能磁共振成像(functional magnetic resonance imaging,fMRI)等成像技术来实现[2-8].

MEG是一种磁源成像技术，在磁屏蔽室中利用超干量子干涉仪(superconducting quantum interference device,SQUID)[9]来测量大脑外部的各个点处的磁场，SQUID阵列有多个传感器通道，用来记录被试头部周围磁场分布并将其绘制成图. MEG有亚毫米级的空间分辨率和毫秒级的时间分辨率，能够解析脑区域之间的活动；MEG源映射的准确性不受头部组织复杂分层造成的信号畸变影响.

MEG相较于fMRI、EEG是一种新的功能成像技术，随着近年来MEG分析方法的涌现，利用MEG分析大脑功能网络的研究也不断增加. Hillebrand等[10]于2012年提出了一种基于图谱的，以大脑分区为节点的MEG功能网络分析框架，利用beamformer技术将MEG信号投影到大脑皮层，然后选择功率最大的信号代表此脑区，以此来构建脑网络矩阵. 2017年，López 等[11]通过研究29名轻度认知障碍患者(mild cognitive impairment,MCI)和29名健康老年被试者的脑磁图网络发现，MCI患者的脑网络效率较低，且验证了患有MCI的受试者得阿尔兹海默症的风险会增高. Briley等[12]于2018年对48名被试者(9～25岁)的MEG数据进行网络研究，发现网络内部和网络之间的振荡的协调与认知技能的成熟相关，网络的连接性随着年龄的增加而增加. 2018年，van Nieuwenhuizen等[13]分析20名经过肿瘤切除的脑膜瘤患者以及20名健康被试者的MEG数据发现，脑膜瘤患者的工作记忆能力较低，同时发现认知功能与默认网络中的功能连接和脑膜瘤患者的中枢病理学相关. 利用MEG成像技术构建功能网络，网络节点的选择是个关键的问题，常见的节点有2种：以传感器位置作为网络节点和以脑区作为网络节点. 前者利用大脑外部传感器采集到的信号直接构建功能网络，这是一种简单直接但却存在较大误差的方法；脑区作为网络节点首先需要根据采集到的脑磁时序信号重建脑内源信号，然后以脑区为单位，通过源信号构建功能网络.

选择代表性神经源信号作为局部脑区信号是构建脑网络的关键，也是本文讨论的主要问题，目前多数研究都是选择脑区中功率最大的神经源信号作为代表性信号[10-11]，这会损失部分区域重要细节信息；也有研究使用脑区质心处的神经源信号作为整个区域的代表性信号[14]，质量中心在几何学上有效性高，但是在信号处理方面其有效性尚未得到证实.

MEG功能网络在源重建方法和网络分析框架上都尚未达到成熟阶段，但多项研究结果已经表明利用MEG研究大脑功能网络的可行性. 本文采用Hillebrand基于图谱的功能网络的分析框架，对大脑分区内信号的提取方法提出了2种改进方案，即叠加平均和聚类的方法，然后对构建的脑网络矩阵进行k均值聚类分析，通过聚类准确率实现与最大功率方法的比较. 最后，通过分析3种方法构建的脑网络特征和小世界属性，对3种方法进行比较和评价.

1 方法与材料

1.1 数据来源

本研究选取了人脑连接组计划(human connectome project，HCP)数据库中51名健康被试者的MEG数据，其年龄在22～35岁，数据来源于圣路易斯大学医学院的全脑MAGNES3600(4D Neuroim-aging,San Diego,CA)系统. 每名被试者的数据均被采集了包括工作记忆、语言处理和运动3种任务状态的MEG数据[15]. MEG系统包括248个磁力仪通道和23个参考通道，在白噪声范围(2 Hz以上)内，磁力仪的均方根噪声为5 fT/sqrt(Hz),数据的采样频率为2 034.51 Hz，心电、眼电、肌电信号与脑磁信号同步采集，所有电极的接触电阻控制在10KOhms之内. 具体任务设计信息可参阅HCP手册(https:∥www.humanconnectome.org/study/hcp-young-adult/article/announcing-1200-subject-data-release/).

1.2 预处理

选用的数据都是经过了HCP数据库的预处理，可参阅HCP手册. 预处理包括数据检查、去除含噪声的数据段、去除有问题的通道、ICA分析. 通过ICA分析可以去除MEG数据中掺杂的心电和眼电伪迹，此外，所有的数据都经过了1.3～150.0 Hz的带通滤波和59～61 Hz、119～121 Hz的Butterworth滤波. 被试的MRI T1结构像与MEG进行配准之后，会被标准化到标准模板上，然后将AAL分区模板逆变换到被试个体空间上，最后得到被试个体空间上的分区信息.

1.3 MEG功能网络构建方法

本研究中MEG功能网络的构建框架是一种基于图谱的脑功能网络. 此网络是以大脑分区作为网络节点，信号之间的相关系数为边构建而成的[10]. 网络构建流程如图1所示.

首先对被测试者进行数据采集，包括：临床数据、脑磁图数据、磁共振 T1像. 将被测试者的MRI T1结构像与MEG数据进行配准，然后将配准后的图像在空间上标准化到MRI标准模板上，标准化后的MRI图像和AAL分区模板的感兴趣区域(region of interest，ROI)通过逆变换转换为个体分区的MRI图像；通过对配准后的MRI结构像分割从而构建得到被试者大脑的头部模型[16]；由脑磁图数据产生的协方差矩阵、被试者的头模型以及个体分区的MRI图像一起用来计算beamformer权重. 之后MEG数据通过对beamformer权重的投影产生时间序列[10]；对每个频段的时间序列进行相关性分析得到连接矩阵，最后对连接矩阵进行图谱分析构建大脑网络.

对于感兴趣区域的选取，本文选择了AAL模板[17]，此模板共有116个脑分区，其中90个属于大脑，在模板空间中选取了前78个大脑分区作为ROI. 每个ROI内包含多个体素，而且体素的数量也是各不相同的，这依赖于空间采样，在本研究中采用的是体素间隔为8 mm的源模型.

源信号的重建过程基于FieldTrip(http:∥www.fieldtriptoolbox.org/)工具包完成，源信号重建流程如图2所示.

源重建流程图中调用的函数都来自FieldTrip工具包. 每个被试者的MRI图像都是不同的，利用ft_read_mri函数读取被试者的磁共振图像，然后通过ft_prepare_headmodel函数构建源重建过程中需要的体积传导模型，也称为头模型，如图3所示. 利用ft_timelockanalysis函数计算出协方差矩阵. 最后在ft_sourceanalysis函数的method属性中选择线性约束最小方差(linearly constrained minimum variance,LCMV)方法[18]计算得到beamformer权重矩阵.

当确定beamformer的权重后，就可以重建每个体素的时间序列，传感器的非均匀投射可以在源信号上得到显示，权值随深度增加，而传感器的噪声保持不变. 为了消除这固有的偏差，在计算源信号之前需要对beamformer的权重利用它本身的向量范数进行归一化. 文献[19]表明α频带最活跃的区域是视觉皮层和顶叶区域，β频带的功能连接在感觉运动区域最强[20],γ频带在听觉皮层区域[21-22]和视觉皮层区域[23]有功能性作用. 因本研究是基于运动、工作记忆和语言处理3种任务态的，而这3种任务状态与视觉皮层、听觉皮层以及运动皮层区域最为相关，因此，在α(8～13 Hz)、β(13～30 Hz)、γ(30～48 Hz)以及全频带(1～48 Hz)共4个频带对体素的时间序列进行分析.

通过体素间距为8 mm的标准源模型进行源重建后会产生约10 000个体素源信号，对每一个体素的时间序列进行研究工作量巨大且也没有必有，因此多数研究中都会选取代表性的时间序列来表示整个ROI的时间序列，对于每一个频带和ROI，通常选取在此频带中功率最大的体素作为整个ROI代表，并使用这个体素的时间序列进行下一步研究分析[10-11].

1.4 改进方案

在用MEG构建脑网络时，常用最大功率方法提取代表性时间序列作为整个ROI的时间序列，但此方法会忽略掉其他体素对整个ROI的作用，导致信息量损失. 因此，提出了2种改进方案：基于叠加平均的方法和基于聚类的方法.

1.4.1 基于叠加平均的方法

脑区源信号具有方向性，直接对其在矢量度量下进行叠加平均会导致一些信号的抵消现象. 本研究的目的是通过计算脑区之间的一致性，构建脑网络，所以源信号的方向性可以忽略，只考虑其在主要方向上的投影即可. 选取某ROI的任意体素源信号作为基准信号，分别计算各体素源信号和基准信号的相关性，并根据其正负性对各体素源信号进行符号修改. 最后对修改后的所有信号进行叠加平均得到ROI的时间序列信号

S=〈sign(Cθ)Sθ〉

(1)

式中：S为ROI的时间序列信号；Cθ为基准信号与θ处体素源信号之间的相关性；Sθ为θ处的体素源信号；〈〉为取平均.

1.4.2 基于聚类的方法

人类大脑在处理某个任务时不止是单个脑区在工作，而是多个脑区共同协调作用，大脑的某个脑区产生的源信号也不仅仅只代表一种生理活动. 聚类根据各个源信号之间的距离把信号分成一簇簇的类，使得类内的源信号的相关性较强，而类间的相关性较弱. 不同的类别内包含的体素数量是不同的，某一类中体素的数量越多则说明其在脑区活动中的重要程度越高. 因此，可以选用包含体素数量最多的类别代表整个脑区的生理活动，而把此类别的聚类中心作为整个脑区的代表信号.

某些聚类方法需要事先设定分类数，例如在k均值聚类中需要知道k的具体值，才能进而将要聚类的对象分为k类. 无法事先确定每个脑区信号分类数目，所以需要事先设定分类数目的聚类方法在此处不适用，本研究利用信号间的相关性进行聚类. 具体方法如下：

假设某个脑区内所有体素信号的集合用N表示，其元素ni表示某一体素的信号，Mi表示聚类后得到的第i类的信号的集合.

在集合N中任取元素ni作为第一类信号的一个元素,则

N=N-{ni}

(2)

M1={ni}

(3)

取集合N中的任意元素nj，令N=N-{nj}，分别计算nj与已有分类Mk中所有信号ml的相关性

c=Cov(nj,ml),ml∈Mk

(4)

max(c)>θ，则将nj划为第k类

(5)

max(c)<θ,则将nj划为新的一类k1类

(6)

重复此过程直到遍历所有体素信号，每个脑区中会聚类出多种不同类别，而这不同的类别代表着不同的功能活动，选取此脑区中体素数最多的那一类别的聚类中心的时间序列信号代表整个脑区即ROI的时间序列信号.

1.5 脑网络构建

通过权重矩阵就可得出体素的时间序列，然后分别通过最大功率、叠加平均和聚类的方法获取代表ROI的时间序列，通过量化时间序列之间的相位关系，可捕捉到活动源之间的函数相关关系. 本研究使用了相关系数来估算ROI之间的功能连接.

相关系数量化了信号之间相关性关系，且

(7)

1.6 脑网络属性分析

对于每种任务状态，都有多次试验，经过筛选，本研究选取了其中50次试验，ROI数为78，因此对于每种任务状态下的每个频段利用每种信号提取方法都可以得到50个78×78的网络矩阵. 3种任务状态在相同频段相同信号提取方法下共得到150个网络矩阵，将工作记忆任务状态下的网络矩阵标号为1～50，运动任务状态下的网络矩阵编号为51～100，语言处理任务状态下的网络矩阵编号为101～150. 对上述150个矩阵利用k均值聚类分为3类，在相同任务状态下产生的网络矩阵具有较高的相似程度，若编号为1～50的矩阵为一类，51～100的矩阵为一类，101～150的矩阵为一类，那么就可以认为利用此种特征提取的方法得到的矩阵k均值聚类准确率为100%.

除此之外，本文还对脑网络特征和小世界属性进行研究分析. 小世界网络是一种介于规则网络和随机网络之间的一种网络属性[24]，自从提出了小世界网络的概念，人们开始研究各种复杂网络的小世界属性，有关脑网络小世界属性的研究也越来越多[25-26]，本研究针对3种信号提取方法得到的网络矩阵进行网络特征的分析，分别采用2个网络特征(聚类系数、最短路径长度)以及3个小世界特征(λ、γ、σ)来分析3种信号提取方法得到的功能网络特征之间的差异. 网络特征的分析用到工具包brain connectivity toolbox(BCT,http:∥www.brain-connectivity-toolbox.net). 特征参数定义及计算公式如表1所示.

表1 网络特征参数定义及计算公式

2 实验结果与分析

2.1 脑网络矩阵k均值聚类准确率的比较

在进行脑网络分析之前，以在工作记忆状态下ID为104012的被试为例，对其α频段以叠加平均方法得到的脑网络矩阵进行了二值化，阈值设置为0.8，并对二值化后的脑网络矩阵进行了可视化，如图4所示.

在表2的k均值聚类准确率结果中，每一个提取方法在每一个频段都对应2个百分数，分别为平均聚类准确率和标准差. 可以看出，叠加平均信号提取方法的聚类准确率在4个频率段上均为最高，且对应的标准差相对较低；最大功率方法在γ频段的准确率与叠加平均方法相当，而在其他3个频段比叠加平均方法低. 从频段上对比3种信号提取的方法可以看出，在α频段，3种方法的准确率均为最低，而标准差却相对较高；在γ频段上，3种方法都相对较好的表现.

表2 k均值聚类准确率

2.2 网络特征差异

以标号为104012被试者为例，该被试在工作记忆状态下的脑网络的网络特征及小世界属性如图5所示，可以看出，3种信号提取方法在4个频段均表现出了小世界属性，即λ的值接近1、γ的值远大于1，且σ的值大于1，但是不同的信号提取方法在不同的频段小世界属性强弱是不同的.σ是由Humphries等[27]提出的一种衡量小世界属性的标量，它被定义为σ=γ/λ，当σ>1时网络具有小世界属性，且σ越大，网络的小世界属性越强. 在全频段(1～48 Hz)，最大功率的信号提取方法得到的脑网络的小世界属性最强；叠加平均方法得到的脑网络的小世界属性次之；聚类方法得到的脑网络的小世界属性最弱；在α(8～13 Hz)频段，叠加平均方法得到的脑网络的小世界属性强度略高于最大功率方法得到网络的小世界属性，聚类方法得到的脑网络的小世界属性最弱；在β(13～30 Hz)频段，3种信号提取方法得到的小世界属性的强弱依次为，最大功率最强，叠加平均次之，聚类最弱；在γ(30～48 Hz)频段，依旧是最大功率的小世界属性最强，叠加平均次之，聚类最弱.

2.3 讨论分析

用脑磁图构建脑网络时，在ROI中选取代表性神经源信号是关键步骤. 本研究根据3种信号提取的方法分别构建大脑网络，然后通过对脑网络进行均值聚类分析，会获得3种信号提取方法的k均值聚类准确率，通过准确率的比较发现，在各个频段，通过叠加平均和最大功率的方法获取的时间序列构建的大脑网络，k均值聚类准确率都在70%左右，而通过聚类的方法获取的时间序列构建的大脑网络的准确率在65%左右，尤其是在α频段，聚类方法的准确率不到50%. 由准确率的比较，可以初步得出结论：利用叠加平均和最大功率的方法获取的时间序列在一定程度上更能代表整个ROI的活动情况，而利用聚类的方法获取的时间序列不能很好地代表整个ROI的活动情况.

最大功率、叠加平均以及聚类的信号提取方法获取时间序列构建的脑网络具有小世界属性，这也验证了之前研究中脑网络具有小世界属性的结论. 在对脑网络进行网络特征及小世界属性分析时，对矩阵进行了二值化，为了使网络中的节点都连通，此过程需要设置阈值，阈值设置得越大，所有节点都连通的可能性越大，同时网络的小世界属性越弱. 图5中的网络特征及小世界属性均是在最小阈值下得到的，其中最大功率和叠加平均方法得到的网络矩阵的最小阈值为0.3,而聚类方法得到网络矩阵的最小阈值为0.55. 从图5中也可以看到，在各个频段，叠加平均和最大功率方法对应的σ值都相差不多且在值的大小在1.5以上，而聚类方法对应的σ值在1.2上下. 聚类方法获取时间序列构建大脑网络的小世界属性较弱.

本研究中采用的聚类方法是通过比较新体素和已有类别中体素的相关性来划分类别的，理论上其比基于距离来划分类别的聚类方法有更高的准确度. 选取体素数最多的类别的聚类中心作为整个ROI的时间序列，如果各个类别的体素数量都相差不多，这也将会产生较大的误差；在同一脑区中的各个神经源信号的功率是不同的，选取功率最大的信号作为代表信号固然会损失部分信息，但是相比与其他信号功率最大的信号更能代表脑区，其有效性也得到了相关研究的证实；叠加平均的方法是把整个脑区的信号求和取均值，这样得到的信号在理论上能代表局部脑区. 上述实验结果也验证了基于最大功率和叠加平均作为信号提取方法的有效性.

3 结论

1) 从k均值聚类准确率和网络特征及小世界属性的分析比较中可以发现，聚类方法不适合获取ROI代表性时间序列，而最大功率和叠加平均作为获取ROI代表性时间序列的方法有较好的表现，在某些频段叠加平均方法优于最大功率方法，在另外一些频段最大功率方法优于叠加平均方法. 从大量的神经源信号中获取代表性时间序列对于构建大脑网络至关重要，不同的神经源信号作为ROI的时间序列构建出的大脑网络截然不同，所以性能良好的信号提取方法对于脑磁图构建网络及网络的分析有着很大的作用.

2) 对于如何从众多神经源信号选取代表性时间序列没有金标准，有研究使用ROI的质量中心作为整个区域的代表性时间序列[11]. 信号提取方法作为脑磁图网络构建中的关键步骤，对于整个网络构建的质量有很大影响，本研究针对这一点，对比了常用的最大功率方法以及本文提出的叠加平均方法，发现2种方法均可作为获取代表性时间序列的方法，为脑磁图网络的构建及分析提供良好的支持.