CN113558643A - 基于vmd和nltwsvm的多特征癫痫信号分类方法 - Google Patents

Abstract

本发明属于涉及一种基于VMD和NLTWSVM的多特征癫痫信号分类方法。包括：对脑电信号数据进行预处理；脑电信号数据包括健康人脑电数据集、癫痫患者发作间期脑电数据集和癫痫患者发作期脑电数据集；对脑电信号进行变分模态分解，得到相对应的本征模态函数；从得到的一系列本征模态函数中提取多特征参数；将提取的多特征参数作为特征组合输入到非线性双子支持向量机中进行训练和分类。针对癫痫脑电信号具有非平稳、非线性等特点，利用变分模态分解时频分析方法获取信号分量，对分量进行多特征参数提取，结合非线性双子支持向量机机器学习方法，实现正确分类癫痫脑电信号和非癫痫脑电信号。

Description

基于VMD和NLTWSVM的多特征癫痫信号分类方法

技术领域

本发明属于非平稳非线性信号的时频分析、模式分类和机器学习，属于生物医学工程信号处理和模式识别技术领域，具体涉及一种基于VMD和NLTWSVM的多特征癫痫信号分类方法。

背景技术

癫痫是大脑神经元突发性超同步化异常放电导致短暂大脑功能障碍的一种慢性非传染性神经***脑疾病，据世界卫生组织统计全世界目前大约有5000万癫痫患者，其中大约80％的癫痫患者生活在中低收入国家，癫痫疾病对癫痫患者的生活、学习、工作和精神造成了严重影响，甚至危害癫痫患者的生命安全。目前对癫痫患者进行神经电生理检查时，脑电图(Electroencephalography，EEG)作为癫痫诊断的“金标准”，通过放置于大脑皮层处的干电极记录患者脑电信号，医师利用肉眼观测脑电信号数据判断脑电图中是否出现表征大脑过度放电的棘波、尖波、棘(尖)-慢综合波和其他发作性节律波。但是仅凭专家肉眼观测判断非常耗时耗力，并且在观测过程中可能存在由于人为主观因素导致造成误诊的情况。因此，癫痫脑电信号和非癫痫脑电信号进行分类识别对癫痫发作预测和临床诊疗领域具有重要意义。

脑电信号是典型的非线性、非平稳的生物弱信号。脑电信号的分析方法有时域、频域、时频域和非线性动力学。其中时频域方法中由Huang等人提出的经验模态分解(Empirical Mode Decomposition，EMD)方法在脑电信号分析领域应用广泛，但是此方法由于通过极值点的包络进行多次递归分解后包络估计误差逐渐被放大，导致相近的频率分量无法被正确分离而存在模态混叠和端点效应等问题，使信号在处理过程中丢失部分细节信息，影响信号的完整性和分析的准确性。针对EMD存在模态混叠问题，Wu和Huang提出了改进的利用噪声辅助EMD的集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)方法，该方法利用白噪声在频域内功率谱密度分布均匀的特点，消除了模态混叠现象并且提高了算法的抗噪性，但是此方法增加了计算量，同时会分解得到超出信号真实组成的多个分量。

发明内容

本发明所要解决的技术问题在于提供一种基于VMD和NLTWSVM的多特征癫痫信号分类方法，针对癫痫脑电信号具有非平稳、非线性等特点，利用变分模态分解时频分析方法获取信号分量，对分量进行多特征参数提取，结合非线性双子支持向量机机器学习方法，实现正确分类癫痫脑电信号和非癫痫脑电信号。

本发明是这样实现的，

一种基于VMD和NLTWSVM的多特征癫痫信号分类方法，该方法包括：

对脑电信号数据进行预处理；脑电信号数据包括健康人脑电数据集、癫痫患者发作间期脑电数据集和癫痫患者发作期脑电数据集；

对脑电信号x(t)进行变分模态分解，得到相对应的本征模态函数u_k(t)，k为K里某一个层数的代表，其中K为本征模态函数的阶数；

从得到的一系列本征模态函数中提取多特征参数；

将提取的多特征参数作为特征组合输入到非线性双子支持向量机中进行训练和分类。

