CN110584598B - 基于时频分析的自动分离高频振荡节律中尖峰活动方法 - Google Patents

Abstract

本发明公开了基于时频分析的自动分离高频振荡节律中尖峰活动方法。首先，采集癫痫病人脑电信号，并提取疑似高频振荡节律所对应的原始信号；然后，对疑似高频振荡节律所对应的原始信号进行尖峰信号检测，并将检测到的尖峰信号进行剔除，得到高精度的高频振荡节律；其中，尖峰信号的检测采用了基于时频分析的方法，使用箱型图获得信号功率的上四分位和下四分位值，检测出高频振荡节律中混有的尖峰活动，提高癫痫脑电信号的高频振荡节律的检测精度；一种高频振荡节律检测设备及存储设备，实现了高频振荡节律检测方法。本发明的有益效果：提高癫痫脑电信号的高频振荡节律的检测精度，帮助医生进行癫痫诊断和致痫灶切除。

Description

基于时频分析的自动分离高频振荡节律中尖峰活动方法

技术领域

本发明涉及一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，属于癫痫脑电信号处理领域。

背景技术

癫痫是一种常见的神经系统疾病，具有自发性和不可预测性等特征。癫痫发作时，患者常表现出运动和行为等异常。全世界癫痫患病率为0.5％-1％，其中我国为0.7％，且仍在逐年增长。大部分癫痫患者可以通过抗癫痫药物进行治疗，但仍有约30％的患者不能通过药物控制癫痫发作，从而被确诊为难治性癫痫患者。

难治性癫痫是指经过常规的、系统的抗癫痫药物治疗，其血液内抗癫痫药物浓度保持在有效的范围内，但仍不能控制病人的癫痫发作，并严重影响病人的工作、学习或正常生活者；其发作频率在每月2-4次以上，病程在4年以上。对于难治性癫痫，控制及治愈的有效措施之一是手术切除癫痫病灶，而手术成功的关键是病灶的精确定位。目前，癫痫病灶定位的主要手段有脑电图检查、神经电生理学检查、核医学检查、脑磁图检查等。其中脑电图检查是最常用的癫痫病灶定位手段，它通常需要通过植入颅内或硬膜下大电极来长时间记录脑电信号，捕捉癫痫患者发病过程中的脑电变化，形成脑电图，然后由有经验的医学专家进行即时观测分析。医学专家通过检查癫痫脑电信号，实现癫痫诊断。

致痫灶是患者癫痫发作期或者发作间期的病变脑皮质区。通常致痫灶分为发作始发区，激惹灶，致痫病变和功能缺失区，其中发作始发区是最为有效和显著的标志，可以看做是致痫灶的替代指标。

传统的脑电图方法重点关注的是40Hz以下的脑电信号，通过提取该频段脑电频率成分来定位癫痫发作始发区。但是该方法检测的低频频率成分很容易受到其他信号的干扰，使得定位结果出错，导致切除手术失败。同时，该方法非常耗时，定位时间大约为24-72小时，加大了手术风险。因此，寻找一种新的标志物或者快速定位方法迫在眉睫。

近二十年来，越来越多的研究者们开始关注80Hz以上癫痫脑电信号的高频振荡节律。在癫痫发作潜伏期，患者大脑开始发生病理性变化。在病变前后的大脑中，高频振荡节律有着显著的不同。大量研究表明高频振荡节律可以作为癫痫发作始发区的生物标志，它们在癫痫发作始发区的发生率更高，具有明显的特异性，所以可以通过检测高频振荡节律来实现发作始发区的快速精确定位。通过发作始发区的精确定位，可以明确癫痫病灶切除手术目的，降低术中因损伤功能区而造成不可恢复的脑神经损伤的风险。

发明内容

本发明的目的是提出一种基于时频分析的自动分离高频振荡节律中尖峰活动方法。一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，包括以下步骤：

S101：采集癫痫病人的原始脑电信号，并对所述原始脑电信号进行50Hz倍频陷波和80-500Hz带通滤波，得到去除50Hz倍频并且频率在80-500Hz的预处理后的脑电信号，并提取所述预处理后的脑电信号中的疑似高频振荡节律所对应的原始信号；具体如下：

