曹春雨 韩 凌 李春胜
(沈阳工业大学电气工程学院生物医学工程系,沈阳 110870)
癫痫是一种常见神经系统疾病,全球约有5000多万患者。脑电图是癫痫诊断的最重要工具。癫痫发作短暂,难以预测,长时程视频脑电图常被用于癫痫发作的记录;与此同时,带来数天甚至数周的连续记录,数据量大,分析耗时耗力。基于以上因素,癫痫脑电信号的自动分类研究,尤其是机器学习方法,是大幅提高分析效率和准确性的重要途径。
脑电图按频率可分为不同节律,低频节律通常表现在不同脑区间的协调与交互,高频节律更多地包含着局部活动信息。通常,不同频率神经网络振荡间能相互影响,不同频率之间的相互作用称为跨频率耦合[1]。相位-幅度耦合(phase amplitude coupling, PAC)是跨频率耦合的一种类型,即高频节律的幅度被低频节律的相位调制。研究表明,高频与低频信号间协调工作是完成大脑功能的重要基础,如低频θ与高频γ间的相位耦合构成了记忆过程基础[2]。还有研究表明,癫痫发作过程伴随着低频节律相位和高频节律幅度间耦合。通过颅内电极获得的脑电信号中低频相位调制高频幅度的特征,可获得很好的感兴趣区和手术切除区的一致性[3]。在头皮脑电图中,相位-幅度调制信息可以定位到致痫区[4]。在癫痫发作区,高频(γ和Ripple)与低频(δ,θ,α和β)节律的跨频率耦合明显高于正常组织[5]。
癫痫脑电信号的分析主要有时域、频域、时频域和非线性动力学等几个方面[6-10]。国内外很多学者对癫痫脑电自动分类进行了深入研究,取得了很多进展[11-13]。本研究将相位-幅度调制的方法应用到癫痫发作状态的自动分类中,对癫痫脑电信号进行了高低频节律间耦合研究,结果显示,发作期脑电图与发作间期相比,多个节律间存在显著的相位-幅度调制,利用量化的调制指数特征对发作期和发作间期数据分类的准确率可达97%,是区分癫痫发作间期和发作期的一种有效特征。
本研究采用的数据来自德国波恩癫痫研究室临床采集的短程脑电数据库,该数据库的数据被广泛应用于癫痫脑电信号分析与癫痫脑电的识别研究中[14]。在本实验中,采用癫痫患者致痫区的发作间期(F数据集)和发作期(S数据集)的脑电信号。每个数据集包含有100个单通道、时长23.6 s的数据,即F和S分别有100段数据,样本量为200个,信号采样频率为173.61 Hz。
癫痫脑电信号幅度和相位信息提取,可以在滤波后进行希尔伯特变换,或者进行复数小波变换。对癫痫脑电信号进行低频段滤波得到低频信号,再对低频信号进行希尔伯特变换可得相位信息;同理,将脑电信号进行高频段滤波得到高频信号,再对高频信号进行希尔伯特变换得到高频幅值[15]。图1是癫痫发作期S数据集的第81通道数据(S81)的相位幅值的提取过程,其中(a)是原始信号,(b)是从(a)中黄色标记区域截取并放大的部分,(c)是对截取信号进行FIR滤波,频率范围为(30~70)Hz,其包络线是滤波后信号进行希尔伯特变换得到的幅值,(d)是低频信号(2~4)Hz进行希尔伯特变换得到的相位信息,(e)是(-π,π]区间的相位。
图1 癫痫脑电数据低频相位和高频幅值的提取。(a)原始信号;(b)截取的信号; (c)高频节律及幅度; (d)低频节律; (e)低频节律的相位Fig.1 Example of amplitude and phase of high and low frequency in epileptic EEG. (a) The original signal; (b) Selected signal; (c) High frequency rhythm and amplitude; (d) Low frequency rhythm; (e) Low frequency phase
为了进一步分析低频相位和高频幅值的关系,对图1(e)相位进行分段,将整个相位范围均分成18个区间,每个区间π/9弧度。如图2所示,对每个相位区间下所有时间点的幅值取均值,并对其均值进行归一化,最终得到相位幅值关系,在相位π/9~3π/9之间取得最大值,红色实线是高斯拟合结果(μ=0.13,σ=0.36)。
图2 相位幅度关系Fig.2 Phase-amplitude correlation
提取相位幅度的另一种方法是采用复数小波变换(complex wavelet transform, CWT),它是连续小波变换的一种特殊形式,复数小波系数由Morlet母小波在带宽为5 Hz和中心频率为0.8 125 Hz时获得。
如果低频率的相位和高频率的幅度之间没有耦合,则在相位幅值关系上每一小段相位所对应的幅值是相等的(相位幅度分布P和均匀分布一致),若每一小段相位对应的幅度有差异,则存在耦合。基于这个原因,定义一种方法去量化相位幅度分布P与均匀分布U之间的差异,这种方法是通过计算Kullback-Leibler (KL)距离来确定差异的,KL距离是一种广泛应用在统计学和信息理论,目的是来推断两种分布之间差异的数量[16]。而调制指数(modulation index, MI)是一个常数乘以相位幅度分布P和均匀分布U之间的KL距离,其值在0~1之间,它的具体计算方法如下。
相位幅度分布P和分布Q的KL距离的定义为
(1)
式中,DKL(P,Q)≥0,当P=Q时,DKL(P,Q)=0,j代表第j个相位分段。
相位幅度分布P的KL距离的公式格式类似于香农熵(H)的定义,有
(2)
KL距离和香农熵H相关,关系如下:
DKLP,U=logN-H(P)
(3)
式中,U代表均匀分布,logN是最大可能熵值,N是相位分段的个数。
利用KL距离定义的MI值如下:
(4)
支持向量机(support vector machine, SVM)是根据统计学习理论提出的一种机器学习方法。它通过非线性映射将向量映射到高维的特征空间,通过选择最优分类面得到一个超平面分割,该分类面能将两类模式分开,并保证其间隔最大化[17]。
设分类面方程为ωx+b=0,则在k维空间中,使得样本(xi,yi),(i=1,…,n),x∈Rd,y∈{+1,-1},满足
(5)
s.t.yi(ωxi+b)-1≥0 (i=1,…,n)
在线性不可分的情况下,可引入ξi≥0,式(5)变为
(6)
s.t.yi(ωxi+b)-1+ξi≥0 (i=1,…,n)
式中,C为误差惩罚因子,用于平衡松弛变量ξi和分类边界的参数,即综合考虑最少错分样本和最大分类间隔。
利用Lagrange乘子方法解决最优问题后,最优决策函数为
(7)
式中,K是SVM的核函数,ai、b分别表示支持向量系数和分类阈值。
径向基是应用最为广泛的非线性核函数,它有两个参数:惩罚因子C和核函数参数g,其表达式为
(8)
本研究提取了发作间期和发作期数据的相位幅值耦合特征。有研究认为,δ节律(0.5~4 Hz)中的高频部分δH(2~4 Hz)与γ节律(30~80 Hz)存在耦合关系,与癫痫发作关系密切[18-19]。本研究将低频分为3个频段,分别是δL(0.5~2 Hz)、δH(2~4 Hz)和θ(4~8 Hz),高频分为3个频段,分别是α(8~13 Hz)、β(13~30 Hz)和γ(30~80 Hz),得到9组频率的组合,如表1所示。对每一对频率进行相位幅度的提取,并得到相位幅度关系,计算MI值,对每一个通道数据的MI值平均得到一个特征,因有9组频率组合,每个通道数据对应9个特征。差异性分析采用单因素方差分析(ANOVA)[20]。
表1 低频和高频节律划分Tab.1 Range of low and high frequency rhythms
进一步计算不同特征的受试者工作特征曲线(receiver operating characteristic curve, ROC),对曲线下面积(area under the curve, AUC)进行分析。
为了分析不同分类方法对结果的影响,分别采用支持向量机线性核函数、非线性核函数和随机森林,对发作期和发作间期的脑电数据MI特征进行分类。随机选取样本的50%作为训练集、50%作为测试集,采用五折交叉验证,得到随着惩罚因子C变化的训练样本准确率,优化C取值,然后计算测试集分类准确率。为了与近似熵和样本熵特征比较,对6个频段对应特征计算ROC曲线、AUC值和分类准确率。
对S81进行CWT变换,得到相位和幅度信息,图3(a)为高频(30~70 Hz)幅度信息,图3(b)为低频(2~4 Hz)相位信息,分别与图1(c)、(e)中的幅值和相位信息相对应。
图3 复数小波变换的幅值信息和相位信息。(a)幅度信息;(b)相位信息Fig.3 Amplitude and phase obtained by complex wavelet transform. (a) Amplitude information; (b) Phase information
数据S81的MI如图4所示,蓝色区域为低MI值的区域,由蓝色区域到红色区域,调制指数在逐渐增强,此区域所对应的低频和高频范围可以反映出相位幅值的分布不是均匀的,存在相位和幅度的耦合。
图4 S81发作期的MI图Fig.4 MI map of S81 during seizure
图5 发作期、发作间期统计分析Fig.5 The statistical analysis results of F and S dataset
F和S数据集的MI值如图5所示,绿色柱状图和黄色柱状图分别代表发作间期和发作期。结果显示,发作期δH-γ的MI值(0.009 9±0.009 6)相比发作间期(0.003 6±0.008 7)显著增加(P<0.01)。θ-γ的发作期MI值(0.008 7±0.006 2)相比发作间期(0.001 4±0.003 2)也有显著增加(P<0.01);θ-β耦合强度在发作期(0.002 2±0.001 3)与发作间期(0.000 5±0.000 7)也存在显著差异,如图5星号标注,而δL-α、δL-β和δL-γ的差异不明显(P>0.01)。
图6是发作期与发作间期各频段幅频调制特征的ROC曲线,支持向量机线性核函数惩罚因子lbC=16。结果显示,用θ-β、θ-γ和δH-γ三组频率所得的MI值求得的特征有较好的准确率,其对应AUC值分别为0.945、0.948和0.873。
图6 发作期与发作间期各频段幅频调制特征的ROC曲线Fig.6 The ROC curve of phase-amplitude modulation in each frequency band
支持向量机线性核函数分类结果如图7所示,lbC=16,准确率为97%。
图7 交叉验证下的训练集和测试集的分类准确率Fig.7 The classification accuracy of training and testing set through cross validation
在非线性核函数条件下,采用网格搜索方式确定最大准确率的C和g值,lbC=4,lbg=11,如图8所示,分类准确率为95%。
图8 径向基核函数对训练集的参数寻优过程。(a)非局部参数下的分类准确率;(b)局部参数下分类准确率Fig.8 Radial basis kernel function parameter optimization of the training set.(a) Classification accuracy under non-local parameters;(b) Classification accuracy under local parameters
随机森林方法的分类结果如图9所示,在决策树数量为16时,准确率为97%。
采用线性核SVM、非线性SVM及随机森林3种方法,使用MI特征的分类准确率分别为97%、95%和97%,MI特征相对不同方法的分类结果差异不明显。与非线性SVM参数的网格搜索相比,随机森林和SVM线性核参数选择简单,更适合应用。当前数据集结果显示SVM线性核分类结果略优于非线性核。
利用6个频段的样本熵特征得到的ROC曲线如图10所示。结果显示,θ、β和γ节律对应的AUC值分别为0.908 2、0.902 8、0.957 5。进一步利用线性核函数对样本分类的准确率为92%,其中lbC=16,如图11所示。
在采用近似熵提取特征时,分类结果略有下降,AUC较大值的分别为0.812、0.789、0.819,对应δH、α和β节律,整体分类准确率为89%。
图11 交叉验证下样本熵的训练集和测试集的分类准确率Fig.11 The classification accuracy of sample entropy training and testing set through cross validation
γ节律被认为与癫痫发作密切相关[21]。本研究显示,γ节律与不同低频节律的相位幅度耦合在发作期都有增强,如图5所示。结果表明,在癫痫发作期γ节律并不是单一出现,其产生与低频信号密切相关。除低频δ与α、β和γ的耦合无显著变化外,其他高低频节律耦合性均有增强,如图5所示。高低频节律相位幅度调制是发作期脑电图的普遍形式,可能反映了局部癫痫网络和全局网络间的调节过程。本研究根据文献[19]将δ分为δL和δH,结果显示二者在发作期和发作间期差异明显,与文献[19]中描述的结果相似。脑电高频振荡(high frequency oscillation, HFO)被认为与致痫区密切相关,而HFO与低频节律的相位耦合可用于定位致痫区。本研究由于所使用数据的采样率较低,对HFO与低频节律耦合的分类情况还需进一步研究。
图6的ROC曲线显示θ-β和θ-γ的AUC值分别为0.945和0.948,表明对应的MI值能对发作期和间期进行较好区分。文献[2]显示,θ-γ相位-相位耦合与记忆认知活动相关,而在癫痫发作期,其相位-幅度耦合有显著增强(见图5),说明高低频节律间这两种耦合形式可能分别与生理和病理过程相关。θ-β的MI值在两种状态下有显著差异,未见到其他研究报道。由于β节律与肢体运动功能相关[22],推测θ-β耦合可能与发作期患者肢体异常活动相关。
文献[11]用样本熵和自回归模型提取癫痫发作期和发作间期的脑电信号特征,支持向量机的分类结果为88.25%。仅利用近似熵作为特征,对癫痫发作期和发作间期的分类结果为88.25%[13]。在对6个频段分解的基础上,本研究利用近似熵和样本熵得到的准确率分别为89%和92%,较上述文献有所提高。从图10可以发现,γ节律样本熵特征的AUC值大于θ-β相位-幅度耦合方法提取的MI特征的AUC值,但所得到的分类准确率低于MI特征的分类准确率,其原因可能是MI特征在单一频率特征的基础上还包含了多个频率之间的相互作用,更能够反映癫痫发作的不同状态。本方法的分类准确率高于单一特征提取方法的近似熵和样本熵。基于以上结果,还可以将MI特征与其他特征相结合来进一步提高准确率。
在本文中,介绍了两种不同的计算相位和幅度的方法,与复数小波方法相比,滤波方法容易产生谐波分量,不正确滤波范围则会破坏原节律的完整性,影响结果。在FIR滤波方法之外,可以采用其他的信号分解方法(如经验模态分解),提高计算结果可靠性。
脑电图高频与低频节律间协同工作是大脑完成功能的重要基础,病理状态下高低频节律之间也存在显著耦合关系。本研究对癫痫患者脑电图分析显示:与发作间期相比,发作期γ节律与不同低频节律的调制指数均显著增强;发作期θ-β相位幅度耦合显著增强;癫痫脑电信号的MI特征可以准确区分发作期和发作间期数据,且MI特征相对不同分类方法的结果差异不明显,稳定性较好,可获得较高的分类准确率。本方法的应用可提高临床视频脑电图的分析效率,同时为研究基于癫痫脑电低频节律提取病理性高频信号提供依据。
[1] Jensen O, Colgin LL. Cross-frequency coupling between neuronal oscillations [J]. Trends Cogn Sci, 2007,11(7): 267-269.
[2] Belluscio MA, Mizuseki K,Schmidt R, et al. Across frequency phase-phase coupling between theta and gamma oscillations in the hippocampus [J]. J Neuroscience, 2012, 32(2): 423-435.
[3] Guirgis M, Chinvarun Y, Del Campo M, et al. Defining regions of interest using cross-frequency coupling in extratemporal lobe epilepsy patients [J]. J Neural Eng, 2015, 26(11): 1-19.
[4] Li Chunsheng, Jacobs D, Hilton T, et al. Epileptogenic source imaging using cross-frequency coupled signals from scalp EEG [J]. IEEE Trans Biomed Eng, 2016, 63(12): 2607-2618.
[5] Amri M, Frauscher B, Gotman J. Phase-Amplitude coupling is elevated in deep sleep and in the onset zone of focal epileptic seizures [J]. Front Hum Neurosci, 2016, 10: 1-11.
[6] Ubeyli ED. Combined neural network model employing wavelet coefficients for EEG signals classification [J]. Digital Signal Processing, 2009, 19(2): 297-308.
[7] 宦飞,郑崇勋,黄远桂.基于时频分布的EEG中棘波的相关检测[J]. 仪器仪表学报, 1999,20(5): 446-450.
[8] 李春胜,王宏.脑电Alpha波的混沌关联维数分析[J]. 仪器仪表学报, 2009, 30(3): 477-480.
[9] Mohseni HR, Maghsoudi A, Kadbi MH, et al. Automatic detection of epileptic seizure using time-frequency distributions[C]//IET 3rd International Conference MEDSIP 2006. IEEE, 2006:29.
[10] 崔林阳,田心.基于颞叶癫痫脑电相位互相关分析法的癫痫超同步放电的初步研究[J]. 中华神经医学杂志, 2008,7(3): 219-222.
[11] 王杰,李牧潇.改进的极限学习机在癫痫脑电分类中的应用[J].计算机仿真, 2014, 31(6): 343-347.
[12] 张涛,陈万忠,李明阳. 基于频率切片小波变换和支持向量机的癫痫脑电信号自动检测 [J]. 物理学报, 2016, 65(3): 038703.
[13] 袁琦,周卫东,李淑芳,等.基于ELM和近似熵的脑电信号检测方法[J]. 仪器仪表学报, 2012, 33(3): 514-519.
[14] Andrzejak RG, Lehnertz K, Mormann F, et al. Indications of nonlinear deterministic and finite dimensional structures in time series of brain electrical activity: Dependence on recording region and brain state [J]. Physical Review E, 2001, 64: 061907.
[15] Tort AB, Komorowski R, Eichenbaum H, et al. Measuring phase-amplitude coupling between neuronal oscillations of different frequencies [J]. J Neurophysiol, 2010, 104(2): 1195-1210.
[16] Kullback S. The Kullback-Leibler distance [J]. Am Statist Assoc, 1987, 41: 340-341.
[17] Chang CC, Lin CJ. A library for support vector machines [J]. ACMTIST, 2011, 2(3):1-27.
[18] Canolty RT, Edwards E, Dalal SS, et al. High gamma power is phase-locked to theta oscillations in human neocortex [J]. Science, 2006, 313(5793): 1626-1628.
[19] Nonoda Y, Miyakoshi M, Ojeda A, et al. Interictal high-frequency oscillations generated by seizure onset and eloquent areas may be differentially coupled with different slow waves [J]. Clin Neurophysiol, 2016, 127(6): 2489-2499.
[20] 李红利.癫痫脑电信号的非线性分析 [D]. 天津: 天津大学, 2012.
[21] Ren L, Kucewicz MT, Cimbalnik J, et al. Gamma oscillations precede interictalepileptiform spikes in the seizure onset zone [J]. Neurology, 2015, 84(6): 602-608.
[22] Brown P, Oliviero A, Mazzine P, et al. Dopamine dependency of oscillations between subthalamic nucleus and pallidum in parkinson′s disease [J]. J Neurosci, 2001, 21(3): 1033-1038.