进一步地，预处理包括：将数据集中每组脑电信号前23s脑电数据以1s为一个单位划分为23个数据片段样本，将每个数据集随机选取50组脑电数据并对其分段获得的数据片段样本作为分类器的训练样本，剩余每个数据集的数据片段样本作为分类器的测试样本。

进一步地，对脑电信号x(t)进行变分模态分解包括：

建立变分模型；

引入二次惩罚因子α和拉格朗日乘法算子λ(t)，使约束性的变分问题转变为非约束性的变分问题，构造扩展拉格朗日函数；

利用乘法算子交替方向法求解变分问题，通过交替更新

和λⁿ⁺¹(t)，寻求扩展拉格朗日函数的“鞍点”，此“鞍点”为变分模型的最优解；

其中，建立变分模型包括：

设u_k(t)为IMF分量，对每个IMF分量u_k(t)作Hilbert变换，构造解析信号，得到其单边频谱：

对各IMF分量的解析信号混合一预估中心频率

将每个IMF分量的频谱调制到相对应的基频带：

计算式(2)所表示的解调信号梯度的平方L²范数，估计出各个IMF分量的信号带宽，引入约束条件，构造如式(3)的变分模型：

式中：δ(t)为冲激函数，*为卷积符号，{u_k}＝{u₁,u₂,…,u_K}为分解得到的各IMF分量，{ω_k}＝{ω₁,ω₂,…,ω_K}为各IMF分量的中心频率，K为IMF分量的数目，

扩展拉格朗日函数的表达式如下：

利用乘法算子交替方向法求解变分问题采用的迭代公式如下：

迭代停止条件为：

式中，τ为拉格朗日乘子更新参数，ε为收敛容限。

进一步地，对本征模态函数提取多特征参数时，脑电信号x(t)通过VMD分解后得到的信号u_k(t)为一个长度为L的时间序列，提取的特征通过以下表达式求取：

最大值：Max＝max(u_k(t)) (9)

最小值：Min＝min(u_k(t)) (10)

平均值：

方差：

偏度：

峰度：

变化系数：

波动指数：

进一步地，将提取的多特征参数作为特征组合输入到非线性双子支持向量机中进行训练和分类具体包括：

步骤5.1：设有训练样本，其中有p个正类样本和q个负类样本，表达式如下：

T＝{(x₁,+1),(x₂,+1),…(x_p,+1),(x_p+1,-1),(x_p+2,-1),…(x_p+q,-1)} (17)

其中x_i∈Rⁿ,i＝1,2,…,p+q；

步骤5.2：对于非线性双子支持向量机，通过引入核函数将原空间线性不可分的样本映射到线性可分的高维空间中，非线性双子支持向量机寻求的仍然是两个非平行的决策超平面，如下：

其中C^T＝[A；B]^T∈R^n×l，A＝(x₁,x₂,…,x_p)^T∈R^p×n，B＝(x_p+1,x_p+2,…,x_p+q)^T∈R^q×n，l＝p+q；

步骤5.3：双子支持向量机要求每条超平面满足距离对应类别样本点最小，以及足够远离另一类样本点的条件，据此构造两个凸二次规划问题如下：

和

其中e₊和e_-是元素全为一的列向量，ξ₊和ξ_-为松弛向量，c_i,i＝1,2是惩罚参数；

步骤5.4：式(19)和式(20)的对偶问题分别为

和

其中S＝[K(A,C^T)e₊]，R＝[K(B,C^T)e_-]；

步骤5.5：通过求解式(21)和式(22)，原始问题的解通过求解以下式子得到：

对于一个新的输入x∈Rⁿ，它的类别判断为

其中|·|是点x到超平面K(x^T,C^T)u_k+b_k＝0,k＝-,+的垂直距离。

本发明与现有技术相比，有益效果在于：

本发明采用了变分模态分解方法，可对信号多尺度时频分析分解，有效避免了模态混叠问题，具有较好的复杂信号分解精度及较好的抗噪声干扰等优点；

利用统计学可以客观表征信号总体数量特征和数量关系，对各IMF分量应用统计学方法提取共8个特征参数并构造箱式图，特征组合方式有效，可准确表征癫痫脑电信号和非癫痫脑电信号的差异性，有利于数据量大时快速计算特征值；使用非线性双子支持向量机进行训练和分类，在分类效果和计算速度上更加具有优势。本发明使用的双子支持向量机作为支持向量机(Support Vector Machine，SVM)的延伸，与传统SVM方法最大的不同在于，双子支持向量机针对每一类别构造相对应的支持超平面，让每一个超平面尽量靠近对应类别的训练点，而远离另一类训练点，而不限定于超平面是否平行和间隔大小。双子支持向量机在分类效果和计算速度上更加具有优势。