首先对所述预处理后的脑电信号加长度为L的时间窗，得到处理后的脑电信号；其中，L为预设值；

然后对所述处理后的脑电信号使用均方根算法，得到滤波后的癫痫脑电的均方根信号RMS(t)，并记录其持续时间T及波峰个数m；

最后设定均方根信号阈值Threshold_RMS，并判断条件RMS(t)>Threshold_RMS&&T>30ms&&m≥6是否成立？若是，则认为RMS(t)为疑似高频振荡节律，进行提取，从而得到疑似高频振荡节律所对应的原始信号；反之则认为RMS(t)为非高频振荡节律，予以遗弃；

S102：对所述疑似高频振荡节律所对应的原始信号进行尖峰信号检测，并将检测到的尖峰信号进行剔除，进而得到高精度的高频振荡节律；具体如下：

首先采用时频分析方法对所述疑似高频振荡节律所对应的原始信号进行50Hz倍频陷波和0-500Hz带通滤波，得到去除50Hz倍频并且频率在0-500Hz的事件信号x(τ)；其中，所述疑似高频振荡节律所对应的原始信号为步骤S101中的所述原始脑电信号；

对x(τ)进行S变换得到其时频信息S(t,f)，并对时频信息S(t,f)进行奇异值分解，得到分解后的奇异值；再使用k-mediods算法对所述分解后的奇异值进行处理，得到有效的奇异值；结合有效奇异值，对S(t,f)进行逆S变换，得到重构后的信号x′(τ)；

最后设定功率下界异常阈值Threshold，对x′(τ)进行阈值处理，检测出尖峰信号并剔除，进而得到高精度的高频振荡节律。

进一步地，在步骤S101中，获取疑似高频振荡节律具体过程如下：

所述均方根信号RMS(t)的计算公式如式(1)所示：

式(1)中，s(k)为所述处理后的脑电信号，N_lL为所述处理后的脑电信号的长度；所述均方根阈值Threshold_RMS的设定公式如式(2)所示：

Threshold_RMS＝3*SD+E (2)

式(2)中，

为方根信号的平均值，

为方根信号的方差。

进一步地，在步骤S102中，对所述事件信号x(τ)进行S时频变换和奇异值分解的具体计算步骤如下：

S201：事件信号x(τ)进行S变换，得到时频信息矩阵S(t,f)；其计算公式如式(3)所示：

式(3)中f表示频率，t表示时间，S(t,f)为所得的时频信息矩阵；

S202：采用奇异值分解方法对S201中得到的时频信息矩阵S(t,f)进行奇异值分解，其计算公式如式(4)所示：

S(t,f)＝UΓV^T (4)

式(4)中，S(t,f)是M×N的矩阵，U和V是M×M和N×N的正交矩阵，Γ是M×N的矩阵，其形式可以表示为：

其中，Λ_r×r(r≤min(M,N))是对角矩阵，形式如下：

Λ_r×r＝diag(σ₁,σ₂,...,σ_r)

σ_i(i＝1,2,...,r)为S(t,f)的分解后的奇异值，且满足σ₁≥σ₂≥...≥σ_r≥0。

进一步地，在步骤S102中，使用k-medoids算法确定有效奇异值的具体步骤为：

S301：从S202所得的奇异值σ_i(i＝1,2,...,r)中，选择o₁＝σ₁和o₂∈{σ₂,σ₃,...σ_r}作为k-medoids算法的中心；

S302：按公式(5)分别计算剩下的奇异值到中心o₁和o₂的距离

式(5)中，s_j∈{σ₂,σ₃,...,σ_r}且s_j≠o₂；若d₁(s_j,o₁)≤d₂(s_j,o₂)，将s_j分配到以o₁为中心的聚类C₁中；否则将s_j分配到以o₂为中心的聚类C₂中；

S303：从聚类C₂中，随机选择一个值o₂′，作为聚类中心，聚类C₁中的中心为o₁′＝o₁；按S302中方法计算新的以o₁′和o₂′聚类C₁′和C₂′；并计算如下值：

