CN101259016A - 实时自动检测癫痫特征波的方法 - Google Patents

Abstract

本发明公开了一种利用经验模态分解、希尔伯特变换和平滑的非线性能量算子方法来实时自动检测癫痫特征波的方法。实时自动检测癫痫特征波要求先采用移动窗口技术对脑电数据进行分段，然后对每段脑电数据采用如下方法计算：首先采用经验模态分解将原始脑电信号分解为一系列的固有模态函数IMF；然后选取信噪比高的固有模态函数IMF，对其采用自动小波去噪方法进行预处理；接着采用希尔伯特变换计算每个固有模态函数的暂态能量并求和得到暂态总能量S；最后对暂态总能量S应用平滑的非线性能量算子来检测棘波。本发明可用于脑电癫痫特征波的实时自动提取，对于患者的诊断以及减轻医生的繁重劳动都具有重要的意义。

Description

实时自动检测癫痫特征波的方法

技术领域

本发明涉及一种实时自动检测癫痫发作的方法，特别涉及一种大脑癫痫棘波的检测方法。

背景技术

脑电图是诊断癫痫最为重要的辅助检查手段，可以帮助医生确定病人的发作是不是癫痫，尤其对那些不典型发作或相似发作的疾病，可以帮助鉴别诊断。临床上许多疾病，如呼吸暂停、运动障碍、晕厥、心律失常、睡眠障碍、偏头痛及各种神经症状，有时极象癫痫发作，凭借脑电图检查，我们往往可以作出准确的诊断。临床上有典型癫痫发作的患者，大部分可以在脑电图检查中发现癫痫样特征波。癫痫样特征波主要有以下几种：棘波(或尖波)、棘(或尖)慢复合波、以及多棘波群等。因此，癫痫脑电的检测，主要是脑电信号中棘波的检测。从60年代起，癫痫脑电的自动检测技术就受到了注意，许多人在这一领域做了大量的工作，提出了多种检测方法。当前棘波检测算法主要分为五种：第一种是利用经典的信号处理技术，通过非线性滤波器对脑电信号进行预测，将预测结果与输入信号作比较，若两者之差较大，即为癫痫特征波。其优点是对高斯噪声有较好的抑制，缺点是当信号中存在伪差(如眼动，心电信号，肌电信号等)时，误检率比较严重，如果脑电信号特性随时间改变而滤波器参数保持不变，则滤波器的性能会受到影响。第二种是脑电信号与由癫痫波形构成的模板集合计算互相关函数，超过阅值即为癫痫特征波。其缺点是：由于癫痫波形的复杂性，定义一个适合普遍病例的模板集合比较困难，而且，为了在误检率和漏检率之间找到折衷的方法，对阈值的设置也是很困难的。第三种是将脑电信号分解成“半波”(两个相邻幅度极值之间的部分)，提取相对于背景活动的特征参数，将这些特征参数与所设定的阈值比较，超过阈值的视为棘波。其特点是直观，较之前两种方法效果好，但棘波定义的模糊性及伪差波形的不规则性，给棘波的特征提取和判别带来较大的困难。第四种是利用小波变换直接计算脑电信号的时频特征获得脑电图中的锐变点，即棘波。将小波变换技术应用于棘波检测的主要问题是母小波的选择，一般选取形状近似于棘波的母小波效果较好。第五种是通过神经网络学习训练来进行检测。比如CN99124210.6专利就是应用了神经网络，该方法需要训练数据来获取知识，将″学到″的只是蕴含在各神经元与权值中，好像人脑的学习过程。因此需要大量的样本，算法复杂和工作量大，同时不能很好的适应不同的对象。请参考《基于脑电信号分析的癫痫特征检测方法及研究进展》(生物医学工程学杂志，22(3)：606-609，2005)。

脑电图中癫痫波形的实时自动检测是临床上很有意义的工作。然而，自动检测这种癫痫样瞬态特征却非常困难，自动检测的困难在于：第一，所谓的特征波形并没有一个明确的标准，癫痫样特征波形与正常脑电之间的差别基本上是相对的、模糊的；第二，脑电信号是一个非线性时变信号。事实上，无论是从时域还是频域上都很难严格定义一个棘波。因此，不论是基于脑电波形分析的时域自动检测方法还是利用参数模型的频域自动检测方法，实际应用都存在不可避免的困难。