附图说明

图1为本发明方法流程图；

图2为变分模态分解模型求解流程图；

图3为数据库中具有代表性的正常、癫痫患者发作间期和癫痫患者发作期三种对象状态的原始脑电信号图；

图4为健康脑电A6、癫痫发作间期脑电C3和癫痫发作期脑电E10具有代表性的脑电数据片段样本及其VMD分解后各阶IMF分量图；

图5为IMF1、IMF2、IMF3、IMF4在正常、癫痫发作间期、癫痫发作期5个数据集中的最大值、最小值、平均值、方差、偏度、峰度、变化系数、波动指数箱式图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，一种基于VMD和NLTWSVM的多特征癫痫信号分类方法，其中，VMD为Variational Mode Decomposition的缩写，中文为变分模态分解，NLTWSVM为Non-linearTwin Support Vector Machine，为非线性双子支持向量机的缩写，本发明所述方法包括以下步骤：

步骤1：选取德国波恩癫痫研究室临床采集的脑电数据库，该数据库包括5组脑电数据集，每组各有100段时长为23.6s的脑电信号，每段信号采集4096个点，采样频率为173.6Hz。数据集A、B分别是健康人在睁眼和闭眼时采集到的大脑皮层脑电数据，数据集C、D分别是癫痫患者发作间期在致痫灶外和致痫灶内采集到的颅内脑电数据，数据集E是癫痫患者发作期在致痫灶内获得的颅内脑电数据。实验选取5种脑电数据各100组，将每组前23s脑电数据以1s为一个单位划分为23个数据片段样本，每个样本长度为173点数据，对每个数据集随机选取50组脑电数据并对其分段获得了1150个数据片段样本，使用每个数据片段样本计算所得的特征值作为分类器的训练样本，剩余每个数据集的数据片段样本对应特征值用于测试。

步骤2：对预处理后的脑电信号x(t)进行变分模态分解，得到本征模态函数u_k(t)，K为变分模态函数的阶数；

步骤2.1通过以下步骤建立变分模型：

步骤2.1.1、对每个IMF分量u_k(t)作Hilbert变换，构造解析信号，得到其单边频谱：

步骤2.1.2、对各IMF分量的解析信号混合一预估中心频率

将每个IMF分量的频谱调制到相对应的基频带：

步骤2.1.3、计算式(2)所表示的解调信号梯度的平方L²范数，估计出各个IMF分量的信号带宽，引入约束条件，构造如下形式的变分模型：

式中：δ(t)为冲激函数，*为卷积符号，{u_k}＝{u₁,u₂,…,u_K}为分解得到的各IMF分量，{ω_k}＝{ω₁,ω₂,…,ω_K}为各IMF分量的中心频率，

步骤2.2、引入二次惩罚因子α和拉格朗日乘法算子λ(t)，使约束性的变分问题转变为非约束性的变分问题。构造的扩展拉格朗日函数表达式如下：

步骤2.3、利用乘法算子交替方向法(Alternate Direction Method ofMultipliers，ADMM)求解上述变分问题，通过交替更新

和λⁿ⁺¹(t)，寻求扩展拉格朗日函数的“鞍点”，此“鞍点”即为变分模型的最优解。迭代公式如下：

迭代停止条件为：

式中，τ为拉格朗日乘子更新参数，ε为收敛容限。

公式(5)所述模态u_k更新求解方法包括如下步骤：

1将式(5)改写为如下等值表达式：

2基于Parseval/Plancherel傅里叶等距变换，将式(25)转换到频域：

3、将式(26)中带有惩罚因子项中的ω用ω-ω_k代替，可得：

4、由于实信号具有埃尔米特对称性，可将式(27)改为非负频域上的半空间积分形式：

5、把式(28)以一个二次优化问题进行求解，可得：

公式(6)所述中心频率ω_k更新求解方法包括如下步骤：

1、将式(6)改写为如下等值表达式：

