基于去趋势波动分析(DFA)的脑卒中后抑郁症静息脑电特征提取与识别

王春方 张力新 刘 爽 孙长城 王勇军 赵 欣綦宏志 周 鹏 万柏坤 杜金刚 明 东*

1(天津大学精仪学院生物医学工程系，天津 300072)

2(天津市康复医学研究所，天津市人民医院，天津 300121)

引言

根据世界卫生组织最新报道，脑卒中(stroke)是全球范围内仅次于冠心病和癌症的位于第三位死亡原因的疾病，并且发病率呈逐年上升趋势。我国每年新增脑卒中患者200万以上，每年死于脑卒中的人数超过150万，目前有700万左右的脑卒中幸存者，其中3/4留有不同程度的后遗症，其中重度致残者约占40%，给患者和家庭带来了痛苦和沉重的经济负担［1－2］。

脑卒中后抑郁症(poststroke depression，PSD)作为脑卒中后发病率最高的精神类并发症，其病死率是无抑郁症状的3～4倍，严重威胁着广大脑卒中患者的身心健康［3］。到目前为止，PSD的诊断标准尚未统一，评估方面也无统一的特异性标准可资借鉴，大多借助功能性抑郁症的各种诊断标准作为依据，并借鉴其量表、参数等。国外研究者主要依赖DSM－IV(美国精神障碍诊断与统计手册第4版)诊断标准，国内目前多使用CCMD－Ⅲ(中国精神疾病分类方案与诊断标准第3版)诊断标准。临床抑郁测评量表的选择包括Hamilton抑郁症量表、Beck抑郁自评量表、老年抑郁症量表等，但由于PSD患者常伴有一定程度意识认知障碍或失语等，在配合量表自述或患者家属代述病情时，往往无法准确描述相关情绪改变等信息，使得医师难于准确判断，因而造成诊断准确率较低、漏诊率偏高。归其原因，主要是PSD在诊断方面没有特异性的客观诊断指标。

目前，研究学者普遍认为，去甲肾上腺素(NE)、5－羟色胺(5－HT)等单胺类神经递质水平低下是造成PSD的主要神经生物学原因［4］。脑电信号是皮层细胞群同步活动时突触后电位的总和，作为突触前膜向突触后膜起信息传递作用的化学物质，神经递质的改变必然会引起神经电生理信号(脑电信号)的异常，这在相关文献中也得到了证实［5－9］。

静息状态的脑电信号中包含了大量的信息。已有文献表明，临床诊断抑郁症患者的静息态EEG信号与正常人EEG信号相比，增加的α频段和θ频段功率值与半球间具有不对称性［10－11］。以上发现都是基于时域和频域的传统分析方法。由于EEG信号属于非稳态的时变非线性信号，越来越多的研究学者倾向于用非线性的方法加以分析。已有文献表明，临床诊断抑郁症患者静息状态EEG信号的复杂度、分形维数、瑞利熵等非线性参数要高于正常人［12－13］。但是，PSD 患者的 EEG 信号不仅具有复杂的非线性特性，而且在许多情况下非平稳性和随机性突出，通常的非线性参数(如相关维数、李亚普诺夫指数)与方法皆难以适用(都要求平稳或准平稳条件)。

已有多项研究表明，EEG信号存在长时程相关性的动力学振荡特征［14－16］。去趋势波动分析(detrended fluctuation analysis，DFA)以幂律方式量化时间序列的时域波动，用标度指数描述时间序列的标度行为或长程相关性，适宜研究各种非稳态时间序列的长程幂律形式函数的相关性，其优势在于能够系统地去除数据中由外部因素造成的不同阶次趋势，减少由于不完善的测量方法引起的噪声污染［17－18］。近几年，该方法已被应用于生命科学、地质学，气象学以及经济领域，成为检测非平稳时间序列的有效工具。应用DFA方法分析EEG信号，可为分析其在不同时间尺度上的数据相关性提供定量描述指标，已用于睡眠、麻醉及癫痫发作等状态下EEG信号的特征提取与模式识别。Lee等利用DFA方法，对睡眠呼吸暂停综合征患者的各睡眠分期EEG信号进行了研究，结果发现患者从清醒期到睡眠1期、2期、3期、4期，EEG信号的标度指数值分别为 0.95、1.00、1.12、1.22［19］。随着睡眠程度的加深标度指数值逐渐增大，说明深度睡眠状态下比清醒时的EEG信号具有更为平滑的振荡模式。Gifani等利用DFA方法，研究了不同麻醉深度下的EEG信号，发现从清醒期、浅度麻醉到深度麻醉，标度指数值从0.8上升至2.0［20］。Parish等对癫痫患者的颅内EEG信号进行了去趋势波动分析，发现致癫痫的海马区部位EEG信号标度指数值在0.68～0.88范围内波动，而非癫痫海马区的则在0.60－0.69范围内波动［14］。Lee等对普通抑郁症患者和正常人的EEG信号进行了DFA分析，结果显示抑郁症患者各导联EEG信号标度指数均值为0.90，而与之对应的健康人标度指数均值为0.78，统计学显示存在显著性差异［21］。