其中，若J<0，用o₂′作为新的聚类中心代替o₂；否则，保留原始的中心o₂；

S304：循环步骤S302和S303，直到新生成的聚类中心不再变化；包含σ₁的聚类是所需要的，获取其聚类的个数为n(1≤n≤r-1)，获得的聚类的个数即为有效奇异值的个数，即最终获得有效奇异值为σ′_i(i＝1,2...,n)；所述有效奇异值σ′_i(i＝1,2...,n)为步骤S202中中得到的奇异值σ_i(i＝1,2,...,r)的有效部分。

进一步地，在步骤S102中，所述重构信号x′(τ)的计算具体步骤为：

S401：根据获取的n个有效奇异值σ_i，按照公式(6)得到重构的时频矩阵

S′(t,f)＝UΓ′V^T (6)

其中Γ′为

Λ′_n×n＝diag(σ₁,σ₂,...,σ_n)

S402：使用S逆变换获得重构的时域信号x′(τ)，其表达式如式(7)所示。

进一步地，在步骤S102中，尖峰信号的检测具体包括以下过程：

S501：选取重构的时域信号x′(τ)在80-500Hz内的功率P_80-500；

S502：采用箱形图对所述功率P_80-500进行处理，得到上四分位和下四分位值Q₁和Q₃；所述功率下界异常阈值Threshold的计算公式如式(8)所示：

Threshold＝Q₁-1.5(Q₃-Q₁) (8)

S503：对重构信号x′(τ)进行阈值分割，高于阈值Threshold的重构信号所对应的疑似高频振荡节律为非尖峰信号，进行保留；低于阈值Threshold的重构信号所对应的疑似高频振荡节律为尖峰，并将之从疑似高频振荡节律中去除，从而得到去除尖峰活动的高精度高频振荡节律。

本发明基于其技术方案所具有的有益效果是：提高癫痫脑电信号的高频振荡节律的检测精度，帮助医生进行癫痫诊断和致痫灶切除。

附图说明

下面将结合附图及实例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于时频分析的自动分离高频振荡节律中尖峰活动方法的流程图；

图2是本发明实施例中高频振荡节律和尖峰的时域波形示意图；

图3是本发明实施例中高频振荡节律和尖峰的时频示意图。

图4是本发明实施例中高频振荡节律和尖峰的频谱示意图。

图5是本发明实施例中高频振荡节律和尖峰所对应的0-500Hz信号重构后的频谱示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种基于时频分析的自动分离高频振荡节律中尖峰活动方法。

请参考图1，图1是本发明实施例中一种基于时频分析的自动分离高频振荡节律中尖峰活动方法的流程图，具体包括以下步骤：

S101：采集癫痫病人的原始脑电信号，并对所述原始脑电信号进行50Hz倍频陷波和80-500Hz带通滤波，得到去除50Hz倍频并且频率在80-500Hz的预处理后的脑电信号，并提取所述预处理后的脑电信号中的疑似高频振荡节律所对应的原始信号；具体如下：

首先对所述预处理后的脑电信号加长度为L的时间窗，得到处理后的脑电信号；其中，L为预设值；

然后对所述处理后的脑电信号使用均方根算法，得到滤波后的癫痫脑电的均方根信号RMS(t)，并记录其持续时间T及波峰个数m；

最后设定均方根信号阈值Threshold_RMS，并判断条件RMS(t)>Threshold_RMS&&T>30ms&&m≥6是否成立？若是，则认为RMS(t)为疑似高频振荡节律，进行提取，从而得到疑似高频振荡节律所对应的原始信号；反之则认为RMS(t)为非高频振荡节律，予以遗弃；

S102：对所述疑似高频振荡节律所对应的原始信号进行尖峰信号检测，并将检测到的尖峰信号进行剔除，进而得到高精度的高频振荡节律；具体如下：

首先采用时频分析方法对所述疑似高频振荡节律所对应的原始信号进行50Hz倍频陷波和0-500Hz带通滤波，得到去除50Hz倍频并且频率在0-500Hz的事件信号x(τ)；其中，所述疑似高频振荡节律所对应的原始信号为步骤S101中的所述原始脑电信号；

