CN113317803A

CN113317803A - 一种基于图论和机器学习的神经症特征提取方法

Info

Publication number: CN113317803A
Application number: CN202110625900.5A
Authority: CN
Inventors: 宁兆龙; 孙兰芳; 王小洁; 胡希平; 郭毅
Original assignee: Chongqing University of Post and Telecommunications
Current assignee: Chongqing University of Post and Telecommunications
Priority date: 2021-06-04
Filing date: 2021-06-04
Publication date: 2021-08-31
Anticipated expiration: 2041-06-04
Also published as: CN113317803B

Abstract

本发明公开了一种基于图论和机器学习的神经症特征提取方法，首先基于EEG数据构建神经症下躯体维度和情感维度的脑功能网络，以62通道作为网络的节点，分别选取相位滞后指数和加权相位滞后指数作为网络中边的度量。提取两指标对应的网络拓扑属性——全局效率，进行融合后作为识别躯体维度和情感维度的特征向量，并利用机器学习分类模型进行分类，且该特征相应的脑功能网络映射到3D脑模型上可定位两个维度的网络异常连接位置。本发明可进一步用于为神经症的识别提供了标识物，并进一步研究疾病的靶点。

Description

一种基于图论和机器学习的神经症特征提取方法

技术领域

本发明针对神经症下的躯体维度及情感维度，提出一种基于图论构建EEG脑功能网络和机器学习方法相结合的神经症特征提取方法。

背景技术

慢性病的发病率、死亡率及致残率居高不下，且趋于年轻化，其具有极强的隐匿性，早期症状并不明显，诸如焦虑症、抑郁症、神经衰弱等各种精神障碍均统称为神经症，属于慢性病的一种，寻求针对疾病且有效的标识物成为亟待解决的问题之一。人类大脑本身就是一个极为复杂的系统，神经元之间错综复杂的连接，基于图论^[1]的脑功能网络^[2]分析方法能够深入了解大脑的结构、各个功能区域的工作机制以及彼此之间是如何相互作用与协调的。目前，大多是基于功能磁共振成像^[3]及磁共振成像^[4]进行脑功能网络的研究。但无论是在设备的购买、维护，还是所需的检查费用方面都极为昂贵。脑电图^[5](Electroencephalogram，EEG)不仅价格便宜、便携性强，而且时间分辨率极高，可提取不同频段的信息进行分析，为临床疾病的诊断提供多种衡量指标。传统的基于时域、频域或时频域进行独立通道或脑区的EEG特征提取，如信息熵^[6]、功率谱^[7]等仍被广泛研究，但疾病的产生通常并不是独立通道或脑区的功能损伤，而是不同功能脑区间的相互作用与协调产生了异常，且这些特征不能有效解决共同源问题。因此，需要研究人员提出一种基于EEG的神经症特征提取与识别方案。

发明内容

本发明的目的主要是针对现有研究的不足之处，将神经症细分为躯体维度及情感维度，设计了一种针对神经症的特征提取算法，将基于图论的EEG脑功能网络构建与机器学习算法相结合提取神经症特征。

为了实现上述目的，本发明采用的技术方案是：一种基于图论和机器学习的神经症特征提取方法，包括以下步骤：

(1)采集EEG数据，对其进行预处理。

(2)针对步骤(1)预处理后的数据，使用傅里叶变换将时域信号转换成频域信号，以获取感兴趣的频段。

(3)提取基于相位滞后指数(PLI)和加权相位滞后指数(WPLI)的特征矩阵。

(4)根据步骤(3)所得的特征矩阵构建脑功能网络：以62个通道作为网络的节点，分别以步骤(3)所得的两种特征矩阵度量节点间边的连接程度，在稀疏度为15％-45％，步长为5％的阈值下分别构建对应的脑网络，脑网络需满足小世界性。

其中，经过预处理的EEG信号，能够按照通道规则分布且没有明显伪迹，这利于后续的分析。相位滞后指数能够有效最小化共同源问题(容积传导或活跃参考电极)，而加权相位滞后指数在相位滞后指数的基础上引入了权重的概念，从而减少对其他无关噪声的敏感性，同时，能够提高其检测相位同步性变化的能力。二者相结合能够很好的对神经症下的两个维度作出区分。人的大脑实际上是由多个功能脑区相互连接形成的复杂网络，基于图论方法构建脑功能网络，是研究脑疾病网络异常拓扑结构的重要手段之一。