发明内容

本发明针对现有技术的不足，提供了一种利用经验模态分解、希尔伯特变换和平滑的非线性能量算子等来实时自动检测癫痫发作的方法。

本发明自动检测病人脑电图癫痫特征波的方法要求实时性，为此首先采用移动窗口技术对信号进行分段，接着用如下方法进行检测：(1)采用经验模态分解将原始脑电信号分解为一系列的固有模态函数(IMF)；(2)计算固有模态函数(IMF)的信噪比，挑选出高信噪比的固有模态函数(IMF)，对其采用自动小波去噪方法进行预处理；(3)采用希尔伯特变换计算每个固有模态函数的暂态能量并求和得到暂态总能量S；(4)对暂态总能量S应用平滑的非线性能量算子来检测棘波。结果表明这种方法具有较高的精确性，能够有效的检测癫痫棘波。

本发明在临床上主要应用于两个方面：一方面是对长期脑电信号的实时监视，其作用是减少待分析的数据量和脑电信号分析者的工作量，目的是检测出可能是棘波的脑电信号片段，尽可能减少误检率，然后由脑电信号专家对提取出的脑电数据段进行复查，作进一步证实；另一方面是对日常临床数据的检测，其目的是减少视觉检查的过程，即在检测出棘波的同时，保证误检率为零。

所述的经验模态分解采取经验筛法将脑电信号x(t)分解成固有模态函数，表示成

$\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i\left(t\right)+r_N\left(t\right)$

其中，N一般为有限整数，r_N(t)称为残余函数，它代表整个信号的平均趋势。

所述的信噪比计算采取以下的方法：

滤除第一个IMF和最后的残余函数r_N(t)，再计算其余每个IMF傅立叶变换后的幅值谱和信噪比：

A_i＝|fft(IMF_i)|；

${\mathrm{SNR}}_i=\max \left(A_i\right)./\underset{i}{\mathrm{\Σ}}{A}_i$

所述的自动小波去噪预处理采取以下的方法：

对选中的高信噪比的IMF进行多级小波分解，计算多级小波分解系数；接着对多级细节部分的小波系数选取合适的阈值，高于阈值的保留，最后使用初始的近似部分小波系数和经过处理后的细节部分小波系数重构信号。

所述的希尔伯特变换及暂态总能量的计算采取以下的方法：

对于给定函数x(t)，其希尔伯特(Hilbert)变换可以定义为函数x(t)与1/πt的卷积

$H\left[x\left(t\right)\right]=\frac{1}{\mathrm{\π}}p.v.{\∫}_{-\∞}^{\∞}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}$

式中，p.v.——积分的主值。

那么这样就得到预处理后的各固有模态函数IMF_i(t)的解析信号，可表示为

$z_i\left(t\right)={\mathrm{IMF}}_i\left(t\right)+\mathrm{jH}\left[{\mathrm{IMF}}_i\left(t\right)\right]=a_i\left(t\right)e^{{\mathrm{j\θ}}_i\left(t\right)}$

其中

$a_i\left(t\right)=\sqrt{{\mathrm{IMF}}_i^2\left(t\right)+H^2\left[{\mathrm{IMF}}_i\left(t\right)\right]},$ ${\mathrm{\θ}}_i\left(t\right)=\arctan \frac{H\left[{\mathrm{IMF}}_i\left(t\right)\right]}{{\mathrm{IMF}}_i\left(t\right)}$

显然，a_i(t)和θ_i(t)分别表示解析信号z_i(t)的瞬时包络和相位。每个固有模态函数IMF_i(t)的暂态能量为|a_i(t)|²，那么总的固有模态函数的暂态能量为

$S\left(t\right)=\underset{i}{\mathrm{\Σ}}{\left|a_i\left(t\right)\right|}^2.$

所述的平滑的非线性能量算子采取以下的方法：