2、根据同样的过程，将式(30)进行傅里叶变换到频域，并最终变换成非负频域上的半空间积分形式：

3、求解式(31)的二次优化问题，求得中心频率的更新方法如下：

因此，如图2所示，使用乘法算子交替方向法具体求解包括如下步骤：

1、初始化

和n为0；

2、根据n←n+1更新n；

3、根据式(29)和式(32)迭代更新{u_k}和{ω_k}；

4、根据

更新

5、满足迭代停止条件(8)则停止迭代，否则返回步骤2。

步骤3：从得到的本征模态函数中提取多特征参数，假设脑电信号x(t)通过VMD分解后得到的信号u_k(t)为一个长度为L的时间序列，提取的特征可通过以下表达式求取：

最大值：Max＝max(u_k(t)) (9)

最小值：Min＝min(u_k(t)) (10)

平均值：

方差：

偏度：

峰度：

变化系数：

波动指数：

步骤4：将特征提取得到的8种特征组合输入到非线性双子支持向量机进行训练，非线性双子支持向量机具体实现步骤包括：

设有训练样本，其中有p个正类样本和q个负类样本，表达式如下：

T＝{(x₁,+1),(x₂,+1),…(x_p,+1),(x_p+1,-1),(x_p+2,-1),…(x_p+q,-1)} (17)

其中x_i∈Rⁿ,i＝1,2,…,p+q；

对于非线性双子支持向量机，通过引入核函数将原空间线性不可分的样本映射到线性可分的高维空间中，非线性双子支持向量机寻求的仍然是两个非平行的决策超平面，如下：

其中C^T＝[A；B]^T∈R^n×l，A＝(x₁,x₂,…,x_p)^T∈R^p×n，B＝(x_p+1,x_p+2,…,x_p+q)^T∈R^q×n，l＝p+q；

双子支持向量机要求每条超平面满足距离对应类别样本点最小，以及足够远离另一类样本点的条件，据此构造两个凸二次规划问题如下：

和

其中e₊和e_-是元素全为一的列向量，ξ₊和ξ_-为松弛向量。c_i,i＝1,2是惩罚参数；

式(19)和式(20)的对偶问题分别为

和

其中S＝[K(A,C^T)e₊]，R＝[K(B,C^T)e_-]；

通过求解式(21)和式(22)，原始问题的解可以通过求解以下式子得到。

对于一个新的输入x∈Rⁿ，它的类别判断为

其中|·|是点x到超平面K(x^T,C^T)u_k+b_k＝0,k＝-,+的垂直距离。

本发明实验数据来自德国波恩癫痫脑电数据库。正常、癫痫患者发作间期和癫痫患者发作期三种对象状态对应的数据集中具有代表性的原始脑电信号图如图3所示。本方法对数据片段样本进行4层VMD分解，得到各阶IMF分量。健康脑电A6、癫痫发作间期脑电C3和癫痫发作期脑电E10具有代表性的脑电数据片段样本及其VMD分解后各阶IMF分量图分别如图4所示。从图4可以看出，VMD分解所得各阶IMF分量代表着原始脑电信号中不同特征成分，癫痫发作期脑电比健康状态和癫痫发作间期脑电幅度变化剧烈。

通过VMD分解得到5种脑电数据集脑电数据片段样本的各阶IMF分量后，对每个IMF分量分别计算其相应的最大值、最小值、平均值、方差、偏度、峰度、变化系数和波动指数8个特征。VMD分解所得IMF1、IMF2、IMF3、IMF4在正常、癫痫发作间期、癫痫发作期5个数据集中的8个特征指数箱式图如图5所示。从图5可看出，本方法所提取的每个特征指数中，癫痫发作期的各阶IMF分量相较于正常状态和癫痫发作间期的IMF分量基本上存在显著统计差异，这些统计差异为准确分类癫痫脑电信号和非癫痫脑电信号奠定了良好的基础。本方法使用的非线性双子支持向量机分类准确率可以达到98.59％，所提方法适用于对癫痫脑电信号和非癫痫脑电信号进行分类检测。

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

Claims (5)

1.一种基于VMD和NLTWSVM的多特征癫痫信号分类方法，其特征在于，该方法包括：

对脑电信号数据进行预处理；脑电信号数据包括健康人脑电数据集、癫痫患者发作间期脑电数据集和癫痫患者发作期脑电数据集；