附图说明

图1为EEG数据的采集和预处理流程；

图2为基于EEG的脑功能网络构建、特征提取流程；

图3为基于PLI的躯体维度和情感维度脑功能网络属性显著性差异结果；

图4为基于WPLI的躯体维度和情感维度脑功能网络属性显著性差异结果；

图5为基于PLI所得特征，并利用三种分类器进行分类的AUC面积；

图6为基于WPLI所得特征，并利用三种分类器进行分类的AUC面积；

图7为基于不同组合特征的三种分类器的分类结果，证明了本方法提取特征的有效性；

图8、图9为基于不同度量指标所得特征，躯体维度和情感维度映射于3D脑模型的脑功能网络连接情况；其中，图8基于PLI，图9基于WPLI；且(a)为躯体维度组脑功能网络连接；(b)为情感维度组脑功能网络连接；且与躯体维度组相比，(c)为情感维度组连接增强的部分，(d)为情感维度组连接减弱的部分。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面将对本发明的具体实施方式作进一步的详细描述。

本发明实例提供了一种基于图论和机器学习的神经症特征提取方法，参见图2和图2，该方法包括基于EEG数据的脑功能网络构建、拓扑特征提取和识别，以及神经症下躯体维度和情感维度的脑功能网络连接定位。

步骤1：EEG数据的采集

EEG数据的采集均在一间室内噪声小于20dB，能够隔绝电磁波和交流传导干扰/AM/FM电波干扰的电声屏蔽室内完成。在采集信号时以REF为默认参考电极，GND为地电极，其余62个电极按照国际标准10-20系统的位置安放。一些参数设置如表1所示。

表1参数设置

确定各项参数后，开始EEG记录。整个记录过程共持续19分钟，包括闭眼状态8分钟，睁眼状态8分钟以及深呼吸3分钟。

步骤2：EEG数据的预处理

EEG数据的预处理包括剔除无用电极、重参考、滤波、剔除坏段、分段和基线校正、降低采样率以及基于独立成分分析的伪迹去除，最终得到高质量EEG数据。

Step1：剔除无用电极：以REF作为参考电极，所有其他电极与REF的电位差作为该电极的记录数值。由于REF和GND的信息不会用于后续的分析，将其剔除掉。

Step2：重参考：以双侧乳突作为重参考的位置，该部位不易受头部运动的影响，对应的电极是TP9和TP10。

Step3：滤波：采集数据时是在0.1Hz-50Hz间进行带通滤波，但一般是使用0.1Hz-45Hz，在此基础上，分别提取不同节律对应的频段进行分析。

Step4：剔除坏段：手动对每个数据进行检查，剔除掉明显的肌电以及其他因素造成的伪迹。

Step5：分段和基线校正：使用2s的时间窗对数据进行分段。

Step6：降低采样率：利用EEGLab工具包中的pop_resample方法将采样率降为100Hz，即1s内记录100个数据点。

Step7：独立成分分析：通过传感器获得的N个观测信号x_i(k),i＝1,2,···,N构成向量X(k)＝[x₁(k),x₂(k),···x_N(k)]^T，每个观测信号又由M个独立未知源信号s_i(k),i＝1,2,···M线性混合，且其构成向量S(k)＝[s₁(k),s₂(k),···s_M(k)]^T，则：

X(k)＝AS(k)，N≥M

其中，A为一N×M维的混合矩阵，为保证该问题有解，观测信号个数N必须大于或等于源信号个数M。最终是希望寻找到一个能够使源信号S(k)从观测信号X(k)分离出来的分离矩阵W，即：

Y(k)＝WX(k)＝WAS(k)＝CS(k)

C＝WA

其中，Y为问题求解后的成分信号，C为分离-混合的复合矩阵，不断优化W使得Y以最大程度去逼近S。

以联合熵H(y)作为目标函数：

H(y)＝H(y₁)+H(y₂)+···H(y_n)

其中，y_n为非线性输出，使H(y)最大，除各输出分量相互独立外。同时边缘熵H(y_i)的值也要达到最大，