本研究尝试利用DFA方法，对正常健康人、脑卒中后无抑郁症患者和脑卒中后抑郁症患者三类被试者静息状态下的EEG信号进行长程幂函数相关性分析，希望通过统计学分析，能提取到相应特征参数，寻找到区分各类被试者的模式识别方法，从而为脑卒中后抑郁症诊断提供具有特异性的客观诊断指标与方法。

1 脑电数据采集

1.1 实验对象

实验对象包括脑卒中患者组11人和健康对照组10人。11名脑卒中患者是天津市南开大学人民医院康复科门诊或住院患者，年龄在43～70周岁之间、(57.8±10.5)岁，其中7名男性、4名女性。11名脑卒中患者中有7名脑卒中后抑郁症患者，4名脑卒中后无抑郁患者。所有脑卒中后抑郁症患者均符合CCMD－Ⅲ抑郁症诊断标准，Hamilton抑郁量表评分均大于15分;所有脑卒中后无抑郁症患者的Hamilton抑郁量表评分均小于15分。10名健康对照组年龄在39～63周岁之间，(56.3±7.2)岁，其中3名男性、7名女性，该组被试的Hamilton抑郁量表评分均小于5分。所有实验对象均为右利手。

1.2 脑电数据采集

EEG数据采集使用NicoletOne 32通道数字视频脑电图仪。电极按照国际标准导联10～20电极系统位置安放(电极阻抗小于10 kΩ)，以Cz作为参考电极，记录16导联头皮EEG信号(Fp1、Fp2、F3、F4、C3、C4、P3、P4、O1、O2、F7、F8、T3、T4、T5、T6);信号采样率为250 Hz，A/D转换精度为12 bit，，滤波范围为0.1～70 Hz。实验条件为远离电磁干扰的弱光环境，要求受试者在保持清醒的状态下闭目休息。在获取EEG数据后，进行变参考、滤波、消除干扰伪迹和独立成分分析(independent component analysis，ICA)等处理，最终选择240 s长度数据进行离线分析。

2 数据分析与分类方法

2.1 DFA分析与标度指数α

去趋势波动分析(DFA)是一种改进的均方根分析方法，它相对于常用的分形分析方法具有以下两个优点:一是能够检测出时间序列信号中表面看似不平稳但其内在深含的自相似性;二是能够避免误检由于外在因素驱动而导致明显的自相似性趋势［22］，具有去伪存真的长处。以每1 s数据为单位时间序列 x(k)，k={1，2，…，N}(N=250)进行DFA分析，具体算法如下:

1)构造时间序列x(k)的去均值的和序列y(i)

式中，〈x〉为时间序列项x(k)的均值。

2)数据长度分割:将和序列y(i)分割为长度为s的不重叠片段，片段个数Ns=N/s，本研究中的分割时间长度取11个值:0.1，0.15，0.2，0.25，…，0.55，0.6 s(对应数据长度分别为:25，38，50，63，75，88，100，113，125，138，150)。

3)去除数据趋势:通过最小二乘法拟合每一个子区间 v(v=1，2，…，Ns)上的局部趋势函数Pv(j)(Pv(j)为第v个片段的二次曲线拟合多项式)，则区间长度为 s的去趋势波动时间序列为Ys(i)，有

4)计算DFA波动函数:分别计算Ns个区间去趋势时间序列的平方均值为

求所有区间均值的平方根，得到DFA波动函数

若原始时间序列x(k)是长程幂函数相关的，则波动函数以幂函数增加，即F(s)∝sa。

5)计算标度指数α

以双对数坐标画出DFA波动函数F(s)和区间长度s的函数关系曲线，并通过线性拟合计算曲线的斜率，即为标度指数α值。