对脑电信号S(t)＝S₁(t)+S₂(t)，S₁(t)和S₂(t)分别表示背景脑电信号和棘波，且二者互不相关，那么脑电信号S(t)的非线性能量算子为：

Ψ[S(t)]＝Ψ[S₁(t)]+Ψ[S₂(t)]

        ＝[S₁ ²(t)-S₁(t+1)S₁(t-1)]+[S₂ ²(t)-S₂(t+1)S₂(t-1)]

E{Ψ[S(t)]}＝K_S(t)R_S(t，t)＝K_S1(t)R_S1(t，t)+K_S2(t)R_S2(t，t)

其中E{·}表示数学期望，R_S(t，t)表示自相关函数，K_S1和K_S2分别表示背景脑电信号和棘波信号高频带能量与整个频带能量的比值，其中能量的计算采用各信号的Wigner分布得到。

当棘波存在时，上式的第二项较大，因为棘波比背景脑电信号含有更多的暂态能量R_S2(t，t)和暂态高频信息K_S2(t)。然而众所周知，时间序列的Wigner分布会产生交叉项和负频率，因此采用平滑的非线性能量算子，即用一个时域窗和非线性能量算子做卷积运算：

${\mathrm{\ψ}}_{\mathrm{SNEO}}\left[S\left(t\right)\right]=\mathrm{\ψ}\left[S\left(t\right)\right]\⊗w_t$

对改进的平滑的非线性能量算子运用这个非负的能量系数得到E{Ψ_SNEO[S(t)]}，就能解决上述问题，设置一个合理的阈值就可以去检测棘波是否存在。

本发明与现有技术相比有益的效果是：本发明采用的时频分析方法是经验模态分解——希尔伯特变换，该方法能够准确的获得脑电信号在时域和频域的信息，为能量的提取作准备；由本发明获得的脑电信号的非线性能量系数，实现了对癫痫发作特征棘波的实时自动检测，能够实时反映所测脑电信号的非线性特性变化，从而自动地检测出脑电癫痫棘波，克服了人工观察复杂而耗时且需要使用额外的经验判据的缺陷，并且极大地提高了癫痫发作检测的准确率，对于患者的诊断以及减轻医生的繁重劳动都具有重要的意义。

附图说明

下面结合附图和具体实施例对本发明作进一步的详细说明。

图1为本发明的工作流程示意图：

图2A为原始脑电信号和经验模态分解结果；

图2B是图2A中挑选的的高信噪比的固有模态函数经过自动小波去噪预处理的结果；

图2C是图2B中去噪后的脑电信号的暂态总能量S；

图2D是图2C中暂态总能量S应用平滑的非线性能量算子后的结果；

图3A和图3B是癫痫病人的脑电棘波的检测结果；图3A含棘波的脑电信号检测结果；图3B含多棘波群的脑电信号检测结果。

具体实施例

图1为本发明的工作流程示意图。

在步骤101，对得到的原始脑电信号采用移动窗口技术进行分段以便实时处理，这里取每段约1200点。接着对每段数据都按照图1阴影部分内容进行处理。

在步骤102，对每段脑电信号进行经验模态分解计算。经验模态分解是基于信号局部特征、具有自适应性，特别适用于分析大量频率随时间变化的非线性、非平稳信号。经验模态分解采用经验筛法将任意复合信号分解成固有模态函数(IMF)的和，如以下公式：

$x\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i\left(t\right)+r_N\left(t\right)$

其中，N一般为有限整数，r_N(t)称为残余函数，它代表整个信号的平均趋势。

经验筛法计算过程如下：

对于采集到的脑电信号x(t)，首先确定出x(t)上的所有极大值点和极小值点，然后将所有极大值点和所有极小值点分别用一条光滑的曲线联接起来，使两条曲线间包含所有的信号。将这两条曲线分别作为x(t)的上下包络线，计算出它们的平均值曲线m(t)，用x(t)减去m(t)得

x(t)-m(t)＝c₁(t)