H(y)＝-E{lnp_y(y)}

其中，p_y(y)是输出量y的概率密度函数。又由于y＝g(u)＝g(Wx),因此：

其中，p_x(x)为输入量x的概率密度函数，u_i无限逼近真实源信号。

使用梯度下降法对H(y)进行优化，即：

其中，T表示转置，x为观测信号，且

其中，g_i(·),i＝1,2,···,n为单调且非线性映射函数。

因此：

其中，η为学习率，通常为一个小于1的正数。

引入了自然梯度法对其进行优化，得到分离矩阵W的迭代公式：

其中，I为单位矩阵，W_(k)为迭代分离矩阵。

非线性函数g(·)使用的是Sigmoid函数：

最终得到分离矩阵W，进而得到源信号u＝Wx，从而对伪迹进行去除。

步骤3：使用傅里叶变换来将时域信号转换成频域信号，以获取感兴趣的频段(包括delta：频率范围为1-3Hz和theta：频率范围为4-7Hz频段)，即：

其中，w为频率，t为时间，e^-iwt为复变函数。F(ω)为象函数，f(t)为象原函数。

步骤4：提取基于相位滞后指数(Phase Lag Index，PLI)和加权相位滞后指数(Weighted Phase Lag Index，WPLI)的特征矩阵，使用PLI及WPLI作为衡量信号间耦合强度的指标。

PLI是一种能够有效最小化共同源问题(容积传导或活跃参考电极)的用于评估大脑不同区域间耦合性强度的度量，通过量化两个信号的相位差分布来估计其相位同步性，且分布通常表现为基于0的相位差，两个时间序列不存在相位耦合关系或是由共同源引起的相位同步(相位差接近或者等于0modΠ)，则相位差呈对称分布，若分布具有不对称性(相位差在范围(-Π,0)和(0,Π)的概率不同)，则两个时间序列间存在一个持续且非零的相位差，即相位滞后，与PLI不同的是，WPLI在PLI的基础上引入了权重的概念，即对于检测到的相位超前与滞后的贡献程度均由交叉谱的虚部大小进行加权，从而减少对其他无关噪声的敏感性，同时，能够提高其检测相位同步性变化的能力。

(1)PLI特征矩阵提取

假定两个EEG信号x_i(t)和x_j(t)在t时刻的相位差

满足：

其中，n、m均为整数，且在神经科学的研究中，通常取n＝1，m＝1，

和

为两个时间序列i和j的相位，当相位差接近于一个常数const时，则称i和j是相位同步的。信号x(t)的瞬时相位为：

其中，

是x(t)的希尔伯特变换：

其中，PV为柯西主值，能够避免t＝ξ时的奇异值。Π为圆周率，ξ为积分变量。

定义信号同步性指标，相位滞后指数PLI为：

其中，sgn为符号函数，N_t为时间序列分段个数，

为两个信号的瞬时相位差，PLI值的变化范围介于0到1之间，当PLI为0时，则相位差呈对称分布，即两个时间序列无相位同步，当PLI为1时，则相位差呈非对称分布，即两个时间序列完美相位同步，且相位同步越强烈，PLI就越接近于1。

(2)WPLI特征矩阵提取

其中，X为两个时间序列的交叉谱，ξ(X)则代表交叉谱的虚部分量，E{·}为数学期望，仅基于虚部分量的计算，增强了其对噪声的鲁棒性，且WPLI值通常介于0和1之间。

步骤5：脑功能网络连接

基于图论的思想，以62个通道作为网络的节点，分别以提取的两种特征矩阵(62×62)度量节点间边的连接程度，在稀疏度为15％-45％，步长为5％的阈值下分别构建对应的脑网络，脑网络需满足小世界性。一个网络是否是小世界网络通常是由属性值σ来衡量，由此涉及到两项参数值：γ和λ，定义如下：

其中，C_real和L_real分别为实际网络对应的聚类系数和最短路径长度，C_random和L_random分别为随机网络对应的聚类系数和最短路径长度。随机网络的构造与实际网络具有相同的点数和边数。当σ值大于1时，则说明该网络具有小世界属性，且σ值越大，小世界属性越强。

聚类系数通常是指网络中与某节点i的相邻节点间彼此又为邻居的可能性，能够衡量一个网络的局部连通性与集群特性，以下为其计算公式：