标度指数α为长程幂函数相关性提供了一个定量化指标。若α＜0.5，则表明分段时间序列相互独立;若0.5＜α＜1.0，则表示分段时间序列存在长程幂率形式的持续相关性(具有分形维数的自相似性);若α=1，则表明分段时间序列以1/f噪声［23］方式起伏涨落;若1.0＜α＜1.5，则表明分段时间序列不具备长程相关性;若α=1.5，则表明分段时间序列为布朗噪声，即随机无关。

2.2 独立样本t检验

为进一步探讨3组被试者EEG数据所具有长程相关性的差异，将各被试者16导联EEG信号的标度指数α值作为样本，分别比较了健康对照组与脑卒中患者组的各导联α均值和标准差、脑卒中后无抑郁组与脑卒中后抑郁组各导联α均值和标准差，并对其差异的显著性水平进行了独立样本t检验。

独立样本t检验是一种常用的统计分析方法，用于检验相互独立的两个来自正态总体的样本均值和方差是否来源于同一总体。该方法的基本思想是应用t分布的特征，将t作为统计量来进行检验，其统计量的计算为

2.3 模式识别与分类

为寻找具有特异性的脑卒中后抑郁症客观诊断方法，尝试以被试者静息态下16导联EEG信号的标度指数α值作为特征参数，采用支持向量机(support vector machines，SVM)分类器对上述3组被试者EEG数据进行模式识别研究。SVM是建立在统计学习VC维理论和结构风险最小化原理基础上的一种机器学习方法，在解决小样本、非线性和高维度模式识别问题中表现出突出优势［2］，目前在生物电信号识别、医学图像分割等领域都获得了广泛应用。在分类识别中，使用每名被试者所对应的40个样本值，根据其标度指数α值转换为分类标签值，以判别其所属分组类型。选择多项式作为SVM核函数，第n个样本相应的分类函数为

式中，多项式核函数

每名被试者的分类判决值F(x)为

式中，g为每名被试者的40个样本值中标签值为1的个数，h为可调整阈值参数。

对于健康对照组与脑卒中患者组EEG信号间的分类，1代表健康对照组，－1代表脑卒中患者组，结果验证使用10折交叉验证的方法进行。对于脑卒中后无抑郁症患者组和脑卒中后抑郁症患者组EEG信号间的分类，1代表脑卒中后无抑郁症患者组，－1代表脑卒中后抑郁症患者组，结果验证使用留一法进行。在对健康对照组与脑卒中患者组进行分类识别时，健康对照组的Hamilton量表评分(2.6±2.5)与脑卒中患者组的Hamilton量表评分(15.7±6.8)相差较大，只有当40个样本中绝大部分样本判别为健康对照组时，才将此被试判别为健康人，所以h值选取相对较大;在对脑卒中后无抑郁症组与脑卒中后抑郁症组进行分类识别时，脑卒中后无抑郁症组的Hamilton量表评分(9.3±1.9)与脑卒中后抑郁症组的Hamilton量表评分(19.4±5.6)相差较小，当40个样本中有一定部分样本判别为脑卒中后无抑郁症组时，即可将此被试判别为脑卒中后无抑郁症组，所以h值选取相对较小。

3 结果

3.1 独立样本t检验

图1为健康对照组、脑卒中后无抑郁症组及脑卒中后抑郁症组中3个典型被试C3导联处EEG的DFA函数线性拟合结果，标度指数 α值分别为0.602、0.760、0.921，显示出 3 组被试者 EEG 均具有长程幂律形式的持续相关性。

图1 DFA波动函数的线性拟合图Fig.1 Linear fit figure of DFA fluctuation function

图2为健康对照组与脑卒中患者组各导联标度指数α均值与标准差的比较，图3为脑卒中后无抑郁组与脑卒中后抑郁组各导联标度指数α均值与标准差的比较。从图2和图3中可以看出，上述所有受试者EEG的标度指数α值均在0.5～1之间，说明健康人与脑卒中患者EEG信号均满足长程幂函数相关性(具有分形自相似性)。健康对照组各导联标度指数α均值均小于脑卒中患者组，经独立样本 T 检验，C3、C4、P3、P4、O1、T5、T6 等 7 个导联具有显著性差异(*:P＜0.05);脑卒中后无抑郁组标度指数α均值除P4、F8两导联以外，其他导联均小于等于脑卒中后抑郁组，但经独立样本t检验，各导联值均无显著性差异。

3.2 模式识别与分类

图2 健康对照组与脑卒中患者组的脑电信号标度指数α(*:与健康对照组对比，P＜0.05;＊＊:与健康对照组对比，P＜0.01)均值与标准差比较Fig.2 The mean and the standard deviation of the scaling exponent for healthy controls and poststroke subjects(*:compared with healthy controls，P ＜0.05;＊＊:compared with healthy controls，P ＜0.01)