从理论上讲，c₁(t)就为x(t)的第一阶IMF，然而实际上，在计算过程中可能产生新的极值、偏移或者增大现有的极值，c₁(t)不满足固有模态函数的条件(固有模态函数满足2个条件：(1)整个信号中，零点数与极点数相等或至多相差1；(2)信号上任意一点，由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零，即信号关于时间轴局部对称；)，因此要求c₁(t)的上下包络线的均值曲线m₁₁，在这第二次筛选过程中，把c₁(t)当作原始数据，得到

c₁-m₁₁＝c₁₁

判断c₁₁是否满足固有模态函数的条件，如果仍然不满足，需要重复多次(k次)筛选过程，

c_1(k-1)-m_1k＝c_1k

直到满足固有模态函数的条件时，得到第一个固有模态函数IMF₁

IMF₁(t)＝c_1k(t)

从x(t)中减去IMF₁得剩余信号r₁(t)，即

x(t)-IMF₁(t)＝r₁(t)

将r₁(t)视为原信号，重复以上过程，就可以获得x(t)的第二阶IMF₂(t)和剩余信号r₂(t)。

通过这样一次次的筛选，就可以得到所有的固有模态函数，从而可以将x(t)表示成

$x\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i\left(t\right)+r_N\left(t\right)$

其中，N一般为有限整数，r_N(t)称为残余函数，它代表整个信号的平均趋势。

如图2A所示，第一个是癫痫病人的原始脑电信号(含棘波)的示意图；imf1是经验模态分解之后的第一阶固有模态函数(IMF)；imf2是经验模态分解之后的第二阶固有模态函数(IMF)；依此类推；共有九阶固有模态函数(IMF)；res代表残余函数r_N(t)。

如图1，在步骤102经验模态分解之后，进行步骤103，计算固有模态函数(IMF)的信噪比，挑选出高信噪比的固有模态函数(IMF)，对其采用自动小波去噪方法进行预处理。

对于分解得到的IMF，其中第一个IMF包含高频噪声，而最后的残余函数r_N(t)表示信号的平均能量，不能够反应棘波能量，因此均被滤除。然后，计算其余每个IMF傅立叶变换后的幅值谱：

A_i＝|fft(IMF_i)|；

再根据幅值谱计算每个IMF的信噪比：

${\mathrm{SNR}}_i=\max \left(A_i\right)./\underset{i}{\mathrm{\Σ}}{A}_i.$

如果分解后的IMF函数超过3个，那么去掉信噪比最小的3个，其余的都选中；如果分解后的IMF函数小于3个，那么全部选中。对选中的这些IMF进行自动一维小波去噪预处理，具体过程是首先对IMF运用symlet8小波进行3级分解，计算3级小波分解系数；接着对从1到3级细节部分的小波系数选取合适的阈值 $\sqrt{2\·\log \left(\mathrm{length}\left({\mathrm{IMF}}_i\right)\right)},$ 并且阈值根据每一层小波分解的噪声水平进行调整。高于阈值的保留。最后使用初始的近似部分小波系数和经过处理后的细节部分小波系数重构信号。

图2B是图2A中挑选的IMF经过自动小波去噪预处理的结果，即Pimf2是第二阶IMF去噪后的结果，依此类推，Pimf7是第七阶IMF去噪后的结果。

如图1，在步骤103预处理信号之后，进行步骤104，希尔伯特(Hilbert)变换的计算和暂态能量的计算。

对于给定函数x(t)，其Hilbert变换可以定义为函数x(t)与1/πt的卷积

$H\left[x\left(t\right)\right]=\frac{1}{\mathrm{\π}}p.v.{\∫}_{-\∞}^{\∞}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}$

式中，p.v.——积分的主值。

那么这样就得到预处理后的各固有模态函数IMF_i(t)的解析信号，可表示为

$z_i\left(t\right)={\mathrm{IMF}}_i\left(t\right)+\mathrm{jH}\left[{\mathrm{IMF}}_i\left(t\right)\right]=a_i\left(t\right)e^{{\mathrm{j\θ}}_i\left(t\right)}$

其中

