基于脑电信号HilbertHuang变换的睡眠分期研究

【摘要】  目的 研究基于脑电信号HilbertHuang变换的睡眠自动分期方法。方法 对睡眠脑电信号进行HilbertHuang变换，求出具有物理意义的瞬时频率，并得到脑电信号在频率上的能量分布，作为睡眠脑电信号各个时期的特征，终利用近邻模式分类方法对睡眠各阶段进行分期决策。结果 通过对560个睡眠脑电信号样本进行分期，平均正确率达到81.7%。结论 经HilbertHuang变换得到的睡眠脑电信号特征，可以作为睡眠分期的有效分类依据。

【关键词】  HilbertHuang变换 脑电信号 睡眠分期 固有模态函数

    Automatic Sleep Stage Classification ba[x]sed on HilbertHuang Transform Method of EEG.LI Gu,FAN Yingle,LI Yi,PANG Quan.Space Medicine & Medical Engineering,2007,20(6):458~463

    Abstract: ob[x]jective To present a new application of HilbertHuang transform method for automatic sleep stage classification. Methods Instantaneous frequencies of sleep EEG with physical meaning and energyfrequency distribution used as feature parameters for each stage were computed with the HilbertHuang transform. Finally, pattern recognition method of nearest neighbors was applied to optimal classification. Results Five hundred and sixty samples were picked from sleep EEG , and classification was made.The mean rate of accuracy was as high as 81.7%. Conclusion Automatic sleep stage classification can be made effectively from features of sleep EEG by HilbertHuang transform.

    Key words: HilbertHuang transform;EEG; sleep stage classification; intrinsic mode function

    Address reprint requetsts to: LI Gu. Biomedical Engineering & Instrument Institute, Hangzhou Dianzi University, Hangzhou Zhejiang 310018, China

    睡眠作为一种复杂的生理过程，是机体复原、整合和巩固的重要环节。目前国际上普遍使用R&K睡眠EEG(electroencephalogram)分期规则，根据睡眠时脑电信号的表现，将睡眠分为：觉醒期、非快速眼动睡眠期（nonrapid eye movement, NREM）和快速眼动睡眠期（rapid eye movement, REM）。其中NREM又可分为Ⅰ、Ⅱ、Ⅲ、Ⅳ期。睡眠分期研究，对于睡眠状态分析以及睡眠质量的科学评估，具有重要的应用价值。

    在早期临床睡眠分析中，主要是由专家对整夜的连续睡眠数据进行人工目测分析，因此评估质量很大程度上依赖于专家的经验和水平［1］。从20世纪70年代开始，国际上开始了自动睡眠分期的研究，例如Smith应用Walsh谱分析来研究睡眠脑电信号，但是谱分析方法建立在信号是平稳随机过程的假设之上，而睡眠脑电信号并不满足平稳性的假设，因此Walsh谱分析方法在分析脑电信号的应用中具有不可避免的局限性。近期有人用小波来研究睡眠脑电信号［2］，小波是傅立叶分析的发展，是一种多尺度的信号分析方法，具有良好的时频局部化特性，非常适合分析非平稳信号的瞬态特征和时变特性，在脑电信号的分析和研究中取得了很大的进展［3］。但是在小波分析中所使用的小波函数具有多样性，因此分析得到的小波分量和小波谱只相对所选择的小波基有意义。

    HilbertHuang变换是由Huang N等［45］人提出的一种新的针对非平稳非线性信号的时频分析方法，它基于信号的局部特征时间尺度，能把复杂的信号分解为有限的内在模式函数之和，每一个固有模态函数(intrinsic mode function，IMF）所包含的频率成分不仅与采样频率有关，重要的是随信号本身的变化而变化，因此它非常适用于研究非平稳的信号。其已在生物医学信号［6］ 、地震信号及在振动工程领域中的故障检测、参数识别［78］等方面进行了广泛应用。考虑到睡眠脑电信号的频域特征表现比较突出，可以说是由多种频率和振幅的波混合而成，具有相当强的非平稳性和非线性［9］，因此本文采用HilbertHuang变换来研究睡眠脑电信号，提取睡眠脑电信号各个时期的特征，为觉醒期、NREM期和REM期等时期的自动判决提供特征选择的依据。

    方  法

    睡眠脑电的生理特征  从觉醒到深睡眠的过程中，脑电信号的频率变化非常明显，首先α波逐渐消失，频率变慢，然后出现θ波及δ波，并且在睡眠过程中还会出现一些特殊的脑电现象，如出现K复合波和睡眠纺锤波。K复合波是频率在1 Hz左右的低频高幅慢波，而睡眠纺锤波是一种每隔一段时间重复出现的振荡信号，频率范围在11~16 Hz左右，持续时间大于0.5 s。虽然可将脑电信号按频率划分成4种波形，但是实际记录到的脑电信号往往是各种波段同时存在的极不规则的信号。

    HilbertHuang变换原理                    航天医学与医学工程                                    第20卷第6期                      李  谷，等.基于脑电信号HilbertHuang变换的睡眠分期研究HilbertHuang变换引入了基于信号局部特征的固有模态函数，因此适用于分析非线性、非平稳信号。其主要由两部分组成，即经验模式分解（empirical mode decomposition, EMD）和Hilbert变换。

    经验模式分解方法目的是通过对信号进行分解，获得固有模态函数。所获得的固有模态函数必须满足:1) 在整个信号长度上，一个IMF的极值点和过零点数目必须相等或至多只相差一点；2) 在任意时刻，由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零，即IMF的上下包络线对称于时间轴。

    获得信号IMF分量后，就可以对各IMF分量做Hilbert变换，从而求出每个IMF的瞬时频率、幅度函数以及相应的Hilbert谱，汇总所有IMF的Hilbert谱就得到原始信号的Hilbert谱。

    睡眠脑电特征提取  通过HilbertHuang变换提取睡眠脑电信号各期特征的具体步骤如下：首先找到一段脑电信号x(t)的极大值和极小值，通过3次样条拟合，获得信号的上包络曲线u(t)和下包络曲线v(t)，计算上下包络曲线在每一点上的平均值，获得平均值曲线m(t)，即：

    m(t)=［u(t)+v(t)］/2(1)

    用x(t)剪去m(t)后剩余部分h1(t)，即：

    h1(t)=x(t)-m(t)(2)

    从理论上讲，h1(t)即为阶IMF分量，但是由于包络线条样条逼近的过冲和俯冲作用，新的极值点会影响原来极值的位置与大小。分解得到的h1(t)并不完全满足IMF的两个条件，因此用h1(t)代替 ，求出与h1(t)相应的上、下包络线为u1(t)和v1(t)，重复以上移动过程k次，当所得到的hk(t)满足上述IMF的两个条件，即得到阶IMF分量c1(t)。

    c1(t)=hk(t)(3)

    接着从原始信号中减去c1(t)，即可获得信号的残余部分r1(t):

    r1(t)=x(t)-c1(t)(4)

    对r1(t)重复上面的过程，可获得第二阶IMF分量。通过EMD方法对信号进行逐次筛选，可获得信号的多个IMF分量和一个残余函数Rn，从而信号可由式(5)表示：

    x(t)=∑ni=1ci(t)+Rn(5)

    对睡眠脑电信号x(t)进行经验模式分解，求得9个IMF分量的示意图，如图1所示。

    通过EMD方法得到脑电信号各阶IMF后，就可以利用Hilbert变换求每个IMF分量的瞬时频率、幅度函数以及相应的Hilbert谱和边际谱。

    对所求得的每一个IMF分量ci(t)作Hilbert变换：

    H［ci(t)］=1π∑＋∞-∞ci(τ)t-τdτ(6)

    构造解析信号zi(t) ：

    zi(t)=ci(t)+jH［ci(t)］=Ai(t)ejθi(t)(7)

    其中幅值函数Ai(t)=c2i(t)+H2［ci(t)］ ，相位函数θi(t)=arctgH［ci(t)］ci(t)

    进一步得到各IMF分量的瞬时频率：

    wi(t)=dθi(t)dt(8)

    图2是各个IMF分量的瞬时频率。从图2中可以看出，通过EMD分解，阶IMF分量的瞬时频率大，接下来的分量瞬时频率逐渐变小。这是由EMD分解的特点决定的，其先对信号的图1  睡眠脑电信号的IMF分量

    Fig.1  Sleep EEG IMFS

    图2  脑电信号的IMF分量的瞬时频率

    Fig.2  Instantaneous frequencies of sleep EEG IMFS

    高频成分进行提取，然后依次分解出低频和更低频的信号分量。EMD的分解过程可以解释为多尺度滤波过程，每个IMF分量都反应了信号的特征尺度，代表着非线性非平稳信号的内在模态特征。与傅立叶方法的全频率成分分解相比，EMD方法的分解方式由局部频率成分决定，它的基函数直接来源于信号本身，因此分解具有局域性和自适应性。

    省略残余函数Rn，对每一阶IMF作Hilbert变换，并求出相应的解析函数的幅值谱和瞬时频率，从而原始信号可以表示为：

    x(t)=∑ni=1Ai(t)ej∫ωi(t)dt(9)

    将(9)式的右部取实部，定义它为Hilbert幅值谱，简称Hilbert谱，记作:

    H(w,t)=Re∑ni=1Ai(t)ej∫ωi(t)dt(10)

    将Hilbert谱对时间进行积分，可以获得信号的Hilbert边际谱：

    h(w)=∫+∞-∞H(w,t)dt(11)

    边际谱是Hilbert谱对时间的积分，它的含义是：信号中（瞬时）频率的总幅值大小，也可以这么认为，瞬时频率表示信号交换快慢的物理量，任一瞬时频率都有一定的能量，将所有时刻某一瞬时频率的能量加起来就是瞬时信号中该频率的总能量，即边际谱的高度。这一物理意义和傅立叶频谱的物理意义是一致的，只不过在傅立叶频谱中，对任意频率都要求具有相同的幅值，这就容易图3  睡眠脑电各时期边际谱破坏信号中本来的真实频率，不适合运用到脑电这种非平稳的信号中。边际谱比傅立叶频谱更加准确地反映信号的实际频率成分，更加准确地提取脑电信号中的频率特征。通过边际谱，可以看出一定时间内睡眠脑电信号的能量主要集中在什么频率上，进一步可以判断出该时刻所处的睡眠阶段。图3表示各个阶段睡眠脑电的希尔伯特黄变换边际谱幅值。

    通过实验中HHT的边际谱特性，可以获得睡眠脑电信号各个时期的频率组成以及幅值随频率的动态变化情况，得到睡眠各个时期的特征。各睡眠时期脑电信号特征总结如下：

    1） 觉醒期：一般为连续的α波。

    2） Ⅰ期：即思睡——入睡期。是清醒至睡眠之间的过渡阶段。从这个时期边际谱上可以看出，有一部分的α波，但是相对于觉醒期来说α波逐渐减少，所以频率变慢，并且含有θ波。Ⅰ期约占总睡眠时期的5%。

    3） Ⅱ期：浅睡期。这一时期睡眠比Ⅰ期睡眠要深，Ⅱ期的边际谱表明其特点是明显地看见K复合波和纺锤波。频率也在逐渐变慢， δ波增多。这一时期约占总睡眠时期的44%~55%。

    4） Ⅲ期：中度深睡期。Ⅲ期图中的δ波明显增多。

    5） Ⅳ期：深睡期。从Ⅳ期的图中可以看出纺锤波消失， δ波比Ⅲ期的多。Ⅲ期和Ⅳ期这两个时期合起来约占总睡眠时间的15%~23%。

    6） REM期：REM图表明， δ波明显减少，有θ波，有时候还有一些α波。这一时期占总睡眠时间的20%~25%。

    从边际谱中可以看出睡眠脑电信号频率范围主要集中在0~30 Hz。为了能更加有效地对睡眠进行分期，因此在睡眠脑电信号原频带划分的基础上进一步细化。借鉴小波变换中对睡眠脑电信号频带划分的情况［2］，本文将脑电信号分为7类，如表1所示。表1  睡眠脑电信号频带划分 按照上述睡眠阶段的特点，进行特征选择。首先求出如表1所示的各个频带的能量Ei，i=1,2,…7，接着求出频带的总能量E8

    E8=∑7i=1Ei(12)

    睡眠脑电信号较复杂，每个个体各频段脑电信号的能量大小差异较大，但是脑电信号在频率上的能量分布相对比较稳定，因此把各个频段的能量与总能量的比值作为睡眠分期的特征参数。按前述各个时期的特征，本文设置了如表2所示的6个特征参数。表2  睡眠脑电特征参数