图3 脑卒中后无抑郁组与脑卒中后抑郁组的脑电信号标度指数α均值与标准差比较Fig.3 The mean and the standard deviation of the scaling exponent for poststroke depression subjects and poststroke non-depression subjects

对于健康对照组与脑卒中患者组，EEG信号间的10折交叉验证结果显示，当阈值参数h取25时，分类正确率可以达到85.0%;当h取30时，分类正确率可以达到100%。这表明，用标度指数α值作为特征参数来区分健康人和脑卒中患者，具有很好的识别效果。对于脑卒中后无抑郁症患者组与脑卒中后抑郁症患者组EEG信号的留一法验证结果显示，当h取10时，分类正确率可以达到72.7%;当h取9时，分类正确率可以达到90.9%。这表明，用标度指数α值作为特征参数区分脑卒中后无抑郁症患者和脑卒中后抑郁症患者也有较好的识别效果。

4 讨论和结论

表征脑卒中患者和健康人静息态EEG信号去趋势波动分析(DFA)的标度指数值均在0.5～1之间，表明两类人群静息态EEG信号均具有长程的幂律关联性，符合长程幂率形式的分形系统动力学行为特征，EEG信号波动行为在不同时间尺度上具有某种相似性。健康人各导联静息态EEG信号去趋势波动分析(DFA)的标度指数值均小于脑卒中患者，在 C3、C4、P3、P4、O1、T5、T6 等 7 个导联具有显著性统计差异。由分形动力学可知，当标度指数值在0.5～1之间时，值越大，说明信号的自相似规律性较强，意味着信息熵较低(值=1时，表现为1/f噪声);反之值越小，则信号的随机性较强，意味着信息熵较高(值=0.5时，表现为白噪声)。这说明，脑卒中患者静息态EEG信号相对于健康人具有更强的自相似规律性，振荡模式更为平滑，可能是脑卒中患者脑部结构病理性改变所致。关于造成不同人群及不同生理病理状态下EEG信号表现出的长程幂函数相关性的差异性原因，还没有完全清楚，一些研究推测，这与神经振荡产生的特殊机制相关，标度行为或长程相关性的差异反映着神经元振荡可能受到与其起源相关的特殊机制的影响［18］;本研究推测，脑卒中引起的神经元坏死等脑部病理性改变使得EEG信号振荡源特性发生了改变，进而影响了EEG信号的振荡模式。有关研究表明，标度指数值可以有效地量化大脑自发振荡的潜在动力学过程特征［16］，不同节律EEG振荡有其特异性标度指数值。基于此，本研究下一步的工作重点之一为:利用去趋势波动分析(DFA)，研究不同节律下各组人群EEG信号的长程幂函数相关性特征。脑卒中后无抑郁症患者绝大多数导联静息态EEG信号去趋势波动分析(DFA)的标度指数值小于脑卒中后抑郁症患者，表明抑郁状态对脑电生理活动的长程幂相关性也有所影响。抑郁症患者静息态EEG信号振荡模式更为平滑，衰减更为缓慢。本研究推测，脑卒中后抑郁症患者由于特定部位脑组织(尤其是神经细胞)受损引起神经递质释放减少，相对于脑卒中后无抑郁症患者神经元振荡模式有所不同。虽然本研究EEG实验数据未检验出其间显著性统计学差异，但利用支持向量机(SVM)分类工具，将DFA的特征参数值放入高维特征空间进行模式识别，仍可以取得高达90%的分类正确率，正说明DFA具有前述能够检测出时间序列信号中表面看似不平稳但其内在深含的自相似性，即去伪存真的长处。未检验出统计学意义上的显著性差异可能是由于样本量的限制，本研究将继续收集脑卒中后两类人群的静息态EEG信号做进一步分析，以得到更加显著可靠的结论。

通过对11名脑卒中患者(其中7名患抑郁症，4名未患抑郁症)和10名健康对照组受试者静息态EEG信号进行去趋势波动分析(DFA)，提取其特征参数(标度指数值)，利用支持向量机(SVM)分类器进行模式识别，可以得到以下初步结论:一是脑卒中患者和健康人静息态EEG信号均具有长程的幂律关联性;二是相对于健康人，脑卒中患者的静息态EEG信号具有更强的自相似规律性，振荡模式更为平滑，而脑卒中后抑郁症患者较无抑郁症患者此现象更为严重。