$a_i\left(t\right)=\sqrt{{\mathrm{IMF}}_i^2\left(t\right)+H^2\left[{\mathrm{IMF}}_i\left(t\right)\right]},$ ${\mathrm{\θ}}_i\left(t\right)=\arctan \frac{H\left[{\mathrm{IMF}}_i\left(t\right)\right]}{{\mathrm{IMF}}_i\left(t\right)}$

显然，a_i(t)和θ_i(t)分别表示解析信号z_i(t)的瞬时包络和相位。希尔伯特变换定义为函数IMF_i(t)与1/πt的卷积，因而经过希尔伯特变换得到的解析信号z_i(t)强调了原始信号IMF_i(t)的局部特性；而其极坐标表达式进一步地表明了其局部特性，即它表示幅值与相位随时间变化的三角函数对IMF_i(t)的最佳局部拟合。

每个固有模态函数IMF_i(t)的暂态能量为|a_i(t)|²，那么总的固有模态函数的暂态能量为

$S\left(t\right)=\underset{i}{\mathrm{\Σ}}{\left|a_i\left(t\right)\right|}^2.$

然后进入步骤105，对暂态总能量S应用平滑的非线性能量算子来检测棘波。本发明提出一个平滑的非线性能量算子用于估计脑电信号的局部能量。平滑的非线性能量算子的计算如下：

对于一个离散时间序列S(t)，非线性能量算子定义如下

Ψ[S(t)]＝S²(t)-S(t+1)S(t-1)

E{Ψ[S(t)]}＝K_S(t)R_S(t，t)

其中E{·}是时间序列S(t)的数学期望，R_S(t，t)是时间序列S(t)的自相关函数，而

$K_S\left(t\right)={\∫}_{-B}^{B}W(t,\mathrm{\ω})\·(1-\cos 2\mathrm{\ωT})\mathrm{d\ω}/{\∫}_{-B}^{B}W(t,\mathrm{\ω})\mathrm{d\ω}$

表示高频带的信号能量与整个频带信号能量的比值，B表示时间序列S(t)的带宽，W(t，ω)是时间序列S(t)的Wigner分布。

令S(t)＝S₁(t)+S₂(t)，S₁(t)和S₂(t)分别表示背景脑电信号和棘波，那么可以得到

