CN114098763A - 一种脑电去噪方法 - Google Patents

一种脑电去噪方法 Download PDF

Info

Publication number
CN114098763A
CN114098763A CN202111522916.XA CN202111522916A CN114098763A CN 114098763 A CN114098763 A CN 114098763A CN 202111522916 A CN202111522916 A CN 202111522916A CN 114098763 A CN114098763 A CN 114098763A
Authority
CN
China
Prior art keywords
electroencephalogram
signal
components
main
sub
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.)
Granted
Application number
CN202111522916.XA
Other languages
English (en)
Other versions
CN114098763B (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.)
Shenzhen International Graduate School of Tsinghua University
Original Assignee
Shenzhen International Graduate School of Tsinghua University
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 Shenzhen International Graduate School of Tsinghua University filed Critical Shenzhen International Graduate School of Tsinghua University
Priority to CN202111522916.XA priority Critical patent/CN114098763B/zh
Publication of CN114098763A publication Critical patent/CN114098763A/zh
Application granted granted Critical
Publication of CN114098763B publication Critical patent/CN114098763B/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/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]
    • A61B5/372—Analysis of electroencephalograms
    • 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

Landscapes

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

Abstract

一种脑电去噪方法，包括如下步骤：S1、通过奇异谱分析将原始脑电信号分解为不同的子成分，并设定自相关系数阈值来提取脑电信号主体成分相关的子成分，由子成分重构得到对应的时间序列构成所述脑电信号主体成分；S2、从所述原始脑电信号中减去所述脑电信号主体成分得到残留脑电成分；S3、将所述残留脑电成分输入深度卷积神经网络提取细节特征，再将所述深度卷积神经网络的输出结合所述脑电信号主体成分以重建除去噪后的脑电信号。本发明的脑电去噪方法可以很好地将脑电信号中的肌电噪声去除，对比传统去噪方法在信噪比大于‑1dB时具有更好的去噪性能，相对单一路径具有更小的失真。

Description