对x(τ)进行S变换得到其时频信息S(t,f)，并对时频信息S(t,f)进行奇异值分解，得到分解后的奇异值；再使用k-mediods算法对所述分解后的奇异值进行处理，得到有效的奇异值；结合有效奇异值，对S(t,f)进行逆S变换，得到重构后的信号x′(τ)；

最后设定功率下界异常阈值Threshold，对x′(τ)进行阈值处理，检测出尖峰信号并剔除，进而得到高精度的高频振荡节律。

在步骤S101中，获取疑似高频振荡节律具体过程如下：

所述均方根信号RMS(t)的计算公式如式(1)所示：

式(1)中，s(k)为所述处理后的脑电信号，N_lL为所述处理后的脑电信号的长度；所述均方根阈值Threshold_RMS的设定公式如式(2)所示：

Threshold_RMS＝3*SD+E (2)

式(2)中，

为方根信号的平均值，

为方根信号的方差。

在步骤S102中，对所述事件信号x(τ)进行S时频变换和奇异值分解的具体计算步骤如下：

S201：事件信号x(τ)进行S变换，得到时频信息矩阵S(t,f)；其计算公式如式(3)所示：

式(3)中f表示频率，t表示时间，S(t,f)为所得的时频信息矩阵；

S202：采用奇异值分解方法对S201中得到的时频信息矩阵S(t,f)进行奇异值分解，其计算公式如式(4)所示：

S(t,f)＝UΓV^T (4)

式(4)中，S(t,f)是M×N的矩阵，U和V是M×M和N×N的正交矩阵，Γ是M×N的矩阵，其形式可以表示为：

其中，Λ_r×r(r≤min(M,N))是对角矩阵，形式如下：

Λ_r×r＝diag(σ₁,σ₂,...,σ_r)

σ_i(i＝1,2,...,r)为S(t,f)的分解后的奇异值，且满足σ₁≥σ₂≥...≥σ_r≥0。

在步骤S102中，使用k-medoids算法确定有效奇异值的具体步骤为：

S301：从S202所得的奇异值σ_i(i＝1,2,...,r)中，选择o₁＝σ₁和o₂∈{σ₂,σ₃,...σ_r}作为k-medoids算法的中心；

S302：按公式(5)分别计算剩下的奇异值到中心o₁和o₂的距离

式(5)中，s_j∈{σ₂,σ₃,...,σ_r}且s_j≠o₂；若d₁(s_j,o₁)≤d₂(s_j,o₂)，将s_j分配到以o₁为中心的聚类C₁中；否则将s_j分配到以o₂为中心的聚类C₂中；

S303：从聚类C₂中，随机选择一个值o₂′，作为聚类中心，聚类C₁中的中心为o₁′＝o₁；按S302中方法计算新的以o₁′和o₂′聚类C₁′和C₂′；并计算如下值：

其中，若J<0，用o₂′作为新的聚类中心代替o₂；否则，保留原始的中心o₂；

S304：循环步骤S302和S303，直到新生成的聚类中心不再变化；包含σ₁的聚类是所需要的，获取其聚类的个数为n(1≤n≤r-1)，获得的聚类的个数即为有效奇异值的个数，即最终获得有效奇异值为σ′_i(i＝1,2...,n)；所述有效奇异值σ′_i(i＝1,2...,n)为步骤S202中中得到的奇异值σ_i(i＝1,2,...,r)的有效部分。

在步骤S102中，所述重构信号x′(τ)的计算具体步骤为：

S401：根据获取的n个有效奇异值σ_i，按照公式(6)得到重构的时频矩阵

S′(t,f)＝UΓ′V^T (6)

其中Γ′为

Λ′_n×n＝diag(σ₁,σ₂,...,σ_n)

S402：使用S逆变换获得重构的时域信号x′(τ)，其表达式如式(7)所示。

在步骤S102中，尖峰信号的检测具体包括以下过程：

S501：选取重构的时域信号x′(τ)在80-500Hz内的功率P_80-500；

