CN110575164B - 脑电信号伪迹去除方法及计算机可读存储介质 - Google Patents

脑电信号伪迹去除方法及计算机可读存储介质 Download PDF

Info

Publication number
CN110575164B
CN110575164B CN201910893087.2A CN201910893087A CN110575164B CN 110575164 B CN110575164 B CN 110575164B CN 201910893087 A CN201910893087 A CN 201910893087A CN 110575164 B CN110575164 B CN 110575164B
Authority
CN
China
Prior art keywords
signal
pure
target
modal
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910893087.2A
Other languages
English (en)
Other versions
CN110575164A (zh
Inventor
杨道国
耿道双
蔡苗
张国旗
郝卫东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guilin University of Electronic Technology
Original Assignee
Guilin University of Electronic Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guilin University of Electronic Technology filed Critical Guilin University of Electronic Technology
Priority to CN201910893087.2A priority Critical patent/CN110575164B/zh
Publication of CN110575164A publication Critical patent/CN110575164A/zh
Application granted granted Critical
Publication of CN110575164B publication Critical patent/CN110575164B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/369—Electroencephalography [EEG]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems
    • A61B5/7267—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems involving training the classification device
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Signal Processing (AREA)
  • Medical Informatics (AREA)
  • Psychiatry (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Evolutionary Computation (AREA)
  • Fuzzy Systems (AREA)
  • Mathematical Physics (AREA)
  • Psychology (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明提供了一种脑电信号伪迹去除方法及计算机可读存储介质，其中，脑电信号伪迹去除方法，包括：获取含冷痛信息的脑电信号作为目标信号，将目标信号进行迭代降噪以得到纯净模态分量；对纯净模态分量进行短时傅立叶时频分析以得到时频图，根据时频图确定有效特征，构造有效特征对应的模态分量以得到有效信号，从有效信号中筛选出纯净信号；将纯净信号的测试数据集输入到分类模型中，获取纯净信号的精度评价参数；纯净信号的检测精度评价参数满足精度需求确定纯净信号为合格信号。通过对带有冷痛信息的多通道脑电信号进行VMD分解和降噪，改善了模态混叠、端点效应的问题，提高了分解降噪的精度同时也提高了提取有效特征的精度。

Description

脑电信号伪迹去除方法及计算机可读存储介质
技术领域
本发明涉及技术领域，具体而言，涉及一种脑电信号伪迹去除方法及计算机可读存储介质。
背景技术
长期以来，急性疼痛作为一种诱发脑活动信号的方式被广泛研究，并作为一种临床诊断疼痛引起的疾病实验室方法。大多数研究疼痛都是经过多通道采集头皮EEG(脑电图)，然后经过信号处理来评估的。研究中发现，疼痛可以调节大脑皮层的同步性，在神经网络中的组成节点之间有明确的无节律的“振荡”。而且痛觉诱发EEG信号非平稳性强且信噪比低，鲁棒性难以掌握，从EEG信号中提取特征非常困难。特别是在急性冷压痛条件下，诱发的受试神经元活动激活的瞬间振荡，实验室中很难被捕捉到。一般对EEG的处理，都是使用Fourier进行时域到频域上的转换，进行时频分析。然而在相关技术中经常提到的频率，是平稳信号EEG处理中通过Fourier变换直接得到的频率与实际变换无关的参量，实质是信号的局部特性。但EEG的频率是随时间变化的，直接使用Fourier变换是不能很好地分析出局部变化特性，即不能动态地反应EEG的时频特性。很多研究采用增加电极导联数目来提高EEG的识别精度，然而通过研究表明，过多的电极则会采集到附近电极的脑区信号，产生电极间脑源信号的交叉混叠污染。
在实际中，冷痛刺激脑信号的检测容易受到噪声环境多变的影响，更加有效的办法则是采用高效的算法来实现信号特征的提取和降噪。目前常用的EEG信号特征提取和去伪方法有经验模态分解(EMD)、集成经验模态分解(EEMD)和小波变换(WT)等，这些方法抗噪性能差、误判率高、提取精度低。有些方法如EMD分解过程中容易出现模态混叠、伪分量以及端点效应等问题，即信号和噪声混叠在一个或多个固有模态函数(IMF)中或者分解出大量无意的分量。小波变换存在平移变化特性和频带混叠等缺陷，自适应能力差，时频分辨率低。不利于刺激诱发脑电信号瞬时频率频率的提取，以及降噪去伪的效果。且上述方法检测冷痛觉脑活动精确度不高，不利于实际应用和医学临床检测。
发明内容
本发明旨在至少解决现有技术中存在的技术问题之一。
本发明的第一个方面提供了一种脑电信号伪迹去除方法。
本发明的第二个方面提供了一种计算机可读存储介质。
鉴于上述，根据本发明的第一个方面，提出了一种脑电信号伪迹去除方法，包括：获取含冷痛信息的脑电信号作为目标信号，将目标信号进行迭代降噪以得到纯净模态分量；对纯净模态分量进行短时傅立叶时频分析以得到时频图，根据时频图确定有效特征，构造有效特征对应的模态分量以得到有效信号，从有效信号中筛选出纯净信号；将纯净信号的测试数据集输入到分类模型中，获取纯净信号的精度评价参数；纯净信号的精度评价参数满足精度需求确定纯净信号为合格信号。
在该技术方案中，获取脑电图中带有冷痛信息的多通道脑电信号，将含冷痛信息的脑电信号未做目标信号，通过对目标信号进行VMD(变分模态分解，Variational ModeDecomposition)分解实现对脑电信号中的去除眼电、心电、血流等伪迹信号迭代滤波，实现对脑电信号的迭代降噪得到纯净模态分量，纯净模态分量中存在多个模态分量，对纯净模态分量中的各个模态分量进行短时傅立叶时频分析可以得到每个模态分量的时频图，时频图可以反映出模态分量的时间和频率的关系，根据时频图中的时间范围和频率范围确定有效特征，有效特征为可以作为有效信号的模态分量，将有效特征构造成有效信号，并对有效信号再次进行滤波的筛选得到纯净信号，再将纯净信号的平均功率谱密度作为特征集投入到分类模型中进行检测，以得到用于判断信号去噪精度的精度评价参数，用户可以对各类的精度评价参数类型设置响应的精度需求，如果精度评价参数达到了设定的精度需求则将纯净信号作为合格信号输出，如果精度评价参数未达到设定的精度需求则将得到的纯净信号作为目标信号返回执行迭代降噪的步骤进行去除伪迹降噪，直至纯净信号满足精度需求，通过对脑电图中带有冷痛信息的多通道脑电信号进行VMD分解和降噪，改善了模态混叠、端点效应的技术问题，适合由外部刺激引起的脑神经活动检测，以及通过外部刺激诱发的非植入式头皮EEG(脑电图)的特征提取和降噪。提高了分解降噪的精度同时也提高了提取有效特征的精度，对诱发类脑电活动的临床检测和应用具有很好的技术价值。
另外，本发明提供的上述技术方案中的脑电信号伪迹去除方法还可以具有如下附加技术特征：
在上述任一技术方案中，将目标信号进行迭代降噪以得到纯净模态分量的步骤，具体包括：对冷痛信息的脑电信号进行分解确定目标模态分量；确定目标模态分量的数量和目标模态分量的中心频率；根据目标模态分量的数量和中心频率对目标模态分量进行优化得到优化模态分量；优化模态分量达到设定条件，确定优化模态分量为纯净模态分量；优化模态分量未达到设定条件，将优化模态分量作为目标模态分量返回执行确定目标模态分量的数量和目标模态分量的中心频率的步骤。
在该技术方案中，对多通道的脑电信号进行VMD分解以得到多个目标模态分量，确定目标模态分量的数量和中心频率并根据模态分量的数量和中心频率进行优化得到优化后的优化模态分量，并设定相应的设定条件，判断优化后的优化模态分量是否满足设定条件，如果满足设定条件则将优化模态分量作为纯净模态分量输出，如果未满足设定条件则将优化模态分量作为目标模态分量再次进行确定目标模态分量的数量和中心频率的步骤，直至得到的优化模态分量满足设定条件。
具体地，对每个目标模态分量利用Hilbert变换计算目标模态分量的频谱，具体公式如下：
Figure GDA0003531934520000031
Figure GDA0003531934520000032
其中，{uk}：＝{u1，…，uK}为所有目标模态分量，{ωk}：＝{ω1，…，ωK}为所有目标模态分量的中心频率集，uk代表第k个模态，ωk代表第k个模态的中心频率，
Figure GDA0003531934520000041
为所有目标模态分量的和。
将得到的频谱利用二次惩罚因子a和拉格朗日乘子λ(t)获取目标模态分量和其中心频率，从而保证信号重构之后的保真度以及降低噪声的效果。
其中增广拉格朗日表达式为：
Figure GDA0003531934520000042
对获取的目标模态分量和其中心频率通过交替方向乘子方法(ADMM)，进行交替更新，通过计算
Figure GDA0003531934520000043
和
Figure GDA0003531934520000044
来寻求增广拉格朗日表达式的‘鞍点’，‘鞍点’为拉格朗日表达式中既不是最大值也不是最小值的临界点，具体表达式为：
Figure GDA0003531934520000045
利用L2范数下的Parseval/Plancherel傅立叶等距法，将步骤上述表达式转换成谱域内求解，使每个目标模态分量都能获取最新的状态，转换后的表达式具体为：
Figure GDA0003531934520000046
将目标模态分量ω用ω-ωk替换第一项，根据Hermitain对称性质，变换为非负频率区间内的积分形式，具体表达式为：
Figure GDA0003531934520000047
所有的目标模态分量可以在频域中的解中获取，第k个目标模态分量更新的表达式为：
Figure GDA0003531934520000048
对目标模态分量个数k值与二次罚函数α进行定义，并初始化
Figure GDA0003531934520000049
Figure GDA00035319345200000410
n＝0；
获取更新目标模态具体表达式为：
Figure GDA0003531934520000051
其中，i表示第i个模态。
获取更新频率具体表达式为：
Figure GDA0003531934520000052
更新拉格朗日乘子的具体表达式为：
Figure GDA0003531934520000053
其中，τ表示更新参数。
确定响应的设定条件，设定条件具体为：
Figure GDA0003531934520000054
如过满足设定条件，则迭代停止，对获取的优化模态分量作为纯净模态分量进行输出；
如果不满足设定条件，则对优化模态分量继续采用交替方向乘子方法进行交替更新。
表示判别约束条件，且大于0，不满足该条件则返回对获取的目标模态分量和其中心频率进行交替更新；
可以理解的是，通过将Wiener滤波嵌入到VMD算法中，这使得它对采样和噪声的鲁棒性更强；实现根据信号的频率特性完成频带的划分，各模态和中心频率在频域不断循环更新，最终实现信号的自适应分解，以得到所需的纯净模态分量。
在上述任一技术方案中，根据时频图确定有效特征的步骤，具体包括：根据时频图确定纯净模态分量处于的时间范围和纯净模态分量的频率；纯净模态分量处于的时间范围与设定时间范围一致，且纯净模态分量的频率处于设定频率范围内，确定纯净模态分量为有效特征。
在该技术方案中，时频图可以有效反映纯净模态分量的时间和频率的对应关系，根据时频图可以确定纯净模态分量处于的时间范围和纯净模态分量的频率，并对纯净模态分量中时间范围与设定的时间是否一致且频率与设定频率相符的模态分量进行提取，并作为有效特征。
在上述任一技术方案中，设定频率范围为小于50赫兹且大于18赫兹。
在该技术方案中，设定频率范围即为常规冷痛特征的频率，具体为18至50赫兹之间，高于50赫兹可以认为为噪声信号，而低于18赫兹则不是冷痛特征。
可以理解的是，需要在纯净模态分量中提取有用的模态分量，其中设定时间范围为冷刺激时间，由于冷痛特征的模态频率为18到50赫兹，将设定频率范围设置为小于50赫兹且大于18赫兹，实现根据冷刺激的时间和冷痛特征的频率实现在纯净模态分量中筛选出有用的模态分量作为有效特征。
在上述任一技术方案中，从有效信号中筛选出纯净信号的步骤，具体包括：根据含冷痛信息的脑电信号和纯净信号确定有效信号的信噪比；有效信号的信噪比大于等于设定数值确定有效信号为纯净信号。有效信号的信噪比小于设定数值确定有效信号为不合格信号；将不合格信号作为目标信号返回执行将目标信号进行迭代降噪以得到纯净模态分量的步骤。
在该技术方案中，根据纯净信号和初始通过在脑电图中采集到的含冷痛信息的脑电信号计算得到信噪比，并且设置与信噪比相应的设定数值，有效信号的信噪比大于设定数值时，确定纯净信号数量较多则判定有效信号为纯净信号，当有效信号的信噪比小于设定数值时，确定纯净信号数量较少则判定有效信号为不合格信号，并将不合格信号作为目标信号执行将目标信号进行迭代降噪以得到纯净模态分量的步骤。
可以理解的是，将设定数值设置为9，即信噪比大于9时则确定有效信号为纯净信号，当信噪比小于9时确定有效信号为不合格信号。
在上述任一技术方案中，通过以下公式确定有效信号的信噪比：
Figure GDA0003531934520000061
其中，x(ti)为含冷痛信息的脑电信号，
Figure GDA0003531934520000062
为纯净信号。
在该技术方案中，通过含冷痛信息的脑电信号和纯净信号根据上述公式可以计算出信噪比。
在上述任一技术方案中，精度评价参数包括：敏感度、特异性、准确度和阳性预测值。
在该技术方案中，引入四个评价指标敏感度特异性,阳性预测值和准确度，对纯净信号的测试集进行分类，根据四个评价指标可以对纯净信号的精度进行评价。
具体地，对敏感度、特异性、准确度和阳性预测值的计算公式具体为：
Figure GDA0003531934520000071
Figure GDA0003531934520000072
Figure GDA0003531934520000073
Figure GDA0003531934520000074
其中TP表示正类中被正确分类的个数；TN表示负类中被正确分类的个数；FP表示负类中被错分成正类的个数；FN表示正类中被错分成负类的个数。
在上述任一技术方案中，纯净信号的精度评价参数未满足精度需求，将纯净信号作为目标信号返回执行将目标信号进行迭代降噪以得到纯净模态分量的步骤。
在该技术方案中，如果检测到纯净信号的检测精度未满足精度需求，则将未满足精度需求的纯净信号作为目标信号返回对目标信号进行迭代降噪处理的步骤。
在上述任一技术方案中，根据纯净信号获取平均功率谱密度，并将平均功率谱密度组合成特征数据集；按照3比1的比例随机将特征数据集划分为训练数据集和测试数据集；将训练数据集输入到SVM分类器中训练以获取分类模型。
在该技术方案中，将纯净信号平均功率谱密度按3：1比例随机分成训练集和测试集，并使用十折交叉验证策略训练最优分类器，将训练数据集输入到SVM分类器中进行训练得到分类模型，其中径向核函数具体为：
K(x,y)＝exp(-|x-y|2/d2)；
g＝1/k；
k＝10；
分类函数具体为：
Figure GDA0003531934520000081
本发明的第二方面提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行上述任一项的脑电信号伪迹去除方法的步骤，因而具有上述任一项的脑电信号伪迹去除方法的有益技术效果，在此不再进行一一赘述。
本发明的附加方面和优点将在下面的描述部分中变得明显，或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
图1示出了本发明一个实施例中的脑电信号伪迹去除方法的流程示意图；
图2示出了本发明另一个实施例中的脑电信号伪迹去除方法的流程示意图；
图3示出了根据本发明再一个实施例中的脑电信号伪迹去除方法的流程示意图；
图4示出了根据本发明又一个实施例中的脑电信号伪迹去除方法的流程示意图；
图5示出了根据本发明又一个实施例中的脑电信号伪迹去除方法的流程示意图；
图6示出了根据本发明又一个实施例中纯净模态数k＝6时，罚函数α＝10000时的VMD分解和Hilbert变换频谱图；
图7示出了根据本发明又一个实施例中纯净模态数k＝6时，罚函数α＝2000时的VMD分解和Hilbert变换频谱图；
图8示出了根据本发明又一个实施例中纯净模态数k＝4时，罚函数α＝10000时的VMD分解和Hilbert变换频谱图；
图9示出了根据本发明又一个实施例中纯净模态数k＝4时，罚函数α＝10000时的短时傅立叶时频图；
图10示出了根据本发明又一个实施例中Fp1-A1通道和Fp2-A2分解前信号时域图；
图11示出了根据本发明又一个实施例中Fp1-A1通道和Fp2-A2分解后，构造有效信号的时域图；
图12示出了根据本发明又一个实施例中SVM二维平面二分类示意图；
图13示出了根据本发明又一个实施例中分类器训练过程中的二维结构图。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互组合。
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不限于下面公开的具体实施例的限制。
下面参照图1至图13描述根据本发明一些实施例的脑电信号伪迹去除方法及计算机可读存储介质。
如图1所示，本发明的实施例提供了一种脑电信号伪迹去除方法，包括：
S102，获取含冷痛信息的脑电信号作为目标信号，将目标信号进行迭代降噪以得到纯净模态分量；
S104，对纯净模态分量进行短时傅立叶时频分析以得到时频图，根据时频图确定有效特征，构造有效特征对应的模态分量以得到有效信号，从有效信号中筛选出纯净信号；
S106，将纯净信号的测试数据集输入到分类模型中，获取纯净信号的精度评价参数；
S108，纯净信号的精度评价参数满足精度需求确定纯净信号为合格信号。
在该实施例中，获取脑电图中带有冷痛信息的多通道脑电信号，将含冷痛信息的脑电信号未做目标信号，通过对目标信号进行VMD分解实现对脑电信号中的去除眼电、心电、血流等伪迹信号迭代滤波，实现对脑电信号的迭代降噪得到纯净模态分量，纯净模态分量中存在多个模态分量，对纯净模态分量中的各个模态分量进行短时傅立叶时频分析可以得到每个模态分量的时频图，时频图可以反映出模态分量的时间和频率的关系，根据时频图中的时间范围和频率范围确定有效特征，有效特征为可以作为有效信号的模态分量，将有效特征构造成有效信号，并对有效信号再次进行筛选得到纯净信号，再将纯净信号投入到分类模型中进行检测，以得到用于判断信号精度的精度评价参数，用户可以对各类的精度评价参数类型设置响应的精度需求，如果精度评价参数达到了设定的精度需求则将纯净信号作为合格信号输出，如果精度评价参数未达到设定的精度需求则将得到的纯净信号作为目标信号返回执行迭代降噪的步骤进行去除伪迹降噪，直至纯净信号满足精度需求，通过对脑电图中带有冷痛信息的多通道脑电信号进行VMD分解和降噪，克服了模态混叠、端点效应的技术问题，适合由外部刺激引起的脑神经活动检测，以及通过外部刺激诱发的非植入式头皮EEG(脑电图)的特征提取和降噪。提高了分解降噪的精度同时也提高了提取有效特征的精度，对诱发类脑电活动的临床检测和应用具有很好的技术价值。
如图2所示，在上述任一实施例中，将目标信号进行迭代降噪以得到纯净模态分量的步骤，具体包括：
S202，冷痛信息的脑电信号进行分解确定目标模态分量；
S204，确定目标模态分量的数量和目标模态分量的中心频率；
S206，根据目标模态分量的数量和中心频率对目标模态分量进行优化得到优化模态分量；
S208，判断优化模态分量是否达到设定条件，判断结果为是则执行S210，判断结果为否则将优化模态分量作为目标模态分量并执行S204；
S210，确定优化模态分量为纯净模态分量。
在该实施例中，对多通道的脑电信号进行VMD分解以得到多个目标模态分量，确定目标模态分量的数量和中心频率并根据模态分量的数量和中心频率进行优化得到优化后的优化模态分量，并设定相应的设定条件，判断优化后的优化模态分量是否满足设定条件，如果满足设定条件则将优化模态分量作为纯净模态分量输出，如果未满足设定条件则将优化模态分量作为目标模态分量再次进行确定目标模态分量的数量和中心频率的步骤，直至得到的优化模态分量满足设定条件。
具体地，对每个目标模态分量利用Hilbert变换计算目标模态分量的频谱，具体公式如下：
Figure GDA0003531934520000111
Figure GDA0003531934520000112
其中，{uk}：＝{u1，…，uK}为所有目标模态分量，{ωk}：＝{ω1，…，ωK}为所有目标模态分量的中心频率集，uk代表第k个模态，ωk代表第k个模态的中心频率，
Figure GDA0003531934520000113
为所有目标模态分量的和。
将得到的频谱利用二次惩罚因子α和拉格朗日乘子λ(t)获取目标模态分量和其中心频率，从而保证信号重构之后的保真度以及降低噪声的效果。
其中增广拉格朗日表达式为：
Figure GDA0003531934520000114
对获取的目标模态分量和其中心频率通过交替方向乘子方法(ADMM)，进行交替更新，通过计算
Figure GDA0003531934520000115
和
Figure GDA0003531934520000116
来寻求增广拉格朗日表达式的‘鞍点’，‘鞍点’为拉格朗日表达式中既不是最大值也不是最小值的临界点，具体表达式为：
Figure GDA0003531934520000117
利用L2范数下的Parseval/Plancherel傅立叶等距法，将步骤上述表达式转换成谱域内求解，使每个目标模态分量都能获取最新的状态，转换后的表达式具体为：
Figure GDA0003531934520000118
将目标模态分量ω用ω-ωk替换第一项，根据Hermitain对称性质，变换为非负频率区间内的积分形式，具体表达式为：
Figure GDA0003531934520000121
所有的目标模态分量可以在频域中的解中获取，第k个目标模态分量更新的表达式为：
Figure GDA0003531934520000122
对目标模态分量个数k值与二次罚函数α进行定义，并初始化
Figure GDA0003531934520000123
Figure GDA0003531934520000124
n＝0；
获取更新目标模态具体表达式为：
Figure GDA0003531934520000125
其中，i表示第i个模态。
获取更新频率具体表达式为：
Figure GDA0003531934520000126
更新拉格朗日乘子的具体表达式为：
Figure GDA0003531934520000127
其中，τ表示更新参数。
确定响应的设定条件，设定条件具体为：
Figure GDA0003531934520000128
如过满足设定条件，则迭代停止，对获取的优化模态分量作为纯净模态分量进行输出；
如果不满足设定条件，则对优化模态分量继续采用交替方向乘子方法进行交替更新。
表示判别约束条件，且大于0，不满足该条件则返回对获取的目标模态分量和其中心频率进行交替更新；
可以理解的是，通过将Wiener滤波嵌入到VMD算法中，这使得它对采样和噪声的鲁棒性更强；实现根据信号的频率特性完成频带的划分，各模态和中心频率在频域不断循环更新，最终实现信号的自适应分解，以得到所需的纯净模态分量。
如图3所示，在上述任一实施例中，根据时频图确定有效特征的步骤，具体包括：
S302，根据时频图确定纯净模态分量处于的时间范围和纯净模态分量的频率；
S304，根据纯净模态分量处于的时间范围与设定时间范围一致，且纯净模态分量的频率处于设定频率范围内，确定纯净模态分量为有效特征。
具体地，设定频率范围为小于50赫兹且大于18赫兹。
在该技术方案中，时频图可以有效反映纯净模态分量的时间和频率的对应关系，根据时频图可以确定纯净模态分量处于的时间范围和纯净模态分量的频率，并对纯净模态分量中时间范围与设定的时间是否一致且频率与设定频率相符的模态分量进行提取，并作为有效特征。
可以理解的是，需要在纯净模态分量中提取有用的模态分量，其中设定时间范围为冷刺激时间，由于冷痛特征的模态频率为18到50赫兹，将设定频率范围设置为小于50赫兹且大于18赫兹，实现根据冷刺激的时间和冷痛特征的频率实现在纯净模态分量中筛选出有用的模态分量作为有效特征。其中，设定频率范围即为常规冷痛特征的频率，具体为18至50赫兹之间，高于50赫兹可以认为为噪声信号，而低于18赫兹则不是冷痛特征。
如图4所示，本发明提供的一个实施例中，提出了一种脑电信号伪迹去除方法，包括：
S402，获取目标信号并将目标信号进行迭代降噪以得到纯净模态分量；
S404，对纯净模态分量进行短时傅立叶时频分析以得到时频图，根据时频图确定有效特征，构造有效特征对应的模态分量以得到有效信号；
S406，根据含冷痛信息的脑电信号和纯净信号确定有效信号的信噪比；
S408，判断有效信号的信噪比是否大于等于设定数值，判断结果为是则执行S410，判断结果为否则将有效信号作为目标信号返回执行S402；
S410，确定有效信号为纯净信号；
S412，将纯净信号的测试数据集输入到分类模型中，获取纯净信号的精度评价参数；
S414，纯净信号的精度评价参数满足精度需求确定纯净信号为合格信号。
具体地，将设定数值设置为9，即信噪比大于9时则确定有效信号为纯净信号，当信噪比小于9时确定有效信号为不合格信号。
在该实施例中，根据纯净信号和初始通过在脑电图中采集到的含冷痛信息的脑电信号计算得到信噪比，并且设置与信噪比相应的设定数值，有效信号的信噪比大于设定数值时，确定纯净信号数量较多则判定有效信号为纯净信号，当有效信号的信噪比小于设定数值时，确定纯净信号数量较少则判定有效信号为不合格信号，并将不合格信号作为目标信号执行将目标信号进行迭代降噪以得到纯净模态分量的步骤。
在上述任一实施例中，通过以下公式确定有效信号的信噪比：
Figure GDA0003531934520000145
其中，x(ti)为含冷痛信息的脑电信号，
Figure GDA0003531934520000146
为纯净信号。
在该实施例中，通过含冷痛信息的脑电信号和纯净信号根据上述公式可以计算出信噪比。
在上述任一实施例中，精度评价参数包括：敏感度、特异性、准确度和阳性预测值。
在该实施例中，引入四个评价指标敏感度特异性,阳性预测值和准确度，对纯净信号的测试集进行分类，根据四个评价指标可以对纯净信号的精度进行评价。
具体地，对敏感度、特异性、准确度和阳性预测值的计算公式具体为：
Figure GDA0003531934520000141
Figure GDA0003531934520000142
Figure GDA0003531934520000143
Figure GDA0003531934520000144
其中TP表示正类中被正确分类的个数；TN表示负类中被正确分类的个数；FP表示负类中被错分成正类的个数；FN(表示正类中被错分成负类的个数。
在上述任一实施例中，根据纯净信号获取平均功率谱密度，并将平均功率谱密度组合成特征数据集；按照3比1的比例随机将特征数据集划分为训练数据集和测试数据集；将训练数据集输入到SVM分类器中训练以获取分类模型。
在该技术方案中将纯净信号平均功率谱密度按3：1比例随机分成训练集和测试集，并使用十折交叉验证策略训练最优分类器，将训练数据集输入到SVM分类器中进行训练得到分类模型，其中径向核函数具体为：
K(x,y)＝exp(-|x-y|2/d2)；
g＝1/k；
k＝10；
分类函数具体为：
Figure GDA0003531934520000151
如图5所示，本发明提供的一个实施例中，提出了一种脑电信号伪迹去除方法，包括：
S502，获取含冷痛信息的脑电信号作为目标信号；
S504，定义目标模态分量的个数k和罚函数值α；
S506，初始化目标模态分量和目标模态分量的中心频率集；
S508，更新拉格朗日乘子、目标模态分量和目标模态分量的中心频率集；
S510，判断是否
Figure GDA0003531934520000152
判断结果为是则执行S512，判断结果为否则执行S504；
S512，停止迭代获取纯净模态分量；
S514，对第k个纯净模态分量进行短时傅立叶时频分析；
S516，判断纯净模态频率是否在0-50HZ，判断结果为是则执行S518，判断结果为否则返回执行S514；
S518，确定该纯净模态为有效特征；
S520，对有效特征的模态进行重构得到有效信号；
S522，判断有效信号的信噪比是否大于9dB，判断结果为是则执行S524，判断结果为否则执行S504；
S524，确定该信号为纯净信号；
S526，将纯净信号带入带有核函数的SVM分类器中训练；
S528，生成分类模型；
S530，收集纯净信号的测试数据集；
S532，评价纯净信号的检测精度。
在该实施例中，基于VMD(变分模态分解，Variational Mode Decomposition)技术，对含有冷痛诱发的脑活动的多通道信号进行检测和特征提取。选择合适的模态分量的个数k和罚函数α，对含冷痛信息的脑电信号中的单通道EEG信号分解得到的督工额纯净模态分量分别进行短时傅立叶时频分析，以得到时频图，根据时频图获取纯净模态分量的频率以及所在时间范围，判断纯净模态分量的频率以及所在时间范围是否处于设定的时间范围和频率以确定带有有效特征的模态分量，其中频率范围为冷痛特征的模态频率18-50赫兹，时间范围也为收到冷刺激的时间。根据带有有效特征的模态分量进行构造有效信号，对有效信号的信噪比进行计算，并根据有效信号的信噪比是否大于9dB进行筛选，以得到纯净信号。将有效特征输入到SVM分类器中训练以获取分类模型，对纯净信号选取测试数据集，将测试数据集输入到分类模型，获取得到精度评价参数，并根据精度评价参数评价纯净信号是否合格。VMD不仅可以分解出有效特征，也可以实现信号的信噪分离。
在上述实施例中，如图10和图11所示，针对具体例子进行阐述。15个身体健康的被试，在120内接受每次间隔20s的冷痛刺激，其间共两次刺激，采样频率fs＝1000Hz。试验在明亮、安静、温度控制在24-26℃的隔离空间内进行，采集电极为Fp1,Fp2,F7,F3,Fz,F4,F8,T3,C3,Cz,C4,T4,A1,T5,P3,Pz,P4,T6,A2,O1,O2，参考电极是A1和A2，一个接地电极。在测试期间，参与测试者始终平躺在木床上保持静息，不做任何动作和思考，冷痛刺激的介质是0℃的冰水。全程刺激由3名实验人员辅助，第一人观察时间发布指令；第二人根据上一人指令“实验开始”和“实验结束”负责操作EEG采集设备记录信号；第三人负责根据第一人指令“开始刺激”和“离开刺激”将含有冰水的宽口桶完全浸没受试的左手和拿离冰桶。第一次刺激在第20s，持续时间30s；第二次刺激在第70s，持续时间为30s，到120s时候，整组测试结束。每个被试累计做6组测试，每做完一组测试休息30分钟，再进行第二组测试，以求数据的可靠性。脑电图中的多通道脑电信号采集通道采用国际通用导联10-20国际电极安放标准，采集的脑电图中的多通道脑电信号将逐步按照上述实施例中的脑电信号伪迹去除方法进行分解降噪和特征提取。并将若干获取得到的平均功率谱密度特征，按3：1的比例随机分成训练数据集和测试数据集输入到SVM分类器中，对使用十折交叉验证策略对分类器进行训练以得到最佳的分类模型。将测试数据集输入到分类模型中，采用“一对一”的二分类模式，所获得的19个通道经VMD分解的精度如表1所示。
表1
Electrode locations(No.) ACC(％) SEN(％) SPE(％) PPV(％)
Fp1-A1(2) 96.9 99.8 98.7 97.5
Fp2-A2(4) 98.0 97.8 96.3 97.9
F3-A1(6) 99.1 97.1 96.9 98.7
F4-A2(8) 97.3 98.6 98.1 95.2
C3-A1(10) 96.4 96.1 99.4 96.4
C4-A2(12) 97.5 96.6 94.8 97.7
P3-A1(14) 95.2 99.4 96.6 98.1
P4-A2(15) 98.9 96.4 97.3 98.0
O1-A1(13) 98.2 98.6 98.9 97.3
O2-A2(11) 98.8 96.7 96.8 97.9
F7-A1(9) 97.5 96.4 96.2 97.1
F8-A2(7) 96.9 96.8 99.1 98.3
T3-A1(5) 94.9 97.7 95.8 96.5
T4-A2(3) 98.2 95.9 96.7 98.3
续表
Figure GDA0003531934520000171
Figure GDA0003531934520000181
如图12和图13所示，为了验证纯净模态的个数k和罚函数α的变化对纯净信号的提取与检测精度。本发明根据不同纯净模态的个数k和罚函数α，将测试数据集导入分类模型中进行分类，获取的冷痛特征检测平均精度如表2所示。
表2
Figure GDA0003531934520000182
如图6至图9所示，并请参阅以及实验数据表2，增加罚函数α不会改变中心频率,只会缩小冷痛觉的频带，会增加精确度。过大的α反而又会降低精确度。而增加k降低冷痛觉中心频率，模态超过4之后，精确度开始降低。另外，VMD分解过程中，直接通过改变k和α，会缩小特征向量的中心频率和频率范围，能够获取更精确冲击频率范围。但过高的α也会混有更多的噪声，从而影响分析结果。
可以理解的是，通过对脑电图中带有冷痛信息的多通道脑电信号进行VMD分解和降噪，克服了模态混叠、端点效应，使时频分辨率和特征提取精度提高。并且能够通过改变模态k值和二次罚函数α来缩小中心频率范围，对冷痛诱发脑活动特征提取精度增高。该脑电信号去除伪迹方法适合由外部刺激引起的脑神经活动检测，以及通过外部刺激诱发的非植入式头皮EEG的特征提取和降噪，对诱发类脑电活动的临床检测和应用具有很好的技术价值。
可以理解的是，选择合适的纯净模态个数k和罚函数α既可以提高VMD分解的精度，也可以提高有效特征提取的精度。
具体地，精度评价参数包括敏感度、特异性、准确度和阳性预测值。
可以理解的是，如果纯净信号不合格，将纯净信号作为目标信号返回执行定义目标模态分量的个数k和罚函数α值的步骤，实现将纯净型号作为将目标信号进行迭代降噪以得到纯净模态分量的效果。
本发明的第二方面实施例中提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行上述任一实施例中脑电信号伪迹去除方法的步骤，因而具有上述任一实施例中的脑电信号伪迹去除方法的有益技术效果，在此不再进行一一赘述。
在本发明中，术语“多个”则指两个或两个以上，除非另有明确的限定。术语“安装”、“相连”、“连接”、“固定”等术语均应做广义理解，例如，“连接”可以是固定连接，也可以是可拆卸连接，或一体地连接；“相连”可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
在本说明书的描述中，术语“一个实施例”、“一些实施例”、“具体实施例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或实例。而且，描述的具体特征、结构、材料或特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种脑电信号伪迹去除方法，其特征在于，包括：
获取含冷痛信息的脑电信号作为目标信号，将所述目标信号进行迭代降噪以得到纯净模态分量；
对所述纯净模态分量进行短时傅立叶时频分析以得到时频图，根据所述时频图确定有效特征，构造所述有效特征对应的模态分量以得到有效信号，从所述有效信号中筛选出纯净信号；
将所述纯净信号的测试数据集输入到分类模型中，获取所述纯净信号的精度评价参数；
所述纯净信号的精度评价参数满足精度需求确定所述纯净信号为合格信号；
所述根据所述时频图确定有效特征的步骤，具体包括：
根据所述时频图确定所述纯净模态分量处于的时间范围和所述纯净模态分量的频率；
所述纯净模态分量处于的时间范围与设定时间范围一致，且所述纯净模态分量的频率处于设定频率范围内，确定所述纯净模态分量为有效特征；
所述从所述有效信号中筛选出纯净信号的步骤，具体包括：
根据所述含冷痛信息的脑电信号和所述纯净信号确定有效信号的信噪比；
所述有效信号的信噪比大于等于设定数值确定所述有效信号为纯净信号；
所述有效信号的信噪比小于设定数值确定所述有效信号为不合格信号；
将所述不合格信号作为目标信号返回执行所述将所述目标信号进行迭代降噪以得到纯净模态分量的步骤。
2.根据权利要求1所述的脑电信号伪迹去除方法，其特征在于，所述将所述目标信号进行迭代降噪以得到纯净模态分量的步骤，具体包括：
对所述冷痛信息的脑电信号进行分解确定目标模态分量；
确定所述目标模态分量的数量和所述目标模态分量的中心频率；
根据所述目标模态分量的数量和中心频率对所述目标模态分量进行优化得到优化模态分量；
所述优化模态分量达到设定条件，确定所述优化模态分量为纯净模态分量；
所述优化模态分量未达到设定条件，将所述优化模态分量作为所述目标模态分量返回执行确定所述目标模态分量的数量和所述目标模态分量的中心频率的步骤。
3.根据权利要求2所述的脑电信号伪迹去除方法，其特征在于，
所述设定频率范围为小于50赫兹且大于18赫兹。
4.根据权利要求3所述的脑电信号伪迹去除方法，其特征在于，通过以下公式确定所述有效信号的信噪比：
Figure FDA0003531934510000021
其中，x(ti)为含冷痛信息的脑电信号，
Figure FDA0003531934510000022
为纯净信号。
5.根据权利要求1所述的脑电信号伪迹去除方法，其特征在于，所述精度评价参数包括：
敏感度、特异性、准确度和阳性预测值。
6.根据权利要求5所述的脑电信号伪迹去除方法，其特征在于，
所述纯净信号的精度评价参数未满足精度需求，将所述纯净信号作为目标信号返回执行所述将所述目标信号进行迭代降噪以得到纯净模态分量的步骤。
7.根据权利要求6所述的脑电信号伪迹去除方法，其特征在于，还包括：
根据所述纯净信号获取平均功率谱密度，并将所述平均功率谱密度组合成特征数据集；
按照3比1的比例随机将所述特征数据集划分为训练数据集和所述测试数据集；
将所述训练数据集输入到SVM分类器中训练以获取所述分类模型。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述的脑电信号伪迹去除方法的步骤。
CN201910893087.2A 2019-09-20 2019-09-20 脑电信号伪迹去除方法及计算机可读存储介质 Active CN110575164B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910893087.2A CN110575164B (zh) 2019-09-20 2019-09-20 脑电信号伪迹去除方法及计算机可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910893087.2A CN110575164B (zh) 2019-09-20 2019-09-20 脑电信号伪迹去除方法及计算机可读存储介质

Publications (2)

Publication Number Publication Date
CN110575164A CN110575164A (zh) 2019-12-17
CN110575164B true CN110575164B (zh) 2022-04-12

Family

ID=68813200

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910893087.2A Active CN110575164B (zh) 2019-09-20 2019-09-20 脑电信号伪迹去除方法及计算机可读存储介质

Country Status (1)

Country Link
CN (1) CN110575164B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111543946A (zh) * 2020-05-08 2020-08-18 南京邮电大学 基于改进变分模态分解算法的癫痫脑电信号自动检测方法
CN111956217B (zh) * 2020-07-15 2022-06-24 山东师范大学 一种面向实时脑电信号的眨眼伪迹识别方法及系统
CN112826510B (zh) * 2021-03-11 2023-02-07 清华大学 电生理信号采集系统
CN114469135A (zh) * 2022-01-17 2022-05-13 北京脑陆科技有限公司 多导联脑电信号的伪迹去除方法、装置及脑机接口
CN114869294B (zh) * 2022-05-05 2023-05-30 电子科技大学 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法
CN115500829A (zh) * 2022-11-24 2022-12-23 广东美赛尔细胞生物科技有限公司 一种应用于神经内科的抑郁症检测分析系统
CN116058803B (zh) * 2023-03-10 2024-04-19 南开大学 一种面向多类型噪声伪迹的eeg信号去噪方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101869477B (zh) * 2010-05-14 2011-09-14 北京工业大学 一种自适应脑电信号中眼电伪迹的自动去除方法
CN103610461B (zh) * 2013-10-17 2016-03-02 杭州电子科技大学 基于双密度小波邻域相关阈值处理的脑电信号消噪方法
CN105677035B (zh) * 2016-01-26 2018-08-17 杭州电子科技大学 基于eemd和小波阈值的运动想象脑电信号消噪方法
CN105956623A (zh) * 2016-05-04 2016-09-21 太原理工大学 一种基于模糊熵的癫痫脑电信号分类方法
TWI678188B (zh) * 2017-05-22 2019-12-01 艾德腦科技股份有限公司 腦波檢測裝置
CN109633271A (zh) * 2019-01-17 2019-04-16 长沙理工大学 基于变分模态分解和维格纳威尔分布的行波时频分析方法

Also Published As

Publication number Publication date
CN110575164A (zh) 2019-12-17

Similar Documents

Publication Publication Date Title
CN110575164B (zh) 脑电信号伪迹去除方法及计算机可读存储介质
Riaz et al. EMD-based temporal and spectral features for the classification of EEG signals using supervised learning
Gevins Analysis of the electromagnetic signals of the human brain: milestones, obstacles, and goals
Boroujeni et al. Diagnosis of attention deficit hyperactivity disorder using non‐linear analysis of the EEG signal
Sharanreddy et al. EEG signal classification for epilepsy seizure detection using improved approximate entropy
Bhardwaj et al. A novel genetic programming approach for epileptic seizure detection
CN110338786B (zh) 一种癫痫样放电的识别与分类方法、系统、装置和介质
CN111417342B (zh) 脑网络的活动估计系统、方法和存储介质
CN110942103A (zh) 分类器的训练方法和计算机可读存储介质
CN108784693B (zh) 基于独立成分分析和卡尔曼平滑的p300单次提取技术
Li et al. An EEG-based method for detecting drowsy driving state
CN106236080B (zh) 基于多通道的脑电信号中肌电噪声的消除方法
CN110658851B (zh) 一种基于脑电信号的无人机航迹规划系统
CN105342605A (zh) 一种去除脑电信号中肌电伪迹的方法
CN114533086A (zh) 一种基于空域特征时频变换的运动想象脑电解码方法
CN105266804B (zh) 一种基于低秩和稀疏矩阵分解的脑电信号处理方法
CN109009098B (zh) 一种运动想象状态下的脑电信号特征识别方法
Ren et al. MUAP extraction and classification based on wavelet transform and ICA for EMG decomposition
CN110338760B (zh) 一种基于脑电图频域数据的精神分裂症三分类方法
Moser et al. Classification and detection of single evoked brain potentials using time-frequency amplitude features
CN116898454B (zh) 基于脑电特征融合深度学习模型的癫痫分类方法及系统
CN111887811B (zh) 基于脑电信号特征的大脑异常放电检测方法及系统
CN106859641A (zh) 一种基于自动ica去除eeg信号中核磁伪迹的方法
Shete et al. Comparison of sub-band decomposition and reconstruction of EEG signal by daubechies9 and symlet9 wavelet
Yuwono et al. Feature selection with Lasso for classification of ischemic strokes based on EEG signals

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20191217

Assignee: Guilin yanchuang Semiconductor Technology Co.,Ltd.

Assignor: GUILIN University OF ELECTRONIC TECHNOLOGY

Contract record no.: X2023980046590

Denomination of invention: Method for Removing Artifacts in EEG Signals and Computer readable Storage Media

Granted publication date: 20220412

License type: Common License

Record date: 20231108

EE01 Entry into force of recordation of patent licensing contract