综上所述，笔者认为可以考虑基于受试者静息状态下EEG信号的去趋势波动分析(DFA)，提取其标度指数值作为脑卒中后抑郁症的特征参数，再利用支持向量机(SVM)进行高维特征空间模式识别，这有望为临床脑卒中后抑郁症(PSD)提供客观有效的辅助诊断手段，所以值得深入继续研究。

［1］王文志.中国脑卒中流行病学特征和社区人群干预［J］.中国医学前沿杂志，2009，1(2):49 －53.

［2］丁世飞，齐丙娟，谭红艳.支持向量机理论与算法研究综述［J］. 电子科技大学学报，2011，40(1):2 －10.

［3］Angelelli P，Paolucci S，Bivona U，et al.Development of neuropsychiatric symptomsin poststroke patients:a crosssectional study［J］.Acta Psychiatr Scand，2004，110(1):55 －63.

［4］何舒，林渝峰，李霞.抑郁症的发病机制研究进展［J］.四川生理科学杂志，2006，28(3):126 －128.

［5］Baskaran A，Milev R，McIntyre RS.The neurobiology of the EEG biomarker as a predictor of treatment response in depression［J］.Neuropharmacology，2012，63(4):507－513.

［6］Koek RJ，Yerevanian BI，Tachiki KH，et al.Hemispheric asymmetry in depression and mania A longitudinal QEEG study in bipolar disorder［J］.Journal of Affective Disorders，1999，53(2):109－122.

［7］陈爱美.基于神经振荡的抑郁症脑电信号分析［D］.秦皇岛:燕山大学，2010.

［8］胡巧莉.基于相位同步分析方法的抑郁症脑电信号的研究［J］. 中国医疗器械杂志，2010，34(4):246 －249.

［9］Davidson RJ.Anterior electrophysiological asymmetries，emotion，and depression:Conceptual and methodological conundrums［J］.Psychophysiology，1998，35(5):607 －614.

［10］Pollock V，Schneider L.Quantitative，waking EEG research on depression［J］.Biol Psychiatry，1990，27(7):757 －780.

［11］Henriques JB，Dividson RJ. Regional brain electrical asymmetriesdiscriminatebetween previouslydepressed and healthy control subjects［J］.J Abnorm Psychol 1990，99(1):22－31.

［12］Ahmadlou M，Adeli H，Adeli A.Fractality analysis of frontal brain in major depressive disorder ［J］.Int J Psychophysiol，2012，85(2):206－211.

［13］李颖洁，邱意弘，朱贻盛.脑电信号分析方法及其应用［M］.北京:科学出版社，2009:111－117.

［14］梁晓花.基于脑电心电数据融合的睡眠分期［D］.南京:江苏大学，2008.

［15］Watters PA，Martin F.A method for establishing long－range power law correlations from the electroencephalogram［J］.Biol Psychiatry，2004，66:79 －89.

［16］Nikulin VV，Brismar T.Long－range temporal correlations in electroencephalographic oscillations:relation to topography，frequency band，age and gender［J］.Neuroscience，2005，130:549－58.

［17］Watters PA，Martin F.A method for estimating long－range power law correlations from theelectroencephalogram ［J］. Biol Psychiatry，2004，66(1):79 －89.

［18］Nikulin VV，Brismar T.Long－range temporal correlations in electroencephalographic oscillations:relation to topography，frequency band，age and gender［J］.Neuroscience，2005，130(2):549－558.

［19］Lee JM，Kim DJ，Kim IY，et al.Nonlinear－analysis of human sleep EEG using detrended fluctuation analysis［J］.Med Eng Phys，2004，26(9):773 －776.

［20］Gifani P，Rabiee HR，Hashemi MH，et al.Optimal fractalscaling analysis of human EEG dynamic for depth of anesthesia quantification［J］.J Franklin Inst，2007，344:212 －229.

［21］Lee JS，Yang BH，Lee JH，et al.Detrended fluctuation analysis of resting EEG in depressed outpatients and healthy controls［J］.Clinical Neurophysiology，2007，118:2489 －2496.

［22］Parish LM，Worrell GA，Cranstoun SD，et al.Long－range temporal correlations in epileptogenic and non－epileptogenic human hippocampus［J］.Neuroscience，2004，125(4):1069 －1076.

［23］Bak P，Tang C，Wiesenfeld K.Self－organized criticality:an explanation of the 1/f noise［J］.Phys Rev Lett，1987，59(4):381－384.