其中，k_i为与节点i相邻的节点数目，e_i为k_i个节点间实际存在的连接边的数目，k_i(k_i-1)/2则为其最多可能存在的连接边的数目，C_i表示聚类系数。

网络的平均聚类系数为N_c个节点的聚类系数的平均值，定义如下：

其中，C为聚类系数，V为网络中的节点集合。

最短路径长度能够反映网络内部进行信息传递的效率，定义如下：

其中，N_c为网络中所有节点的个数，l_ij为节点i与节点j之间的最短路径长度，L则为任意两节点间最短路径长度的平均值。

本实施例使用BrainNet Viewer，传入节点文件及边文件。其中，节点文件是一个“node”后缀的ASCII文本文件，主要描述节点的三维坐标，边文件是一个以“edge”为后缀的ASCII文本文件，表示节点之间的一个关联矩阵，该关联矩阵通过计算PLI与WPLI得出，且使用的分别为睁眼状态delta频段下以及闭眼状态theta频段下稀疏度为20％时的功能连接矩阵。

步骤6：根据前面所构建的脑功能网络，分别计算躯体维度和情感维度基于PLI的睁眼状态delta频段下稀疏度为20％时脑功能网络以及基于WPLI的闭眼状态theta频段下稀疏度为20％时脑功能网络的全局效率，参见图3和图4。计算公式如下：

其中，E_global为全局效率。

基于两个指标所得的全局效率进行融合，即为神经症下躯体维度和情感维度的特征向量。进一步可以使用决策树、K近邻或支持向量机进行分类。三种分类器进行分类的AUC面积，如图5和图6所示。

将62电极所在位置映射到大脑3D模型上，以定位躯体维度和情感维度脑功能网络的异常连接，参见图8和图9。其中，脑功能网络分别为基于PLI的睁眼状态delta频段下以及基于WPLI的闭眼状态theta频段下稀疏度为20％时的脑功能网络。躯体维度主要体现为各个脑区与顶叶、颞顶叶及枕叶的连接异常增强，而情感维度主要在额叶、前额叶与中央顶叶的连接异常增强。

通过以上步骤，本方法可以以较高的准确率对神经症下躯体维度和情感维度进行识别，且其脑网络的异常连接定位具有较高的可靠性。如图7所示，验证了本特征提取方法的有效性。

以上方案内容仅表述本发明的技术方案，并非最完美精确的解决方案。随着技术的革新与时代的变迁，方案可能产生更合理更高效的变化。选择示例性实施方式并进行描述是为了解释本发明的原理和应用，为了便于研究人员与技术人员的参考、理解与实践本发明的具体细节。若在本发明构想的基础上所作的改变，其所产生的功能作用仍未超出说明书及附图所涵盖的内容时，仍应属本发明的保护范围。

参考文献

[1]Yu Q.,Y.Du,J.Chen,et al.Application of Graph Theory to AssessStatic and Dynamic Brain Connectivity:Approaches for Building Brain Graphs[J].Proceedings of the IEEE,2018,106(5):886-906.

[2]Gonuguntla V.,K.C-Veluvolu,J.Kim.Recognition of Event-associatedBrain Functional Networks in EEG for Brain Network Based Applications[A]//2020:271-274.

[3]Al-Momani S.,S.Dhou.Spinal functional Magnetic Resonance Imaging(fMRI)on Human Studies:A Literature Review[A]//2019:1-5.

[4]S-Takkar M.,M.Kumar-Sharma,R.Pal.A review on evolution of acousticnoise reduction in MRI[A]//2017:235-240.

[5]Rasheed K.,A.Qayyum,J.Qadir,et al.Machine Learning for PredictingEpileptic Seizures Using EEG Signals:A Review[J].IEEE Reviews in BiomedicalEngineering,2020,1.

[6]F.Karimzadeh et al.,“A distributed classification procedure forautomatic sleep stage scoring based on in_x0002_stantaneouselectroencephalogram phase and envelope features,”IEEE Transactions on NeuralSystems and Rehabilitation Engineering,2018.