S502：采用箱形图对所述功率P_80-500进行处理，得到上四分位和下四分位值Q₁和Q₃；所述功率下界异常阈值Threshold的计算公式如式(8)所示：

Threshold＝Q₁-1.5(Q₃-Q₁) (8)

S503：对重构信号x′(τ)进行阈值分割，高于阈值Threshold的重构信号所对应的疑似高频振荡节律为非尖峰信号，进行保留；低于阈值Threshold的重构信号所对应的疑似高频振荡节律为尖峰，并将之从疑似高频振荡节律中去除，从而得到去除尖峰活动的高精度高频振荡节律。

请参考图2，图2是高频振荡节律和尖峰的时域波形示意图。

请参考图3，图3是在通过S变换方法获得高频振荡节律(左图)和尖峰(右图)的时频信息，从图3可以看出，高频振荡节律和尖峰活动的时频图较为相似，难以区分出来。

请参考图4，图4是高频振荡节律与尖峰活动的频谱图。

请参考图5，图5是高频振荡节律(左图)和尖峰(右图)所对应的重构信号的频谱图，从中可以看出重构后的信号中，高频振荡节律所对应的重构信号的幅值在80-500Hz频段内明显高于尖峰所对应的重构信号。因此，在得到重构信号的80-500Hz频段内的功率时，采用箱型图得到上四分位和下四分位Q₁和Q₃，然后设置功率阈值Threshold＝Q₁-1.5(Q₃-Q₁)，低于重构后的时域信号所对应的疑似高频振荡节律为尖峰，并将之从疑似高频振荡节律中去除，从而得到去除尖峰活动的疑似高频振荡节律。

本发明的有益效果是：提高癫痫脑电信号的高频振荡节律的检测精度，帮助医生进行癫痫诊断和致痫灶切除。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于：包括以下步骤：

S101：采集癫痫病人的原始脑电信号，并对所述原始脑电信号进行50Hz倍频陷波和80-500Hz带通滤波，得到去除50Hz倍频并且频率在80-500Hz的预处理后的脑电信号，并提取所述预处理后的脑电信号中的疑似高频振荡节律所对应的原始信号；具体如下：

首先对所述预处理后的脑电信号加长度为L的时间窗，得到处理后的脑电信号；其中，L为预设值；

然后对所述处理后的脑电信号使用均方根算法，得到滤波后的癫痫脑电的均方根信号RMS(t)，并记录其持续时间T及波峰个数m；

最后设定均方根信号阈值Threshold_RMS，并判断条件RMS(t)>Threshold_RMS&&T>30ms&&m≥6是否成立，若是，则认为RMS(t)为疑似高频振荡节律，进行提取，从而得到疑似高频振荡节律所对应的原始信号；反之则认为RMS(t)为非高频振荡节律，予以遗弃；

S102：对所述疑似高频振荡节律所对应的原始信号进行尖峰信号检测，并将检测到的尖峰信号进行剔除，进而得到高精度的高频振荡节律；具体如下：

首先采用时频分析方法对所述疑似高频振荡节律所对应的原始信号进行50Hz倍频陷波和0-500Hz带通滤波，得到去除50Hz倍频并且频率在0-500Hz的事件信号x(τ)；其中，所述疑似高频振荡节律所对应的原始信号为步骤S101中的所述原始脑电信号；

对x(τ)进行S变换得到其时频信息S(t,f)，并对时频信息S(t,f)进行奇异值分解，得到分解后的奇异值；再使用k-mediods算法对所述分解后的奇异值进行处理，得到有效的奇异值；结合有效奇异值，对S(t,f)进行逆S变换，得到重构后的信号x′(τ)；

最后设定功率下界异常阈值Threshold，对x′(τ)进行阈值处理，检测出尖峰信号并剔除，进而得到高精度的高频振荡节律。

2.根据权利要求1所述的一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于：步骤S101中，获取疑似高频振荡节律具体过程如下：

所述均方根信号RMS(t)的计算公式如式(1)所示：

式(1)中，s(k)为所述处理后的脑电信号，N_lL为所述处理后的脑电信号的长度；所述均方根阈值Threshold_RMS的设定公式如式(2)所示：