Ψ[S(t)]＝Ψ[S₁(t)]+Ψ[S₂(t)]+[{2S₁(t)S₂(t)-S₁(t-1)S₂(t+1)-S₁(t+1)S₂(t-1)]

由于S₁(t)和S₂(t)相互独立，互不相关，因此上式第三项为0。

E{Ψ[S(t)]}＝K_S(t)R_S(t，t)＝K_S1(t)R_S1(t，t)+K_S2(t)R_S2(t，t)

当棘波存在时，上式的第二项较大，因为棘波比背景脑电信号含有更多的暂态能量R_S2(t，t)和暂态高频信息K_S2(t)。然而众所周知，时间序列的Wigner分布会产生交叉项和负频率，因此采用一个时域窗和非线性能量算子做卷积运算，其中选择时域窗为Barlett窗。

${\mathrm{\Ψ}}_{\mathrm{SNEO}}\left[S\left(t\right)\right]=\mathrm{\Ψ}\left[S\left(t\right)\right]\⊗w_t$

对改进的平滑的非线性能量算子运用这个非负的能量系数得到E{Ψ_SNEO[S(t)]}，就能解决上述问题，步骤106设置一个合理的阈值。当超过这一阈值时，为癫痫特征波，当低于这一阈值时，为正常波形。

图2C是图2B中去噪后的脑电信号的暂态总能量S(t)；图2D是图2C中暂态总能量S(t)与平滑的非线性能量算子卷积后的输出结果。图2C的第2个棘波容易被认为是噪声，经处理后得到的图2D更清楚些。

图3A和图3B是癫痫病人的脑电棘波的检测结果。图3A含棘波的脑电信号检测，其中上图是脑电信号，下图是检测结果；图3B含多棘波群的脑电信号检测，其中上图是脑电信号，下图是检测结果。

本发明方法利用经验模态分解、希尔伯特变换和平滑的非线性能量算子等来实时自动检测癫痫发作。该方法应用了经验模态分解方法的自适应分解的特性，将脑电信号分解成了一系列的经过幅值调制和频率调制的固有模态函数IMF，并对这些IMF进行自动小波去噪预处理得到有用的分量及其能量谱，对所有有用能量谱求和，将其与平滑的非线性能量算子做卷积得到输出结果。设定一个合理的阈值，就可以检测出棘波。

上面已结合附图对本发明的具体实施方式进行了示例性的描述，显然本发明并不限于此，在本发明范围内进行的各种形式的改变均没有超出本发明的保护范围。

Claims (8)

1. 一种实时自动检测癫痫特征波的方法，包括如下步骤：

(1)采用经验模态分解将原始脑电信号分解为多个固有模态函数IMF；

(2)计算固有模态函数IMF的信噪比，挑选出至少三个高信噪比的固有模态函数IMF，对其采用自动小波去噪方法进行预处理；

(3)采用希尔伯特变换计算所述每个固有模态函数的暂态能量并求和得到暂态总能量S；

(4)对暂态总能量S应用平滑的非线性能量算子来检测棘波。

2. 根据权利要求1所述的用于实时自动检测癫痫特征波的方法，其特征在于，在所述步骤(1)前，对所述原始脑电信号采用移动窗口技术进行分段处理。

3. 根据权利要求1或2所述的实时自动检测癫痫特征波的方法，其特征在于，所述经验模态分解得到的所述固有模态函数IMF的零点数与极点数相等或至多相差1；且所述固有模态函数IMF上任意一点，由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零，即信号关于时间轴局部对称。

4. 根据权利要求3所述的实时自动检测癫痫特征波的方法，其特征在于，所述步骤(2)中的信噪比计算采取以下的方法：

计算每个IMF傅立叶变换后的幅值谱和信噪比SNR：

A_i＝|fft(IMF_i)|；

${\mathrm{SNR}}_i=\max \left(A_i\right)./\underset{i}{\mathrm{\Σ}}{A}_i.$

5. 根据权利要求4所述的实时自动检测癫痫特征波的方法，其特征在于，所述步骤(2)中采取一维小波去噪预处理的方法是：

首先对IMF进行多级小波分解，并计算多级小波分解系数；接着对多级细节部分的小波系数选取合适的阈值，最后使用初始的近似部分小波系数和经过处理后的细节部分小波系数重构信号。

6. 根据权利要求5所述的实时自动检测癫痫特征波的方法，其特征在于，在所述步骤(3)中，所述预处理后的各固有模态函数IMF_i(t)的解析信号，表示为

$z_i\left(t\right)={\mathrm{IMF}}_i\left(t\right)+\mathrm{jH}\left[{\mathrm{IMF}}_i\left(t\right)\right]=a_i\left(t\right)e^{j{\mathrm{\θ}}_i\left(t\right)}$

其中H[]是希尔伯特变换计算，

$a_i\left(t\right)=\sqrt{{\mathrm{IMF}}_i^2\left(t\right)+H^2\left[{\mathrm{IMF}}_i\left(t\right)\right]},$ ${\mathrm{\θ}}_i\left(t\right)=\arctan \frac{H\left[{\mathrm{IMF}}_i\left(t\right)\right]}{{\mathrm{IMF}}_i\left(t\right)};$

每个固有模态函数IMF_i(t)的暂态能量为|a_i(t)|²，总的固有模态函数的暂态能量为 $S\left(t\right)=\underset{i}{\mathrm{\Σ}}{\left|a_i\left(t\right)\right|}^2.$

7. 根据权利要求6所述的实时自动检测癫痫特征波的方法，其特征在于，在所述步骤(4)中，对得到的所述暂态总能量与平滑的非线性能量算子采取以下的方法：

E{Ψ_SNEO[S(t)]}＝K_S1(t)R_S1(t，t)+K_S2(t)R_S2(t，t)

上式第一项表示背景脑电信号信息，第二项表示棘波信息。

8. 根据权利要求7所述的实时自动检测癫痫特征波的方法，其特征在于，在所述步骤(4)后，选取一个阈值，高于阈值的即认为是棘波。

Images