[7]S.S.Mohd Radzi,V.S.Asirvadam and M.Z.Yusoff,"Fractal Dimension andPower Spectrum of Electroencephalography Signals of Sleep Inertia State,"inIEEE Access,vol.7,pp.185879-185892,2019,doi:10.1109/ACCESS.2019.2960852.

Claims

1.一种基于图论和机器学习的神经症特征提取方法，其特征在于，包括以下步骤：

(1)采集EEG数据，对其进行预处理；

(2)针对步骤(1)预处理后的数据，使用傅里叶变换将时域信号转换成频域信号，以获取感兴趣的频段；

(3)提取基于相位滞后指数(PLI)和加权相位滞后指数(WPLI)的特征矩阵；

2.根据权利要求1所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：步骤(1)所述预处理包括剔除无用电极、重参考、滤波、剔除坏段、分段和基线校正、降低采样率以及基于独立成分分析的伪迹去除。

3.根据权利要求2所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：所述基于独立成分分析的伪迹去除包括以下步骤：

以寻找一个能够使源信号S(k)从观测信号X(k)分离出来的分离矩阵W为目的，以联合熵H(y)作为目标函数，使H(y)最大：

H(y)＝-E{lnp_y(y)}

其中，p_y(y)是输出量y的概率密度函数，

其中，p_x(x)为输入量x的概率密度函数，u_i为无限逼近于真实的源信号，n为源信号个数。

使用梯度下降法对H(y)进行优化，即：

其中，T表示转置，x为观测信号，且

其中，g_i(·),i＝1,2,···,n为单调且非线性映射函数；

因此：

其中，η为学习率；

其中，I为单位矩阵，W_(k)为迭代分离矩阵，

4.根据权利要求1所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：所述感兴趣的频段包括delta：频率范围为1-3Hz和theta：频率范围为4-7Hz。

5.根据权利要求1所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：所述提取基于相位滞后指数特征矩阵包括以下步骤：

两个EEG信号x_i(t)和x_j(t)在t时刻的相位差

满足：

其中，n、m均为整数，φ_i(t)和φ_j(t)分别为两个时间序列i和j的相位，当相位差接近于一个常数const时，则称i和j是相位同步的，信号x(t)的瞬时相位为：

其中，

是x(t)的希尔伯特变换：

其中，PV为柯西主值，Π为圆周率，ξ为积分变量；

定义相位滞后指数PLI为：

其中，sgn为符号函数，N_t为时间序列分段个数，

为两个信号的瞬时相位差，PLI值的变化范围介于0到1之间，当PLI为0时，则两个时间序列无相位同步，当PLI为1时，则两个时间序列完美相位同步，且相位同步越强烈，PLI就越接近于1。

6.根据权利要求1所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：所述提取加权相位滞后指数特征矩阵包括：

其中，X为两个时间序列的交叉谱，ξ(X)则代表交叉谱的虚部分量，E{·}为数学期望。

7.根据权利要求1所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：所述小世界性由属性值σ来衡量，当σ值大于1时，则说明该网络具有小世界属性，且σ值越大，小世界属性越强，定义：

其中，C_rdal和L_CCCC分别为实际网络对应的聚类系数和最短路径长度，C_random和L_random分别为随机网络对应的聚类系数和最短路径长度；

聚类系数指网络中与某节点i的相邻节点间彼此又为邻居的可能性，能够衡以下为其计算公式：

其中，k_i为与节点i相邻的节点数目，e_i为k_i个节点间实际存在的连接边的数目，k_i(k_i-1)/2则为其最多可能存在的连接边的数目，C_i表示聚类系数；

其中，N为网络中所有节点的个数，l_ij为节点i与节点j之间的最短路径长度，L则为任意两节点间最短路径长度的平均值。

8.根据权利要求1-7任一项所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：还包括根据步骤(3)所构建的脑功能网络，分别计算躯体维度和情感维度基于PLI的睁眼状态delta频段下稀疏度为20％时脑功能网络以及基于WPLI的闭眼状态theta频段下稀疏度为20％时脑功能网络的全局效率，基于两个指标所得的全局效率进行融合，即为神经症下躯体维度和情感维度的特征向量。

9.根据权利要求8所述一种基于图论和机器学习的神经症特征提取方法，其特征在于：还包括将62电极所在位置映射到大脑3D模型上的步骤。