Threshold_RMS＝3*SD+E (2)

式(2)中，

为方根信号的平均值，

为方根信号的方差。

3.根据权利要求1所述的一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于，步骤S102中，对所述事件信号x(τ)进行S时频变换和奇异值分解的具体计算步骤如下：

S201：对所述事件信号x(τ)进行S时频变换，得到时频信息矩阵S(t,f)；其计算公式如式(3)所示：

式(3)中f表示频率，t表示时间，S(t,f)为所得的时频信息矩阵；

S202：采用奇异值分解方法对S201中得到的时频信息矩阵S(t,f)进行奇异值分解，其计算公式如式(4)所示：

S(t,f)＝UΓV^T (4)

式(4)中，S(t,f)是M×N的矩阵，U和V是M×M和N×N的正交矩阵，Γ是M×N的矩阵，其形式可以表示为：

其中，Λ_r×r(r≤min(M,N))是对角矩阵，形式如下：

Λ_r×r＝diag(σ₁,σ₂,…,σ_r)

σ_i(i＝1,2,...,r)为S(t,f)的分解后的奇异值，且满足σ₁≥σ₂≥...≥σ_r≥0。

4.根据权利要求3所述的一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于：在步骤S102中，使用k-medoids算法确定有效奇异值的具体步骤为：

S301：从S202所得的奇异值σ_i(i＝1,2,...,r)中，选择o₁＝σ₁和o₂∈{σ₂,σ₃,...σ_r}作为k-medoids算法的中心；

S302：按公式(5)分别计算剩下的奇异值到中心o₁和o₂的距离

式(5)中，s_j∈{σ₂,σ₃,...,σ_r}且s_j≠o₂；若d₁(s_j,o₁)≤d₂(s_j,o₂)，将s_j分配到以o₁为中心的聚类C₁中；否则将s_j分配到以o₂为中心的聚类C₂中；

S303：从聚类C₂中，随机选择一个值o₂′，作为聚类中心，聚类C₁中的中心为o₁′＝o₁；按S302中方法计算新的以o₁′和o₂′聚类C₁′和C₂′；并计算如下值：

其中，若J<0，用o₂′作为新的聚类中心代替o₂；否则，保留原始的中心o₂；

S304：循环步骤S302和S303，直到新生成的聚类中心不再变化；包含σ₁的聚类是所需要的，获取其聚类的个数为n(1≤n≤r-1)，获得的聚类的个数即为有效奇异值的个数，即最终获得有效奇异值为σ_i′(i＝1,2...,n)；所述有效奇异值σ_i′(i＝1,2...,n)为步骤S202中中得到的奇异值σ_i(i＝1,2,...,r)的有效部分。

5.根据权利要求4所述的一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于：步骤S102中，所述重构信号x′(τ)的计算具体步骤为：

S401：根据获取的n个有效奇异值σ_i，按照公式(6)得到重构的时频矩阵

S′(t,f)＝UΓ′V^T (6)

其中Γ′为

Λ′_n×n＝diag(σ₁,σ₂,...,σ_n)

S402：使用S逆变换获得重构的时域信号x′(τ)，其表达式如式(7)所示：

6.根据权利要求1所述的一种基于时频分析的自动分离高频振荡节律中尖峰活动方法，其特征在于，步骤S102中，尖峰信号的检测具体包括以下过程：

S501：选取重构的时域信号x′(τ)在80-500Hz内的功率P_80-500；

S502：采用箱形图对所述功率P_80-500进行处理，得到上四分位和下四分位值Q₁和Q₃；所述功率下界异常阈值Threshold的计算公式如式(8)所示：

Threshold＝Q₁-1.5(Q₃-Q₁) (8)

S503：对重构信号x′(τ)进行阈值分割，高于阈值Threshold的重构信号所对应的疑似高频振荡节律为非尖峰信号，进行保留；低于阈值Threshold的重构信号所对应的疑似高频振荡节律为尖峰，并将之从疑似高频振荡节律中去除，从而得到去除尖峰活动的高精度高频振荡节律。

Images