一种脑电去噪方法
技术领域
本发明涉及脑电信号处理领域，特别是涉及一种脑电去噪方法。
背景技术
脑电图(Electroencephalography,EEG)是一种记录大脑自发性生物电位的精密电子仪器，被广泛应用于脑科学研究、神经疾病的临床诊断以及脑机接口等领域。然而，EEG在采集过程中往往会不可避免地被参杂入各种其它类型的生理电信号，例如眼电、心电以及肌电等。在大多数脑电等应用场景中，这些生理电信号会干扰脑电信号的分析或处理过程，因而被视为噪声[1]。相比于眼电、心电或者情绪扰动给脑电带来的影响，肌电(Electromyography,EMG)扰动很难直接通过形态学特征或者频谱分析进行检测和分离。
传统的脑电去噪方法往往通过经验特征对含有肌电的异常片段进行检测，并将异常片段整段去除，这带来的问题是造成了脑电时序的不连续并且同时丢失了有价值的脑电信息。理想的肌电消除方式是在被污染的脑电信号上消除加性噪声，保留原始的脑电信息。当存在附加电极通道以获取肌电参考信号时，可以利用附加通道的相关信息来对肌电噪声进行估计[2]。当不存在附加参考信号时，盲源分离技术(Blind Source Separation,BSS)成为了主流的去噪方式。独立成分分析(Independent Component Analysis,ICA)作为一种广泛应用的盲源分离技术，能够将混合前统计独立的非高斯信源进行分离[3]。然而，肌电信源的数量和位置的不稳定性导致ICA在含有肌电的信号盲源分离任务中并不能保持很好的效果。基于二阶统计量优化的典型相关分析(Canonical Correlation Analysis,CCA)作为盲源分离步骤被应用于肌电噪声的去除[4]，其效果优于独立成分分析。
在可穿戴设备与大数据相关技术愈发流行的当下，人们愈加强调移动健康设备的便携性与舒适度。盲源分离技术基于多通道信息，对于未来的便携式可穿戴脑电设备而言，基于单通道信号的去噪方式得到越来越多的关注。奇异谱分析(Singular SpectrumAnalysis)从单个通道EEG中将信号拆解成趋势、振荡和噪声等成分，被用于大脑活动信息的分离和噪声的去除[5，6]。然而，奇异谱分析对于脑电的细节特征并不敏感。深度学习被广泛引用于计算机视觉与信号处理领域，张等人通过搭建卷积神经网络来实现端到端的脑电去噪路径[7]。虽然深度网络能够在极高的污染下也能部分获取脑电信息，但是结果显示在较高信噪比(SNR)的情况下脑电形态较真值相比于奇异谱分析会有更大的失真。
参考文献：
[1]S.Sadiya,T.Alhanai,and M.M.Ghassemi,“Artifact detection andcorrection in EEG data:A review,”in Proc.10th International IEEE-EMBSConference on Neural Engineering(NER),Prague,May.2021,pp.495-498.
[2]Y.Li,P.Wang,and C.Liu,“Electromyogram(EMG)removal by addingsources of EMG(ERASE)-a novel ICA-based algorithm for removing myoelectricartifacts from EEG,”Frontiers in Neuroscience,vol.14,Jan.2021.
[3]L.Albera,A.Kachenoura,and P.Comon,“ICA-based EEG denoising:acomparative analysis of fifteen methods,”Bulletin of the Polish Academy ofSciences-Technical Sciences,vol.60,no.3,pp.407-418,Sept.2012.
[4]W.D.Clercq,A.Vergult,B.Vanrumste,W.Van Paesschen,and S.Van Huffel,“Canonical correlation analysis applied to remove muscle artifacts from theelectroencephalogram,”IEEE Transactions on Biomedical Engineering,vol.53,no.12,pp.2583-2587,Nov 2006.
[5]J.Cheng,L.Li,C.Li,and Y.Liu,“Remove diverse artifactssimultaneously from a single-channel EEG based on SSA and ICA:Asemi-simulatedstudy,”IEEE Access,vol.7,pp.60276-60289,May 2019.
[6]Q.Liu,A.Liu,X.Zhang,X.Chen,R.Qian,and X.Chen,“Removal of EMGartifacts from multichannel EEG signals using combined singular spectrumanalysis and canonical correlation analysis,”Journal of HealthcareEngineering,vol.2019,Dec.2019.
[7]Z.Haoming,W.Chen,Z.Mingqi,W.Haiyan,and L.Quanying,“Anovelconvolutional neural network model to remove muscle artifacts from EEG,”inProc.IEEE International Conference on Acoustics,Speech and Signal Processing(ICASSP),Jun.2021,pp.1265-1269.
需要说明，在上述背景技术部分公开的信息仅用于对本申请的背景的理解，因此可以包括不构成对本领域普通技术人员已知的现有技术的信息。
发明内容
本发明的主要目的在于克服上述背景技术的缺陷，提供一种脑电去噪方法。
为实现上述目的，本发明采用以下技术方案：
一种脑电去噪方法，包括如下步骤：
S1、通过奇异谱分析(SSA)将原始脑电信号分解为不同的子成分，并设定自相关系数阈值来提取脑电信号主体成分相关的子成分，由子成分重构得到对应的时间序列构成所述脑电信号主体成分；
S2、从所述原始脑电信号中减去所述脑电信号主体成分得到残留脑电成分；
S3、将所述残留脑电成分输入深度卷积神经网络提取细节特征，再将所述深度卷积神经网络的输出结合所述脑电信号主体成分以重建除去噪后的脑电信号。
进一步地：
步骤S1包括如下步骤：
S11、构建轨迹矩阵：对于一段长度为N的时间序列X，设置窗口长度设置L，构建轨迹矩阵，其中窗口长度L满足需要分解到的感兴趣最小频率成分fmin；
S12、奇异值分解：对所述轨迹矩阵进行奇异值分解，得到轨迹矩阵子成分；
S13、重构分组：对于每个轨迹矩阵子成分Y1,Y2,…,YL进行反对角线重构得到对应的时间序列x1,x2,…,xL，计算重构后时间序列x1,x2,…,xL对应的自相关系数为r1,r2,…,rL，利用自相关系数阈值Rthr进行分组，得到的时间序列组合构成所述脑电信号主体成分。
步骤S11包括：
对于长度为N的时间序列X，窗口长度为L，所述轨迹矩阵如下：
Figure BDA0003408433250000031
其中，参数K＝N-L+1；对于需要分解到的感兴趣最小频率成分fmin，窗口长度满足L≥fs/fmin，fs对应于采样频率；优选地，窗口长度为32。
步骤S12包括：
求所述轨迹矩阵的奇异值分解转换为求对应协方差矩阵的特征值分解；其中，对于C＝YYT，C的特征值即Y矩阵奇异值的平方，按降序排列为λ1,λ2,…,λL，对应的特征向量即Y矩阵左奇异向量为u1,u2,…,uL,令
Figure BDA0003408433250000041
即Y矩阵右奇异向量，第i个轨迹矩阵子成分为：
Figure BDA0003408433250000042
Figure BDA0003408433250000043
并有Y＝Y1+Y2+…+YL；Yi＝uiui TY，uiui T构成了子带信息。
步骤S13包括：
对于每个轨迹矩阵子成分Y1,Y2,…,YL进行反对角线重构得到对应的时间序列x1,x2,…,xL，其自相关系数如下：
Figure BDA0003408433250000044
其中E[]为期望算符，xo(t)表示为时间序列，xd(t)为对应的一个采样点延时序列；计算重构后时间序列x1,x2,…,xL对应的自相关系数为r1,r2,…,rL。
步骤S13包括：
在训练阶段，用于分组的自相关系数阈值Rthr通过计算不同重建后时间序列的组合与脑电信号真值的相对均方根误差来求得：
Figure BDA0003408433250000045
使得对应自相关系数大于Rthr的重建后时间序列的组合与脑电信号真值Xpure具有最小的均方根误差，所述时间序列组合构成脑电信号主体成分：
xmain＝∑(xj|rj＞Rthr)
在测试阶段，所述自相关系数阈值被设定为一个经验值。
步骤S2中，从被污染脑电信号中减去脑电信号主体成分得到残留脑电成分：
xres＝∑(xj|rj＜Rthr)。
在卷积神经网络的训练阶段，将脑电信号真值Xpure与脑电信号主体成分xmain的差值即残差干净脑电信号xclean作为有监督训练的标签值：
xclean＝Xpure-xmain
训练得到卷积神经网络模型的参数集
Figure BDA0003408433250000046
得到去噪后的脑电细节信号xnl：
Figure BDA0003408433250000051
步骤S3包括：
使用端到端的卷积神经网络对所述残留脑电成分进行细节特征提取和去噪；其中，将原始脑电信号X与脑电信号主体成分xmain的差值即所述残留脑电成分xres作为网络的输入：
xres＝X-xmain
将所述脑电细节信号xnl结合所述脑电信号主体成分xmain重建得到去噪后的脑电信号Xfre：
Xfre＝xnl+xmain。
所述方法还包括：构建自适应自相关系数阈值模块；其中，将原始脑电信号经过奇异谱分析分解后得到的子带信号按照其自相关系数从高至低的顺序排列，从高至低依次计算将对应子带信号自相关系数作为总体自相关系数阈值后的脑电信号主体成分与脑电真值的相对均方根误差，以获得最优的子带信号自相关系数，作为自相关系数阈值的标签值；由此构建足量的原始脑电信号与对应标记的自相关系数阈值，用于自相关系数阈值模块的训练和测试。
优选地，自适应阈值模块包含特征提取和利用神经网络回归得到自相关系数阈值；
特征提取阶段选择对脑电信号和其它噪声具有区分能力的特征；对于一段长度为N的时间序列X，设定参数m得到(N-m+1)个m维迟滞序列Xm(i)＝{xi,xi+1,xi+2,…,xi+m-1}，其中1i≤N-m+1；两序列间的距离为d[Xm(i),Xm(j)]＝maxk＝0,,…,m-1(|x(i+k)-x(j+k)|；根据相似容限参数r来定义两个序列匹配m个点的概率
Figure BDA0003408433250000052
其中Bi为给定Xm(i)的情况下Xm(j)与之距离小于等于r的数量；两个序列匹配m+1个点的概率为
Figure BDA0003408433250000053
其中Ai为给定Xm+1(i)的情况下Xm+1(j)与之距离小于等于r的数量；序列X样本熵的估计值为
Figure BDA0003408433250000054
标准差SD用于评估序列的离散程度，
Figure BDA0003408433250000055
子带的信号的自相关系数用于以评价肌电成分含量；复杂性被定义为原信号的单位延迟信号迁移性与原信号迁移性的比值；
提取的特征按顺序排列作为神经网络的输入，最优子带自相关系数作为网络学习阈值的标签值；隐藏层神经元个数应大于输入层以学习深层特征；根据网络输出值的范围设定最后一层的激活函数为tanh；将训练后的模型参数用于从特征中回归得到优化后的自相关系数阈值以分离脑电主体成分。
一种计算机可读存储介质，存储有计算机程序，所述计算机程序由处理器执行时，实现所述的脑电去噪方法。
本发明具有如下有益效果：
本发明提出了一种联合奇异谱分析和卷积神经网络的脑电去噪方法。针对脑电信号在采集过程中常常遭受的肌电干扰，本发明提出的联合信号处理方法可以很好地将脑电信号中的肌电噪声去除，相比传统去噪方法在信噪比大于-1dB时具有更好的去噪性能。
本发明将奇异谱分析结合深度卷积神经网络，通过一种融合路径去除脑电信号中肌电成分，先通过奇异谱分析提取脑电的主体节律和形态，再将残差通过卷积神经网络去噪并提取细节特征，这种方法相对单一路径具有更小的失真。
附图说明
图1为本发明实施例的脑电去噪方法流程图。
图2为本发明实施例的SNR＝1dB时被污染的脑电、脑电真值和SSA-CNN联合算法去噪后的脑电的时域(a)和功率谱密度(b)波形图。
图3为本发明实施例的SNR＝1dB时脑电真值、SSA-CNN联合算法去噪后的脑电、SSA-CCA算法去噪后的脑电和CNN去噪后的脑电波形对比图。
图4为不同算法的去噪输出与脑电真值在一定范围信噪比内的RRMSE(a)和CC(b)。
图5为本发明实施例利用特征提取与神经网络计算分解子带信号集合的最优自适应自相关系数的示意图。
具体实施方式
以下对本发明的实施方式做详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。
参阅图1，本发明实施例提供一种脑电去噪方法，包括如下步骤：
S1、通过奇异谱分析将原始脑电信号分解为不同的子成分，并设定自相关系数阈值来提取脑电信号主体成分相关的子成分，由所述子成分重构得到对应的时间序列构成所述脑电信号主体成分；原始脑电信号通常是被污染的脑电信号；
S2、从所述原始脑电信号中减去所述脑电信号主体成分得到残留脑电成分；
S3、将所述残留脑电成分输入深度卷积神经网络提取细节特征，再将所述深度卷积神经网络输出的脑电细节信号结合所述脑电信号主体成分以重建除去噪后的脑电信号。
本发明实施例的脑电去噪方法联合奇异谱分析和卷积神经网络进行处理，通过奇异谱分析将信号分解为不同的子成分，并设定自相关系数阈值来提取脑电主体成分相关的子成分，再将残差部分输入深度网络提取细节特征，最后将深度网络的输出结合脑电主体成分重建除去噪后的脑电信号。具体流程如图1所示。
本发明的联合信号处理方法可以很好地将脑电信号中的肌电噪声去除，对比传统去噪方法在信噪比大于-1dB时具有更好的去噪性能，相对单一路径具有更小的失真。
脑电主体成分提取
奇异谱分析将一段固定长度的时域信号分解为不同的节律成分，这种分析仅依赖数据本身的特征。奇异谱分析主要包括3个步骤：构建轨迹矩阵、奇异值分解和重构分组。
构建轨迹矩阵
对于一段长度为N的时间序列X，当窗口长度设置为L时，轨迹矩阵定义如下：
Figure BDA0003408433250000071
其中，参数K＝N-L+1。当需要分解到的感兴趣最小频率成分为fmin时，窗口长度满足L≥fs/fmin，fs对应于采样频率。经过实验比较，窗口长度可以设置为经验值，如32。
奇异值分解
求轨迹矩阵的奇异值分解可以转换为求对应协方差矩阵的特征值分解。对于C＝YYT，C的特征值(Y矩阵奇异值的平方)按降序排列为λ1,λ2,…,λL，对应的特征向量(Y矩阵左奇异向量)为u1,u2,…,uL,令
Figure BDA0003408433250000072
(Y矩阵右奇异向量)，则第i个轨迹矩阵子成分定义为：
Figure BDA0003408433250000081
并有Y＝Y1+Y2+…+YL。由于Yi＝uiui TY，所以uiui T构成了子带信息。
重构分组
对于每个轨迹矩阵子成分Y1,Y2,…,YL进行反对角线重构得到对应的时间序列x1,x2,…,xL。定义自相关系数如下：
Figure BDA0003408433250000082
其中E[]为期望算符，xo(t)表示为时间序列，xd(t)为对应的一个采样点延时序列。鉴于肌电类似于白噪声，具有较低的自相关系数，所以计算重构后时间序列x1,x2,…,xL对应的自相关系数为r1,r2,…,rL。
在算法的训练阶段，用于分组的自相关系数阈值Rthr可以通过计算不同重建后时间序列的组合与脑电信号真值的相对均方根误差来求得。具体的数学表达为：
Figure BDA0003408433250000083
即使得对应自相关系数大于Rthr的重建后时间序列的组合与脑电信号真值Xpure具有最小的均方根误差，这些时间序列组合构成脑电信号主体成分：
xmain＝∑(xj|rj＞Rthr)
在算法的测试阶段，由于脑电信号真值是未知的，因而自相关系数阈值被设定为一个经验值。
从被污染脑电信号中减去脑电信号主体成分得到残留脑电成分：
xres＝∑(xj|rj＜Rthr)
之后将残留脑电成分送入深度卷积神经网络实现脑电细节成分的提取。
脑电细节成分提取
端到端的卷积神经网络被用于脑电信号的去噪任务当中，实施例使用该网络作为脑电细节成分提取的方法，用于对残留脑电成分的去噪训练。若将原始被污染脑电信号作为网络输入，会造成结果中脑电主体成分的失真。本方案将被污染脑电信号X与脑电主体成分xmain的差值即残留脑电成分xres(残差带噪脑电信号)作为网络的输入：
xres＝X-xmain
在网络的训练过程中，将脑电信号真值Xpure与脑电主体成分xmain的差值xclean(残差干净脑电信号)作为有监督训练的标签值：
xclean＝Xpure-xmain
训练得到卷积神经网络模型的参数集
Figure BDA0003408433250000091
去噪后的脑电细节信号xnl表示为：
Figure BDA0003408433250000092
脑电真值仅用于训练阶段。在训练过程中，被污染脑电信号通过已知的干净脑电加上肌电噪音构建。提取脑电主体成分的最优SSA分解阈值可通过最大化SSA分解(分解是对于被污染信号)后的主体成分与干净脑电的相关系数来确定。xclean是xres的深度网络标签值。
重建除去噪后的脑电信号
将深卷积神经网络输出的脑电细节信号xnl结合脑电信号主体成分以重建除去噪后的脑电信号，即得到去噪后的脑电信号Xfre：
Xfre＝xnl+xmain。
性能分析
为验证算法用于脑电信号中肌电噪声去除的效果，训练和测试阶段使用的数据集构建于现有公开数据集。CHB-MIT数据库包含有23例儿科癫痫病例的脑电数据，其采样率为256Hz，记录时长的中位数达10小时。实施例从CHB-MIT数据集中按照1/16的重叠率切片得到2s时长的脑电片段，通道选择为所有病例都共有的18个电极通道(FP1-F7,F7-T7,T7-P7,P7-O1,FP1-F3,F3-C3,C3-P3,P3-O1,FP2-F4,F4-C4,C4-P4,P4-O2,FP2-F8,F8-T8,T8-P8,P8-O2,FZ-CZ,CZ-PZ)。
由于该数据集在采集的同时存在噪声污染，因而需要设定一些准则用于筛选干净的脑电片段。峰度(Kurtosis)被用于筛除眼电污染片段：
Kurtosis＝E[(x-E[x])4]-3(E[(x-E[x])2])2
较高的Kurtosis意味着脑电片段更大的可能受到眼电污染。迁移性(Mobility,M)被用于筛除肌电污染片段：
Figure BDA0003408433250000101
其中di＝xi-xi-1。较低的Mobility意味着脑电片段更大的可能受到肌电污染。鉴于这些筛除参数，共得到有64000条脑电片段。所得脑电片段通过1Hz高通滤波器以滤除基线漂移。
EEGdenoiseNet整理了多个生理电信号相关的公开数据集形成了供脑电去噪算法使用的模拟数据库。实施例使用该数据库中5598例的采样率为256Hz的肌电数据集作为模拟加性污染源。
模拟污染脑电信号通过以下公式构建：
X＝Xpure+λXemg
其中Xpure为脑电真值片段，Xemg为肌电片段，参数λ控制脑电被污染程度。为了方便与单一路径方法的效果进行对比，信噪比被定义为：
Figure BDA0003408433250000102
其中RMS表示信号片段的均方根。
通过随机选取方式扩充肌电片段至64000例，使之与脑电真值片段数量一致。选取两个信号片段集合中各60000例样本，并在-5至4内随机生成信噪比参数以构建训练集。同时将剩余4000对样本均等分成验证集与测试集各2000对，并将每对样本按照信噪比-5至4以单位间隔量生成10倍扩增验证集20000例和测试集20000例。
在验证集上通过相对均方根误差和相关系数来评价与比较本算法和相关算法性能。相对均方根误差(RRMSE)被定义为：
Figure BDA0003408433250000103
相关系数被定义为：
Figure BDA0003408433250000104
其中Cov(X,Y)和Var(X)分别表示协方差与方差。
图2显示了一例在信噪比为1dB时本算法的去噪效果。其中(a)与(b)分别表示时域波形和频域波形。从(a)中可以看到，被污染脑电经过算法处理过后很好地保留了原始的脑电信息，与脑电真值具有很高的相似性。从(b)中可以看到，高频污染成分被很好的消除且不至于完全丧失脑电自身的高频信息。
图3显示了一例信噪比为1dB时不同算法去噪后的信号与脑电真值的比较结果。绿色信号表示经过联合奇异谱分析和典型相关分析(SSA-CCA)算法的去噪输出，可以看到绿色信号较好的还原了脑电的低频主体形状，但是在丢失了大量的细节信息，这一点在局部峰值处更加明显。蓝色信号表示经过卷积神经网络(CNN)的去噪输出，可以看到蓝色信号存在的问题是对于脑电主体形状的还原有较大的失真度，主要表现为峰值出现缺口。相比之下，联合奇异谱分析和卷积神经网络(SSA-CNN)算法输出与真值差距最小。
图4显示了不同算法在一定范围信噪比内去噪输出与真值之间的RRMSE(a)与CC(b)。可以看出，当SNR大于-1dB时，SSA-CNN算法具有最低的RRMSE，即算法输出相比其它算法与脑电真值具有最小的差距；当SNR大于-3dB时，SSA-CNN的输出与脑电真值的相似度已然最高。因而，在SNR大于-1dB时，可以认为联合奇异谱分析和卷积神经网络的脑电去噪算法具有最优性能。
变型实施例
对于残差带噪脑电和脑电主体成分的分离，依据的自相关系数取值为经验值。尽管依据经验值的联合算法相较于独立路径能够实现肌电去噪效果的提升，但是在实际情况中，不同时间不同位置的脑电受到肌电的污染程度不同，并且肌电信号的自相关系数也非定值，这就造成选择经验值作为阈值会使得分离的效果受到影响，即分离后的脑电主体成分较真值有较大的差异，进而影响后续的脑电细节特征提取。针对此问题，本发明优选实施例将自适应自相关系数阈值模块作为上述技术方案的补充扩展。
构建的自适应自相关系数阈值模块如图5中所示。被污染脑电经过SSA分解，得到的子带信号按照其自相关系数从高至低的顺序排列。从高至低依次计算将对应子带信号自相关系数作为总体自相关系数阈值后的脑电主体成分与脑电真值的相对均方根误差，以获得最优的子带信号自相关系数，作为自相关系数阈值的标签值。根据上述路径构建足量的被污染脑电信号与对应标记的自相关系数阈值，用于自适应模块的训练和测试。
自适应阈值模块包含两个步骤，第一步为特征提取，第二步为利用神经网络回归得到自相关系数阈值。
特征提取阶段需要选择对脑电信号和其它噪声具有区分能力的特征。对于生理信号的规律性，较高的样本熵(Sample Entropy,SE)意味着时间序列具有较高的复杂度，反之亦然。肌电信号通常情况下相较脑电信号有更高的复杂度。对于一段长度为N的时间序列X，设定参数m得到(N-m+1)个m维迟滞序列Xm(i)＝{xi,xi+1,xi+2,…,xi+m-1}，其中1i≤N-m+1。定义两序列间的距离为d[Xm(i),Xm(j)]＝maxk＝0,1,…,m-1(|x(i+k)-x(j+k)|。根据相似容限参数r来定义两个序列匹配m个点的概率
Figure BDA0003408433250000121
其中Bi为给定Xm(i)的情况下Xm(j)与之距离小于等于r的数量。同理定义两个序列匹配m+1个点的概率
Figure BDA0003408433250000122
Figure BDA0003408433250000123
其中Ai为给定Xm+1(i)的情况下Xm+1(j)与之距离小于等于r的数量。此时序列X样本熵的估计值为
Figure BDA0003408433250000124
标准差(Standard Deviation,SD)用于评估序列的离散程度，噪声污染程度越大，信号的标准差往往越大，
Figure BDA0003408433250000125
子带的信号的自相关系数(AutocorrelationCoefficient,R)和迁移性(Mobility,M)同时被选为特征以评价肌电成分含量，较低的自相关系数和较低的迁移性意味着脑电中含有更多的肌电。复杂性(Complexity,C)被定义为原信号的单位延迟信号迁移性与原信号迁移性的比值。
提取的特征按顺序排列作为神经网络的输入，最优子带自相关系数作为网络学习阈值的标签值。隐藏层神经元个数应大于输入层以学习深层特征。根据网络输出值范围设定输出层激活函数为tanh。训练后的模型参数可用于从特征中回归得到优化后的自相关系数阈值以更好地分离脑电主体成分。
本发明的背景部分可以包含关于本发明的问题或环境的背景信息，而不一定是描述现有技术。因此，在背景技术部分中包含的内容并不是申请人对现有技术的承认。
以上内容是结合具体/优选的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，其还可以对这些已描述的实施方式做出若干替代或变型，而这些替代或变型方式都应当视为属于本发明的保护范围。在本说明书的描述中，参考术语“一种实施例”、“一些实施例”、“优选实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。尽管已经详细描述了本发明的实施例及其优点，但应当理解，在不脱离专利申请的保护范围的情况下，可以在本文中进行各种改变、替换和变更。

Claims (10)

1.一种脑电去噪方法，其特征在于，包括如下步骤：
S1、通过奇异谱分析将原始脑电信号分解为不同的子成分，并设定自相关系数阈值来提取脑电信号主体成分相关的子成分，由所述子成分重构得到对应的时间序列构成所述脑电信号主体成分；
S2、从所述原始脑电信号中减去所述脑电信号主体成分得到残留脑电成分；
S3、将所述残留脑电成分输入深度卷积神经网络提取细节特征，再将所述深度卷积神经网络的输出结合所述脑电信号主体成分以重建除去噪后的脑电信号。
2.如权利要求1所述的脑电去噪方法，其特征在于，步骤S1包括如下步骤：
S11、构建轨迹矩阵：对于一段长度为N的时间序列X，设置窗口长度设置L，构建轨迹矩阵，其中窗口长度L满足需要分解到的感兴趣最小频率成分fmin；
S12、奇异值分解：对所述轨迹矩阵进行奇异值分解，得到轨迹矩阵子成分；
S13、重构分组：对于每个轨迹矩阵子成分Y1,Y2,…,YL进行反对角线重构得到对应的时间序列x1,x2,…,xL，计算重构后时间序列x1,x2,…,xL对应的自相关系数为r1,r2,…,rL，利用自相关系数阈值Rthr进行分组，得到的时间序列组合构成所述脑电信号主体成分。
3.如权利要求2所述的脑电去噪方法，其特征在于，步骤S11包括：
对于长度为N的时间序列X，窗口长度为L，所述轨迹矩阵如下：
Figure FDA0003408433240000011
其中，参数K＝N-L+1；对于需要分解到的感兴趣最小频率成分fmin，窗口长度满足L≥fs/fmin，fs对应于采样频率；优选地，窗口长度为32。
4.如权利要求2或3所述的脑电去噪方法，其特征在于，步骤S12包括：
求所述轨迹矩阵的奇异值分解转换为求对应协方差矩阵的特征值分解；其中，对于C＝YYT，C的特征值即Y矩阵奇异值的平方，按降序排列为λ1,λ2,…,λL，对应的特征向量即Y矩阵左奇异向量为u1,u2,…,uL，令
Figure FDA0003408433240000021
即Y矩阵右奇异向量，第i个轨迹矩阵子成分为：
Figure FDA0003408433240000022
Figure FDA0003408433240000023
并有Y＝Y1+Y2+…+YL；Yi＝uiui TY，uiui T构成了子带信息。
5.如权利要求2至4任一项所述的脑电去噪方法，其特征在于，步骤S13包括：
对于每个轨迹矩阵子成分Y1,Y2,…,YL进行反对角线重构得到对应的时间序列x1,x2,…,xL，其自相关系数如下：
Figure FDA0003408433240000024
其中E[]为期望算符，xo(t)表示为时间序列，xd(t)为对应的一个采样点延时序列；计算重构后时间序列x1,x2,…,xL对应的自相关系数为r1,r2,…,rL。
6.如权利要求2至5任一项所述的脑电去噪方法，其特征在于，步骤S13包括：
在训练阶段，用于分组的自相关系数阈值Rthr通过计算不同重建后时间序列的组合与脑电信号真值的相对均方根误差来求得：
Figure FDA0003408433240000025
使得对应自相关系数大于Rthr的重建后时间序列的组合与脑电信号真值Xpure具有最小的均方根误差，所述时间序列组合构成脑电信号主体成分：
xmain＝∑(xj|rj＞Rthr)
在测试阶段，所述自相关系数阈值被设定为一个经验值；
步骤S2中，从原始脑电信号中减去脑电信号主体成分得到残留脑电成分：
xres＝∑(xj|rj＜Rthr)。
7.如权利要求1至6任一项所述的脑电去噪方法，其特征在于，在卷积神经网络的训练阶段，将脑电信号真值Xpure与脑电信号主体成分xmain的差值即残差干净脑电信号xclean作为有监督训练的标签值：
xclean＝Xpure-xmain
训练得到卷积神经网络模型的参数集
Figure FDA0003408433240000031
得到去噪后的脑电细节信号xnl：
Figure FDA0003408433240000032
8.如权利要求1至7任一项所述的脑电去噪方法，其特征在于，
步骤S3包括：使用端到端的卷积神经网络对所述残留脑电成分进行细节特征提取和去噪；其中，将原始脑电信号X与脑电信号主体成分xmain的差值即所述残留脑电成分xres作为网络的输入：
xres＝X-xmain
将脑电细节信号xnl结合所述脑电信号主体成分xmain重建得到去噪后的脑电信号Xfre：
Xfre＝xnl+xmain。
9.如权利要求1至8任一项所述的脑电去噪方法，其特征在于，还包括：构建自适应自相关系数阈值模块；其中，将原始脑电信号经过奇异谱分析分解后得到的子带信号按照其自相关系数从高至低的顺序排列，从高至低依次计算将对应子带信号自相关系数作为总体自相关系数阈值后的脑电信号主体成分与脑电真值的相对均方根误差，以获得最优的子带信号自相关系数，作为自相关系数阈值的标签值；由此构建足量的原始脑电信号与对应标记的自相关系数阈值，用于自相关系数阈值模块的训练和测试；
优选地，自适应阈值模块包含特征提取和利用神经网络回归得到自相关系数阈值；
特征提取阶段选择对脑电信号和其它噪声具有区分能力的特征；对于一段长度为N的时间序列X，设定参数m得到(N-m+1)个m维迟滞序列Xm(i)＝{xi,xi+1,xi+2,…,xi+m-1}，其中1≤i≤N-m+1；两序列间的距离为d[Xm(i),Xm(j)]＝maxk＝0,1,…,m-1(|x(i+k)-x(j+k)|)；根据相似容限参数r来定义两个序列匹配m个点的概率
Figure FDA0003408433240000041
其中Bi为给定Xm(i)的情况下Xm(j)与之距离小于等于r的数量；两个序列匹配m+1个点的概率为
Figure FDA0003408433240000042
其中Ai为给定Xm+1(i)的情况下Xm+1(j)与之距离小于等于r的数量；序列X样本熵的估计值为
Figure FDA0003408433240000043
标准差SD用于评估序列的离散程度，
Figure FDA0003408433240000044
子带的信号的自相关系数用于以评价肌电成分含量；复杂性被定义为原信号的单位延迟信号迁移性与原信号迁移性的比值；
提取的特征按顺序排列作为神经网络的输入，最优子带自相关系数作为网络学习阈值的标签值；隐藏层神经元个数应大于输入层以学习深层特征；根据网络输出值的范围设定最后一层的激活函数为tanh；将训练后的模型参数用于从特征中回归得到优化后的自相关系数阈值以分离脑电主体成分。
10.一种计算机可读存储介质，存储有计算机程序，其特征在于，所述计算机程序由处理器执行时，实现如权利要求1至9任一项所述的脑电去噪方法。
CN202111522916.XA 2021-12-13 2021-12-13 一种脑电去噪方法 Active CN114098763B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111522916.XA CN114098763B (zh) 2021-12-13 2021-12-13 一种脑电去噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111522916.XA CN114098763B (zh) 2021-12-13 2021-12-13 一种脑电去噪方法

Publications (2)

Publication Number Publication Date
CN114098763A true CN114098763A (zh) 2022-03-01
CN114098763B CN114098763B (zh) 2023-05-19

Family

ID=80365348

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111522916.XA Active CN114098763B (zh) 2021-12-13 2021-12-13 一种脑电去噪方法

Country Status (1)

Country Link
CN (1) CN114098763B (zh)

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040015894A1 (en) * 1998-12-30 2004-01-22 Lange Daniel H. Method and apparatus for extracting low SNR transient signals from noise
US20140073870A1 (en) * 2012-09-11 2014-03-13 Nellcor Puritan Bennett Llc Methods and systems for determining physiological information based on a peak of a cross-correlation
US20160106331A1 (en) * 2013-04-22 2016-04-21 The Regents Of The University Of California Fractal index analysis of human electroencephalogram signals
CN107688120A (zh) * 2016-11-17 2018-02-13 丹阳华神电器有限公司 基于模糊熵的含噪信号处理方法及迭代奇异谱软阈值去噪方法
CN108309290A (zh) * 2018-02-24 2018-07-24 华南理工大学 单通道脑电信号中肌电伪迹的自动去除方法
CN108836301A (zh) * 2018-05-04 2018-11-20 江苏师范大学 一种基于奇异谱分析和稀疏表示的单次诱发电位提取方法
CN109009092A (zh) * 2018-06-15 2018-12-18 东华大学 一种去除脑电信号噪声伪迹的方法
CN109558863A (zh) * 2019-01-07 2019-04-02 哈尔滨工业大学（深圳） 脑电信号特征表示与提取方法、装置及存储介质
CN110348632A (zh) * 2019-07-11 2019-10-18 广东电网有限责任公司 一种基于奇异谱分析和深度学习的风电功率预测方法
CN110897639A (zh) * 2020-01-02 2020-03-24 清华大学深圳国际研究生院 一种基于深度卷积神经网络的脑电睡眠分期方法
CN112932501A (zh) * 2021-01-25 2021-06-11 上海海事大学 一种基于一维卷积神经网络自动识别失眠方法
CN113229829A (zh) * 2021-04-15 2021-08-10 广东工业大学 一种四元数脑电信号提取方法及系统

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040015894A1 (en) * 1998-12-30 2004-01-22 Lange Daniel H. Method and apparatus for extracting low SNR transient signals from noise
US20140073870A1 (en) * 2012-09-11 2014-03-13 Nellcor Puritan Bennett Llc Methods and systems for determining physiological information based on a peak of a cross-correlation
US20160106331A1 (en) * 2013-04-22 2016-04-21 The Regents Of The University Of California Fractal index analysis of human electroencephalogram signals
CN107688120A (zh) * 2016-11-17 2018-02-13 丹阳华神电器有限公司 基于模糊熵的含噪信号处理方法及迭代奇异谱软阈值去噪方法
CN108309290A (zh) * 2018-02-24 2018-07-24 华南理工大学 单通道脑电信号中肌电伪迹的自动去除方法
CN108836301A (zh) * 2018-05-04 2018-11-20 江苏师范大学 一种基于奇异谱分析和稀疏表示的单次诱发电位提取方法
CN109009092A (zh) * 2018-06-15 2018-12-18 东华大学 一种去除脑电信号噪声伪迹的方法
CN109558863A (zh) * 2019-01-07 2019-04-02 哈尔滨工业大学（深圳） 脑电信号特征表示与提取方法、装置及存储介质
CN110348632A (zh) * 2019-07-11 2019-10-18 广东电网有限责任公司 一种基于奇异谱分析和深度学习的风电功率预测方法
CN110897639A (zh) * 2020-01-02 2020-03-24 清华大学深圳国际研究生院 一种基于深度卷积神经网络的脑电睡眠分期方法
CN112932501A (zh) * 2021-01-25 2021-06-11 上海海事大学 一种基于一维卷积神经网络自动识别失眠方法
CN113229829A (zh) * 2021-04-15 2021-08-10 广东工业大学 一种四元数脑电信号提取方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
CHENG YU, QINGLAN SHAN, ZHIDE LI,ET AL: "A Joint SSA and CNN Algorithm to Remove EMG Artifacts from Single Channel EEG" *

Also Published As

Publication number Publication date
CN114098763B (zh) 2023-05-19

Similar Documents

Publication Publication Date Title
Yang et al. Automatic ocular artifacts removal in EEG using deep learning
Chen et al. Removing muscle artifacts from EEG data: Multichannel or single-channel techniques?
Inuso et al. Wavelet-ICA methodology for efficient artifact removal from Electroencephalographic recordings
Bono et al. Hybrid wavelet and EMD/ICA approach for artifact suppression in pervasive EEG
Molla et al. Artifact suppression from EEG signals using data adaptive time domain filtering
Shukla et al. Wavelet based empirical approach to mitigate the effect of motion artifacts from EEG signal
Romero PCA-based noise reduction in ambulatory ECGs
Li et al. ECG denoising method based on an improved VMD algorithm
Abbaspour et al. Evaluation of wavelet based methods in removing motion artifact from ECG signal
Patil et al. Quality advancement of EEG by wavelet denoising for biomedical analysis
Dora et al. Adaptive single-channel EEG artifact removal with applications to clinical monitoring
Sheoran et al. Methods of denoising of electroencephalogram signal: A review
Rashmi et al. ECG denoising using wavelet transform and filters
Walters-Williams et al. Performance comparison of known ICA algorithms to a wavelet-ICA merger
de Cheveigné Time-shift denoising source separation
CN111631710A (zh) 一种状态相关的动态脑电信号中肌电伪迹的消除方法
Dhull et al. EEG Artifact Removal Using Canonical Correlation Analysis and EMD-DFA based Hybrid Denoising Approach
Agounad et al. Intelligent fuzzy system for automatic artifact detection and removal from eeg signals
Feng et al. A Novel SSA-CCA Framework forMuscle Artifact Removal from Ambulatory EEG
CN114098763B (zh) 一种脑电去噪方法
Khezri et al. Surface electromyogram signal estimation based on wavelet thresholding technique
Nguyen et al. A deep sparse autoencoder method for automatic EOG artifact removal
CN111493864B (zh) Eeg信号混合噪声处理方法、设备和存储介质
Pramudita et al. Removing ocular artefacts in EEG signals by using combination of complete EEMD (CEEMD)—Independent component analysis (ICA) based outlier data
Khatter et al. Study of various automatic eeg artifact removal techniques

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