沈潇童,毕 卉,王苏弘,李文杰,邹 凌
1.常州大学 信息科学与工程学院,江苏 常州213164
2.常州市生物医学信息技术重点实验室,江苏 常州213164
3.常州市第一人民医院,江苏 常州213003
抑郁症(Depressive Disorder)是一种常见的精神疾病,其特征是对自身喜欢的活动失去兴趣,并伴有至少两周无法进行日常活动的行为。近年来青少年抑郁越来越引起社会和家庭的关注[1]。通过计算机辅助技术对抑郁症患者和健康人进行分类诊断,有助于降低临床诊断的误诊率,帮助患者及早得到正确的治疗。
目前,脑电在抑郁症的研究领域已经有了很多重要发现。中国科学院大学的盖淑萍等人提取了静息态下的脑电信号的小波包节点功率谱熵特征,发现其能够作为鉴别抑郁症的有效量化指标[2]。Zhu 等根据静息状态功能磁共振成像研究重度抑郁症患者自发性脑部活动的变化,发现了双侧正中前额叶皮质、楔前叶、角回、右侧海马旁回、右侧颞极自发脑区活动的改变[3]。Dharmadhikari 等人在静息态θ 波段发现了前额叶区域的不对称性,证明前额叶的θ 不对称是一种潜在性的生物标志[4]。此外,Esther 在对成年人抑郁患者的研究中发现,静息状态下抑郁患者产生的心理伤害在脑电δ 波段的左额叶和右额叶存在联系[5]。Lee 等在研究中发现,高α 波段(10~12 Hz)功率能够对抑郁组和正常组进行区分,受试者工作曲线值达到了70%[6]。Mumtaz等使用静息态脑电数据,提取额叶α 不对称和信号频段功率作为特征,发现在顶区、中心区域、颞区、枕区抑郁患者的右半脑特征比左半脑的特征显著[7]。Bachmann 等通过静息态脑电图(Electroencephalogram,EEG)对中年女性抑郁患者进行了分类研究,并通过对单导联电极Pz的特征提取和分类研究,提出单导联检测抑郁的可行性[8]。
对于导联的选择,本质是对导联相应的信号提取的特征进行选择的过程。Wang等人利用共空间模式算法(Common Spatial Pattern,CSP),通过空间模式向量和导联的最大系数的关系对电极进行优化选择[9]。Arvaneh等在此基础上提出了稀疏共空间模式(Sparse Common Spatial Pattern,SCSP)算法,其是在分类精度约束下对最小信道数进行优化,对去除噪声和不相关信道后的结果进行电极优化选择[10]。此外,Li等使用主成分分析法对每个被试的电极对应的数据集进行计算,然后根据不同的任务刺激的发生率进行导联的选择[11]。
机器学习在脑电分类学习中有着广泛的运用,利用决策树和SVM对静息状态脑电信号的波段功率进行分类的精度达到了80%[12],基于此分类方法能够对单导联抑郁症患者和正常人的分类提供有效的支持。Alaa 等人使用了极限学习机(Extreme Learning Machine,ELM)算法对抑郁患者和健康人的脑电信号进行分类,达到了92%的分类精度[13]。Xia 团队使用静息态的脑电数据,使用小波变换、短时傅里叶变换、经验模式分解提取特征,使用逻辑回归算法进行分类,发现小波变换提取的特征精度最高达到了87.5%[14]。
基于EGI 公司64 导脑电采集系统,本文首先采集了16名抑郁患者和16名健康人的静息状态下的脑电信号。然后,通过对全导联的特征提取和特征分析,确定特征显著的脑区,并确定最佳导联,同时使用遗传算法进行整体导联的选择。最后,对单导联和多导联提取特征进行分类,并对时域和频域算法的运行效率进行分析。
频谱不对称分析法是利用功率谱密度计算单个导联信号能量比率的算法。它通过计算α 频段附近高频段和低频段相对性差异来描述脑电信号的频谱不对称性[15]。
计算方法分成三部分,参数如图1所示[8]。
图1 频域波段分布
(1)计算高低频段边界F1n、F2n、F3n、F4n
其中,n 代表人数(n ∈[1,32]),m 代表导联数(m ∈[1,64]),F1、F2 代表低频段,F3、F4 代表高频段,smn为功率谱密度,Slmn和Shmn分别为高低功率谱密度,fcn为α 波段的功率最大值。
这里F1、F2 的确定在靠近θ 波段处,两者频带宽度相差4 Hz,F3、F4 频段的确定在靠近传统意义上的β 段处,两者频带宽度相差24 Hz。
(2)计算高低频率功率谱密度总和
(3)计算选择高低频段的SASI值(即为SASI的特征)
去趋势波动分析算法最早运用于DNA 序列的检测[16],现在逐渐运用于脑电信号分析领域,其能够直接计算时域上信号的特点[17]。为了进行DFA 的均方根系数计算,将4分钟的EEG信号分成均等的段数,每段20 s,使用DFA算法进行处理。计算流程通常分为三个步骤。
(1)计算拟合序列
对于一个时间序列的信号,定义序列x(i),i 为每段的长度,区间为[1,N];通过求离散序列的积分,得到一个新的时间序列。如式(8)所示。
(2)计算均方根系数F(n)
在确定拟合前后的序列后,通过式(9)计算出结果F(n)。
(3)计算F(n)和n 的关系
同时取对数ln F(n)~ln n,得到两者之间的斜率α,α 即为提取的时域特征。
目前,遗传算法(Genetic Algorithm,GA)广泛被运用于求最优解问题和特征选择中。以染色体X 为例,X={0,1,0,1,1,0,0,0,0,0},染色体X 由10个基因以二进制形式组成,1 位表示特征选中,0 位表示特征未选中。在特征选择中,将基因的数量类比为特征的数量。
遗传算法流程分为五部分,分别为初始化、适应度函数、选择、交叉、变异。具体算法伪代码如下所示。
步骤1 初始化参数X、Z、N、Tmax、CN、TC ,设置X 为初始化染色体群体,Z 为最佳的染色体,N 为最优选择染色体数,Tmax为一代的最大数量,CN 为交叉的总数,TC 为交叉次数的两倍;
步骤2 估计染色体的适应度函数F(X);
步骤3 For t=1 to Tmax
For i=1 to CN
使用轮盘选择方法选择两对父母染色体
通过两对父母染色体交叉产生两个子代
Next i
For j=1 to TC
根据突变率进行变异
Next j
估计产生子代的适应度
添加新产生的子代到当前种群X 中
对种群进行排序,选择最优的N 条染色体
如果种群中有更好的染色体,更新Z
Next t
所有被试均在常州市第一人民医院经临床确诊,共筛选出患有抑郁病症的青少年16例,健康青少年16例,年龄范围(16.31±1.25),抑郁组和对照组之间的年龄无统计学意义(p >0.05)。所有被试实验前均使用汉密顿抑郁量表(Hamilton Depression Rating Scales,HADS)进行检测[18],抑郁组分数高于对照组。
两组的青少年实验研究已经通过常州市第一人民医院伦理委员会批准。所有受试者或其监护人均签署了知情同意书,自愿参加本实验。
采用EGI公司64导联脑电采集系统,如图2所示(脑电采集系统帽以61、62号导联为前端,靠近眼部;以35、39、37 号导联为后端,靠近耳部),采样频率为500 Hz,电极阻抗设定在80 kΩ以下,利用Net Station软件实时采集。使用EEGLAB(版本号v14.1.2)工具箱预处理原始脑电信号,参考电极使用平均参考。
图2 64导联脑电采集系统
采集的静息态的数据分为闭眼数据和睁眼数据,闭眼条件下被试受到的外部视觉刺激较少,采集的信号受到噪声干扰较小,提升了信号的质量,减少了对特征提取的干扰。因此,实验采集了4分钟闭眼状态下的脑电数据,被试坐在安静的房间内,尽量保证头部不会大幅度移动,在数据采集过程保持清醒,范式流程图如图3所示。
图3 范式流程图
将采集的数据直接通过Net Station 软件转换成Matlab可读取的格式,在EEGLAB工具箱下对信号进行处理,具体步骤包括:低通和高通滤波(0.5~40 Hz)[8,15]、人工的伪迹检测和坏导联替换、参考点转换,然后使用独立成分分析(Independent Component Analysis,ICA)去除眼动以及其他噪声。最后选择干净的3 分钟内的数据,作为特征提取的数据段。
其中,由于β 频段包含对抑郁评估有用的电生理信息,而θ 频段是稳定不受疾病影响的频段[19],SASI 算法就是通过计算出这两个实际频段的功率,得到频段间功率的差异率。而波段的边界频率是专门为每个个体选择的,并与脑电信号α 频段的最大功率相关。理论上所选择的计算频带接近传统的β(13~30 Hz)和θ(4~8 Hz)频段,但当最大的功率值落在较高或者较低的α 频段(8~13 Hz),那么最高和最低的实际频率的数值存在偏移。为了保证范围,将滤波的频段定为0.5~40 Hz[15]。
首先,本文对抑郁症患者的脑区进行相关电生理机制的查询,发现在抑郁症患者的脑电信号生理机制研究中,抑郁症患者相比于正常人在顶叶、枕叶的区域有较多的改变[20],不对称的现象较为常见,这些改变与患者年龄和抑郁症程度有关[21]。
然后,根据SASI 和DFA 算法,提取了全部32 位被试64 导联的时频特征。根据式(7)计算SASI 的值如图4 所示,抑郁组的整体特征呈现正性相关,对照组则呈负性相关。根据式(9)计算的DFA 的值如图5 所示,发现抑郁组的特征上升趋势比对照组快,52~61导联的数值存在波动,能够发现31~34导联存在一个明显的下降趋势,而SASI的值在31~34更趋于稳定。
图4 SASI各导联均值折线图
图5 DFA各导联均值折线图
再根据抑郁患者在顶区和枕区发现的DFA和SASI特征,对各个导联进行结果统计,并和已发现的顶区(P区)和枕区(O 区)进行生理特性结合,分别对抑郁组和对照组顶区和枕区的导联提取的SASI和DFA的特征进行t 检验,根据两种特征t 检验结果,发现Pz 导联和Oz导联具有较高的统计学意义(P <0.05),如表1所示。而进一步比较发现Pz 导联比Oz 导联在统计学上更具有差异性,最终选择Pz导联进行特征分类。
表1 两导联t检验结果
针对提取32 位被试的64 导联的时频特征,分别构建了SAIS 和DFA 的32×64 维度的特征矩阵,利用遗传算法对每个人每个导联的特征进行选择、交叉、变异,最终将得到的子代的特征对应的导联标签号选出。选取的导联的编号如表2所示。
表2 两特征选取的导联
如表2所示,SASI提取了总共34导联的特征(P <0.05),DFA提取了26导联的特征(P <0.05)。然后,将两者同时选中的导联作为两者的分类导联,如图6 所示,以此减少了导联数量,提升了特征质量,便于后续的分类。
图6 导联选择结果
Pz导联提取的特征作为分类的特征集,构建了单导联的特征值矩阵,利用SVM建立了分类模型,使用留一法进行交叉验证和分类(参数:训练数据X ,标签Y ,核函数Gaussian,惩罚系数1,预期标准化true,核尺度5.8,类别[0;1],折数32)。分类结果如表3 所示。使用SVM分类器对Pz导联特征分类的结果来看,SASI分类的精度都达到了85%,DFA的分类精度达到了86%。
表3 SVM单导联分类结果 %
根据遗传算法选取的每一导联对应的特征,使用SVM对其进行分类,分类结果如表4所示。
表4 SVM多导联分类结果 %
为了能够更好地为脑机接口技术提供范例,本文基于Matlab 2016a 平台的Matlab 语言进行计算,统计了SASI 和DFA 算法的运行和时间,便于后续BCI 技术的参考。本文对于SASI 和DFA 的实时性进行了分析,分别统计了一分钟内10 s为间隔下,计算64导联的算法需要运行的时间(机器配置:Intel®i5-4590CPU,32 GB 内存,NVIDIA GeForce GTS 450 显卡)。DFA 算法计算时间如图7所示,SASI算法计算时间如图8所示。
图7 DFA算法耗时
图8 SASI算法耗时
从多导联选择的角度分析,遗传算法对SASI 方法选取了34导联,DFA方法选取了29导联,通过两者所选取导联的交叉比较分析,将3,4,12,19,20,26,27,30,33,34,45,48,58,62,共计14 个导联进行分类(橙色标记),如图4所示。图4中时频特征选择的导联在左侧较多,前额叶、额叶部分均有涉及,顶区包括了34 号的Pz以及附近的33 号导联,这些区域对后续抑郁症的研究具有参考意义。
经过特征和导联分类结果来看,多导联的分类精度比单导联的分类精度提升36%左右,DFA算法提升得更加显著,达到了90.6%,经过遗传算法提取的导联进行特征分类能对分类准确率进行改善。
从特征的性质来看,本文提取的DFA 特征是直接从时域的角度进行计算,将均方根系数与每个时间段的关系α 作为特征,α 的值在[0.5,1]代表了序列具有相关性,即大波动后可能出现大波动,小波动后可能出现小波动[17]。因此α 作为特征进行分类过程中能够更好地对抑郁症和正常人的时间序列进行评估。SASI计算了最高α 波段功率附近高低频段的频谱差异率作为特征,其正性和负性从频域上表征了抑郁症和正常人的特点[15],但是部分患者存在特异性,其频域的特征表现与正常人一致,这就导致在分类过程中产生了识别率的差异。而DFA 计算的是时间序列的相关性,相对的分类准确率能得到提升。
从算法优缺点看,SASI 提取的特征更容易从表面区分出抑郁症患者和正常人,但是准确程度仍有待提高,DFA 提取的特征更能体现抑郁患者时域特征的变化,患者和正常人的波动趋势更为明显,减少了分类过程的容错率。从算法执行时间看,SASI 的运算效率最快。在图5中,64导联进行频域特征计算所需要的时间在0.3 s 左右,但是DFA 的运算效率不高,从图4 可以看出,随着所选取的数据段的增加,所消耗的时间呈近似于指数量级的增加,消耗的时间甚至需要16 min以上。
从单导联和多导联的区域来看,多导联选出的导联在前额叶、额叶、中心区域、顶区、颞区,如图6 所示,而单导联只反映了一个顶区的特性,多导联能较为宏观地包括整个脑区的生理特性[21]。因此多导联提取的特征质量更高,分类效果更好。
综上,多导联分类的结果相较于单导联更令人满意,对于区分抑郁症患者更有帮助。SASI 算法更适合用于抑郁症的实时反馈阶段以及人机交互部分,而DFA算法则需要对选取的时间窗口进行界定,从单导联的角度进行BCI 技术在抑郁症上的应用,这与DFA 相较于SASI更高的分类准确率一致。
本研究发现多导联的分类结果提升显著,多导联的结果略优于单导联,但实际分类过程中,由于电极帽佩戴及偏移等因素的影响,单导联检测存在一定的偏差,多导联能从数据整体性差异的角度给出一个大脑全局的属性特征,具有一定的稳定性。未来的研究将继续对特征进行优化,改善算法运行效率,提升分类精度,为后续BCI和抑郁症的研究做出贡献。