对脑电信号x(t)进行变分模态分解，得到相对应的本征模态函数u_k(t)，k为K里某一个层数的代表，其中K为本征模态函数的阶数；

从得到的一系列本征模态函数中提取多特征参数；

将提取的多特征参数作为特征组合输入到非线性双子支持向量机中进行训练和分类。

2.按照权利要求1所述的方法，其特征在于，预处理包括：将数据集中每组脑电信号前23s脑电数据以1s为一个单位划分为23个数据片段样本，将每个数据集随机选取50组脑电数据并对其分段获得的数据片段样本作为分类器的训练样本，剩余每个数据集的数据片段样本作为分类器的测试样本。

3.按照权利要求1所述的方法，其特征在于，对脑电信号x(t)进行变分模态分解包括：

建立变分模型；

引入二次惩罚因子α和拉格朗日乘法算子λ(t)，使约束性的变分问题转变为非约束性的变分问题，构造扩展拉格朗日函数；

利用乘法算子交替方向法求解变分问题，通过交替更新

和λⁿ⁺¹(t)，寻求扩展拉格朗日函数的“鞍点”，此“鞍点”为变分模型的最优解；

其中，建立变分模型包括：

设u_k(t)为IMF分量，对每个IMF分量u_k(t)作Hilbert变换，构造解析信号，得到其单边频谱：

对各IMF分量的解析信号混合一预估中心频率

将每个IMF分量的频谱调制到相对应的基频带：

计算式(2)所表示的解调信号梯度的平方L²范数，估计出各个IMF分量的信号带宽，引入约束条件，构造如式(3)的变分模型：

式中：δ(t)为冲激函数，*为卷积符号，{u_k}＝{u₁,u₂,…,u_K}为分解得到的各IMF分量，{ω_k}＝{ω₁,ω₂,…,ω_K}为各IMF分量的中心频率，K为IMF分量的数目，

扩展拉格朗日函数的表达式如下：

利用乘法算子交替方向法求解变分问题采用的迭代公式如下：

迭代停止条件为：

式中，τ为拉格朗日乘子更新参数，ε为收敛容限。

4.按照权利要求1所述的方法，其特征在于，对本征模态函数提取多特征参数时，脑电信号x(t)通过VMD分解后得到的信号u_k(t)为一个长度为L的时间序列，提取的特征通过以下表达式求取：

最大值：Max＝max(u_k(t)) (9)

最小值：Min＝min(u_k(t)) (10)

平均值：

方差：

偏度：

峰度：

变化系数：

波动指数：

5.按照权利要求1所述的方法，其特征在于，将提取的多特征参数作为特征组合输入到非线性双子支持向量机中进行训练和分类具体包括：

步骤5.1：设有训练样本，其中有p个正类样本和q个负类样本，表达式如下：

T＝{(x₁,+1),(x₂,+1),…(x_p,+1),(x_p+1,-1),(x_p+2,-1),…(x_p+q,-1)} (17)

其中x_i∈Rⁿ,i＝1,2,…,p+q；

步骤5.2：对于非线性双子支持向量机，通过引入核函数将原空间线性不可分的样本映射到线性可分的高维空间中，非线性双子支持向量机寻求的仍然是两个非平行的决策超平面，如下：

其中C^T＝[A；B]^T∈R^n×l，A＝(x₁,x₂,…,x_p)^T∈R^p×n，B＝(x_p+1,x_p+2,…,x_p+q)^T∈R^q×n，l＝p+q；

步骤5.3：双子支持向量机要求每条超平面满足距离对应类别样本点最小，以及足够远离另一类样本点的条件，据此构造两个凸二次规划问题如下：

和

其中e₊和e_-是元素全为一的列向量，ξ₊和ξ_-为松弛向量，c_i,i＝1,2是惩罚参数；

步骤5.4：式(19)和式(20)的对偶问题分别为

和

其中S＝[K(A,C^T)e₊]，R＝[K(B,C^T)e_-]；

步骤5.5：通过求解式(21)和式(22)，原始问题的解通过求解以下式子得到：

对于一个新的输入x∈Rⁿ，它的类别判断为

其中|·|是点x到超平面K(x^T,C^T)u_k+b_k＝0,k＝-,+的垂直距离。

Images