　　结  果

    本文所用的睡眠脑电数据来源于MITBIT生理信息库中的EDF(european data format)格式的睡眠数据库。在这个数据库中，共有8位被试者，年龄在21~35岁，健康白种人。所有记录的数据包括心电信号（electrocardiograph, ECG）、脑电信号、肌电信号 (electromyogram，EMG)信号等等，数据的采样频率为100 Hz。每一份记录的数据都附带注释信息，注释是由经验丰富的医生根据R&K规则，以30 s为一个单位所进行的人工分期的结论，以此分期结论作为研究各睡眠时期特征的参考标准。

    一般情况下，人的睡眠时间为6 h，为了更准确地提取特征，本文选择时间从00∶00 am~05∶00 am连续5 h，由C3电极记录下来的脑电数据。根据R&K分期规则，以30 s为一个单位时间，分别对每一个单位时间的脑电数据进行EMD分解，得到IMF分量，接着求出由表2所示的6个特征参数。

    模式分类采用近邻法进行自动睡眠分期。由于Ⅲ期和Ⅳ期波形相似，只是在振幅上略有区别［2］，因此在本文实验中将Ⅲ期和Ⅳ期划分为一类。Ⅰ期是一个过渡阶段，非常敏感，在整个睡眠时期内所占的时间非常短，只占5%左右，从实验数据中发现，根据专家的分期结果，Ⅰ期与Rem期经常交替出现，而且仅通过EEG的频率特征往往无法直接区别出Ⅰ期和Rem期，一般还需要更多的信息如眼动信息等，因此本文实验中将Ⅰ期和Rem期划分为一类。按照人工分期结果，通过近邻法，在被试者 脑电信号中分别选取觉醒期、Ⅱ期、Ⅲ/Ⅳ期和Ⅰ期/Rem期各140个信号，总共560个信号。将本文方法得到的自动分期状态(A)，与专家的人工分期状态(M)进行比较，符合率如表3所示。

    讨  论

    由表3可以看出，通过HilbertHuang变换提取睡眠各个时期的特征，再利用近邻法对睡眠进行分期平均正确率达到了81.7%，与利用功率谱估计［1］和小波变换［2］进行自动分期时的平均正确率79.4%和77.6%相比，利用HilbertHuang变换对睡眠脑电进行分期的分期准确性相对较高。人睡眠时脑电信号在频率上的变化非常明显，各个时期的特征与觉醒期差别很大，因此觉醒期的准确度比较高。Ⅲ期和Ⅳ期属于深睡期，因此相对于Ⅰ期、Ⅱ期和Rem期来说眼动、体动等干扰较小，相对比较平稳，因此正确率也比较高。Ⅱ期中K复合波和纺锤波是周期性突然出现的，但睡眠是一个连续的过程，因此有时虽然处在Ⅱ期中，但是没有K复合波和纺锤波，这就使得Ⅰ期、Ⅱ期和Rem期分期的正确率较低。需要说明的是睡眠分期是为了在研究睡眠过程时更加方便，人为的依据脑电信号和生理表现而规定睡眠所处的阶段，但是睡眠过程本身是连续的，没有明显的分界线，往往是逐渐变化，重叠交错。表3   睡眠分期实验结果

    根据R&K规则，基本上30 s反应一个睡眠阶段，但是由于睡眠各个阶段之间的分界并不明显，有时候30 s内可能并不只存在一个睡眠阶段，只能说在这个30 s有可能处于某个阶段的睡眠状态，因此不同的判断方法往往会有不同的判断结果。另外个体特征存在着差异，可能也是造成错判的原因。由于人体的复杂性，可以增加除睡眠脑电信号以外的信息，如心动周期，呼吸波周期等等来提高分期结果的准确度。

    结  论

    针对睡眠脑电信号的非平稳特性，本文提出基于一种新的非平稳信号处理技术——HilbertHuang变换对睡眠进行自动分期的方法。HilbertHuang变换的特点是引入了EMD分解，经EMD分解得到的IMF分量完全是从脑电信号本身分离出来的，不需要任何先验知识，不受人为的影响，因此IMF分量能够更好反映脑电信号的固有特性，其分解是客观的，自适应的。而通过HilbertHuang变换得到的瞬时频率有较好的时频分辨率以及较为清晰的物理意义。实验结果表明HilbertHuang变换能够有效地提取睡眠各个时期的特征，为睡眠自动分期提供很好的依据，使得分期的准确度比较高。因此HilbertHuang变换将成为生物医学信号处理领域内一个新的有效的分析工具。但是HilbertHuang变换理论本身还存在一些未解决的问题，例如波动问题和端点效应等，都需要进一步研究和解决。

【参考文献】
  ［1］ Mutapcic A,Shimayama T,Flores A. Automatic sleep stage classification using frequency analysis of EEG［C］.Proceedings of XIX International Symposium on Information and Communication Technologies，2003.

［2］ Oropesa E, Cycon HL,Jobert M. Sleep Stage Classification using Wavelet Transform and Neural Network ［R］. ICSI Technical Report , TR99008,1999.

［3］ LI Yong, ZHANG Shengxun. Apply Wavelet Transform to analysis EEG signal［C］.18th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Amsterdam，1996： 10071008.

［4］ Huang NE, ZHEN Shen, Long SR.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［C］. Proceedings of the Royal Society of London, London 1998,454(A):903995.

［5］ Huang NE, Wu ML,Long SR. A confidence limit for the empirical mode decomposition and Hilbert spectral analysis ［C］. Proceedings of the Royal Society of London,London 2003,459(a):23172345.

［6］ Huang W, Shen Z,Huang NE, et al .Nonlinear indicial response of complex nonstationary oscillations as pulmonary hypertension responding to step hypoxia［C］.Proceedings of the National Academy of Sciences USA,1999,96(5):18341839.

［7］ Yang JN, Lei Y, Pan S, et al. System identification of linear structures ba[x]sed on hilberthuang spectral analysis［J］.Earthquake Engineering and Structural Dynamics,2003,32(9):14431467.

［8］ Huang NE. The age of large amplitude coastal seiches on the Caribbean coast of Puerto Rice［J］. Phy Oceanography，2000, 30(8):405409.

［9］ Veltcheva AD. Wave and group transformation by a hilbert spectrum［J］.Coastal Engineering Journal，2002,44(4):283300.