改进的双变量收缩函数模型脑电信号消噪方法*

陈顺飞，罗志增，周镇定

（杭州电子科技大学智能控制与机器人研究所，杭州310018）



改进的双变量收缩函数模型脑电信号消噪方法*

陈顺飞，罗志增*，周镇定

（杭州电子科技大学智能控制与机器人研究所，杭州310018）

摘要：针对传统小波消噪全局阈值处理独立性假设和双变量函数模型对没有父系数的最高层小波系数不做处理的缺陷，提出一种高密度离散小波变换中利用双变量收缩函数对脑电信号进行消噪的方法。子小波系数根据双变量函数实现局部自适应收缩处理。同时根据父系数趋于0时，阈值函数近似于软阈值函数，对最高尺度小波系数进行软阈值法消噪。从实际信号处理效果和客观定量指标两方面进行评价，结果表明这种改进算法都优于软阈值法、硬阈值法以及双变量收缩法。

关键词：脑电信号；信号消噪；高密度小波变换；双变量收缩函数

脑电信号EEG（Electroencephalogram）是由头皮表面大量神经元突触后电位同步综合而形成的，反映大脑运行状态和神经细胞活动情况的生物电信号。因此，对于脑电信号的分析可以获得大量的生理、心理及病理信号［1-2］。脑电信号是一种非平稳非线性极其微弱的随机信号，其幅值是μV级，极易受到噪声干扰，如工频干扰、静电干扰、眼电伪迹、心电伪迹等［3-4］。去除脑电信号中的噪声信号，提升脑电信号的信噪比是脑电信号处理中的首要环节。小波变换由于其具有良好的时频局部化特性、多分辨率特性、解相关性、选基灵活性［5］等特点，在生物医学信号消噪领域中得到广泛应用［6］。小波阈值消噪法［7］将信号通过小波分解转化为一系列的分解系数，对不同尺度下的小波系数作阈值收缩处理，可以起到消噪作用。常用的阈值消噪方法有：软阈值法［8］、硬阈值法［9］等。传统的小波阈值消噪方法都假定不同尺度上的小波系数之间是相互独立的，忽视了父小波系数与子小波系数之间的强相关性［10］。信号与噪声的小波变换在各尺度下的不同传播特性表明，父系数与子系数之间的强相关性恰恰是普遍存在的。因此基于独立性假设的阈值消噪算法将导致原本相关的信息被重复处理，使信号消噪过程中丢失一些重要的细节特征，不利于信号奇异点的检测，且易出现伪吉布斯现象，给进一步分析带来误差［11］。针对传统小波阈值消噪方法的不足，本文提出一种改进的双变量收缩函数模型消噪方法。在高密度离散小波变换出色的信号特征表达能力的基础上，结合双变量收缩函数［12-14］具有考虑父小波系数与子小波系数之间强相关性的特性，实现小波系数阈值的局部自适应收缩处理，避免了由于独立性假设导致的重复处理。考虑到传统的双变量收缩函数模型缺乏对最高层小波系数的处理，而最高层小波系数虽然处于较低频带，但同样含有噪声［15］。本文提出的方法将对最高层小波系数采用软阈值法进行阈值处理。改进算法在脑电信号消噪处理实验应用中取得了较好的消噪效果。

1　双变量收缩函数

假设脑电信号是被方差为σ2n的加性高斯白噪声所污染，则有：

式（1）中，n为高斯白噪声，x为脑电信号。消噪的目的就是要从含噪信号s中估计出脑电信号x，将含噪信号s由时域变换到小波域，分解得到多尺度的小波系数，如式（2）。

式（2）中的y代表含噪信号的小波系数，w为脑电信号所对应的小波系数，n为噪声的小波系数。由于双变量收缩函数需考虑小波系数及其父系数的相关性，故设w1k表示第k层的小波系数，w2k表示下一层中与w1k位置相对应的父节点小波系数，可得向量表示为：

其中yk=（y1k，y2k），wk=（w1k，w2k），nk=（n1k，n2k），m为子小波系数的数目。采用最大后验概率公式估计未受噪声污染的信号的小波系数w^ (yk)，记为式（4）。

由贝叶斯概率估计作进一步推导之后，可将式（4）转换为式（5）。

由式（5）可知，只有分别求得噪声系数的概率密度函数pn（n）和目的信号系数的概率密度函数pw（wk），才可得到真实信号的系数估计(yk)。假设各层噪声信号n1k，n2k为独立同分布的高斯变量，因此n的概率密度函数pn（n）可表示为：

而不同尺度下小波系数的概率密度函数pw（wk）可表示为［12］：

其中(g)+={0 g＜0 g g≥0。式（8）表明真实信号第k层小波系数的估计值不仅与子小波系数y1k有关，还与父小波系数y2k、第k层的信号方差σ2k和噪声信号方差σ2n有关，使系数收缩处理具有较强的局部自适应能力。当父系数较大时，小波系数估计值的收缩比例较小，而当父系数趋于零时，公式计算得到的阈值函数就近似于软阈值函数。估计值的双变量收缩阈值函数如图1所示。结合式（8）与图1的结果表明，当含噪子小波系数和含噪父小波系数的平方和小于阈值时，估计值为0。由于小波变换后较小的小波系数中含有多数噪声，因此双变量收缩阈值处理可以有效的滤除噪声。同时，图1中子小波系数的估计值小于相对应的含噪子小波系数y1k，实现了小波系数的收缩变换。传统的双变量收缩函数模型对于没有父系数的最高层小波系数不做任何处理，最高层小波系数虽然处于较低频段，但是同样含有噪声［15］。从上述的分析中得知，当父系数趋于零时，阈值函数类似于软阈值函数。因此本文提出的算法模型在双变量函数模型的基础上，增加对最高层小波系数的软阈值处理。

图1　?的双变量收缩函数图

2　高密度离散小波变换

高密度离散小波变换HDDWT（High Density Discrete Wavelet Transform）［16-17］是一种改进的小波变换，相比于一代小波分析，它具有近似的平移不变的特性，更出色的信号表达能力，能够在时频域上更加细致地刻画信号，拥有更佳的信号降噪性能。

农业始终是我国重要的产业之一，我国人民就没有放松对农业生产技术的追求和钻研，近些年，规模化、现代化、产业化的农业已经逐渐代替传统农业生产方式，科学技术对于农业生产的重要性日益显著。玉米作为我国重要的粮食作物和经济作物，相较水稻和小麦，玉米在耐旱等方面更具优势，山区气候和地质条件对玉米的产量和质量有更高的需求，我国近些年一直在研究山区玉米高产栽培技术，并进行推广应用。

HDDWT最大特点是具备间尺度分析的特性，即第二个母小波是按1/2的倍数进行平移，这使得HDDWT在时频域上的采样密度为非冗余小波变换的三倍。基于高时频采样密度的HDDWT在降低对信号平移敏感性的同时获得了出色的时频域分辨能力。HDDWT对信号进行分解重构示意图如图2所示。它由三个通道滤波器组来实现。图2中hi（n）（i=0，1，2）分别表示小波的低通滤波器、带通滤波器和高通滤波器。

图2　HD-DWT对信号进行分解重构示意图

HDDWT尺度系数与小波系数分别由式（12）～式（14）计算得到

其中，di，j（n）表示第j级分解后得到的第i个小波系数。同时HDDWT的尺度函数与两个母小波函数需分别满足以下的约束关系。

式中，Φ（t）为基函数，ψ1(t)、ψ2(t)分别代表两个小波函数，值得注意的是ψ2(t)是以1/2的整数倍作平移的。根据多采样率信号处理的关系，HDDWT若要实现重构需满足表达式（18）和式（19）。

其中Hi（z）为hi（n）（i=0，1，2）的z变换，其中H2（z）可通过H0（z）和H1（z）因式分解得到。如果达到上述重建条件，且Φ（t）足够正则，则信号重构的表达式为：

本文基于双变量收缩函数的局部自适应能力和高密度离散小波变换出色的信号描述能力，对脑电信号进行消噪分析。具体算法实现步骤如下：

Step 2对任意的一列小波系数y1k找到其对应的父系数y2k。

Step 3根据公式（8）求得对应小波系数收缩后的估计值，并以代替原始的小波系数，完成消噪处理。同时对最高尺度的小波系数进行软阈值法消噪处理。

Step 4对各尺度的系数自下而上进行重构，即可得到消噪后的脑电信号。

3　实验结果分析

3.1实验设计与结果

本文应用美国NeuroScan公司的scan4.3脑电采集系统对脑电信号进行采集实验，设备采样频率为250 Hz，精度为32 bit，并选择陷波滤波滤除50 Hz工频干扰。本实验以3位身体健康的在校男性研究生作为受试者，在受试者意识清醒的情况下进行脑电实验。脑电电极按照国际标准导联10/20系统放置，以左侧乳突和右侧乳突为参考电极。本文选择典型样本C3通道进行脑电信号消噪实验分析。本实验通过采集想象左手、右手、双脚、舌头所产生的脑电信号进行消噪处理。每次实验的持续时长为9 s：0～2 s，受试者处于想象运动准备状态；2 s时，系统发出“叮”提示音，提示受试者集中注意力；3 s～6 s时，受试者根据屏幕上出现上（↑）、下（↓）、左（←）、右（→）四个方向的诱导图片分别想象双脚脚趾前曲、舌头后卷、左手挥动和右手挥动这4种运动，直至第6 s时系统发出“嘟”的提示音，受试者结束想象，放松调整。一次完整实验的采样时序如图3所示。完成4种想象运动实验各采集4组脑电数据，每组进行10次实验，总共160次实验数据。

图3　实验设计流程时序图

选取一组实验中C3通道运动想象部分（第3 s～7 s）的脑电信号波形进行消噪处理。因为db小波函数在时频域上都有较好的局域性，所以选择db5作为小波分解的母小波进行对比实验分析。实验通过多次对比发现，当小波分解的层数为3层时，信号消噪效果较好。将原始信号分别经过小波软阈值、小波硬阈值以及改进的BSF-HDDWT法进行消噪效果对比分析，分析结果如图4所示。从图4可以看出三种算法都取得一定的消噪效果。虽然经前两者算法处理后的信号噪声得到了明显的消除，但原始脑电信号中一些细节信息与边缘特征亦被过滤，出现过平滑的现象。而采用改进的BSF-HDDWT方法消噪后的脑电信号波形相对清晰，更重要的是原始信号的细节特征也被很好地保留下来。

为了定量分析改进算法的优越性，利用两个常用的去噪指标均方根误差RMSE与信噪比SNR来评价信号的消噪效果。对160次实验数据求得3种消噪算法的RMSE和SNR，并取均值和标准差如表1所示。根据RMSE与SNR的评判标准可知，RMSE的值越小，SNR的值越大表示信号的消噪效果越好。从表1可以看出，改进的BSF-HDDWT法在SNR和RMSE方面都得到较好的提高，从定量上说明改进的BSF-HDDWT算法优于其他2种算法。

图4　多种消噪方法消噪效果对比图

表1　3种消噪算法的RMSE和SNR结果对比

3.2仿真实验与结果分析

为进一步验证改进BSF-HDDWT算法的消噪效果，用Matlab构造一个采样频率为250 Hz，采集时长为6 s的仿真信号S（t），如式（21）所示。

在仿真信号中添加信噪比为5的高斯白噪声，对混合后的信号利用软阈值法、BSF-HDDWT法、以及改进的BSF-HDDWT法进行消噪分析，对消噪后的信号进行傅里叶变换转化到频域如图5所示。从图5中可以看出，软阈值法虽然滤除了大部分噪声，但是由于软阈值法阈值选取的全局性缺陷，忽略小波系数层与层之间的相关性，造成对细节信号多次重复处理，使得仿真中的高频成份40 Hz的信号被严重削弱，这会造成信号的失真现象严重。BSF-HDDWT法相比于软阈值法在细节保持方面得到明显改善，从图中可以看出信号的三个峰值，但是由于该方法没有对最高层的小波系数进行处理，因此在低频信号存在较多的噪声，如14 Hz～40 Hz之间的噪声较多，而在实际脑电信号分析中，14 Hz～40 Hz频率段含有丰富的脑电信号特征。而本文对BSF-HDDWT的改进算法，相比于前两者，既可以保留更多的细节信号，又可以滤除最高层小波系数中的低频噪声信号。

对信噪比为1、5、10加噪后的信号进行分析，用小波软阈值法、小波硬阈值法、小波BSF法、BSFHDDWT法、改进的BSF-HDDWT法分别进行消噪处理，评价结果如表2所示。

图5　多种消噪方法的消噪效果频域对比图

表2　5种消噪方法的去噪评价结果对比

结合图5和表2，小波软阈值法和硬阈值法由于全局性的阈值处理，使得幅值为5 V，频率为40 Hz的信号分量被严重削弱，因此SNR和RMSE的评价效果较差。而以双变量收缩函数为基础的3种算法的消噪结果则明显优于上述两种方法。双变量收缩阈值函数利用各层小波系数的局部方差对消噪阈值大小进行调整，同时考虑父系数与子系数之间的强相关性对小波系数作收缩处理，使小波阈值消噪具有局部自适应的能力，进而更好的保留信号的细节特征。从表2中可以看出，相比于一代离散小波，具备近似平移不变特点的HDDWT能更细致地刻画信号的时频域特点，使处理之后的信号更接近于目标信号的真实波形特征，具有更佳地信号消噪能力。对比BSF-HD⁃DWT法与改进的BSF-HDDWT法的SNR和RMSE评价结果，当信噪比为1和5时，由于噪声干扰较强，改进的BSF-HDDWT法能够改善对最高层小波系数的阈值消噪处理，因此具有相对BSF-HDDWT法较好的消噪效果。而当信噪比为10时，BSF-HDDWT法的消噪效果则优于改进算法，这是因为此时噪声干扰较小，改进的BSF-HDDWT法会对信号进行过度处理，导致消噪效果变差。从中可以得出结论，本文改进的BSF-HDDWT法比较适合频带复杂的强噪声信号的消噪处理。由于实验采集到的脑电信号是一种极其微弱且噪声较强的信号，因此可以采用本文改进的BSF-HDDWT法对其进行消噪处理。

4　结论

本文针对脑电信号微弱且频域特征复杂、噪声干扰强的特点，引入对不同尺度的小波系数有局部自适应能力，能有效刻画信号细节特征的双变量收缩阈值函数，并结合拥有出色信号表达能力的高密度离散小波变换对脑电信号进行消噪处理。同时对含有低频噪声的最高尺度小波系数进行软阈值法消噪处理。利用多种消噪方法进行实验对比分析，通过实际信号与仿真信号的处理效果观察和SNR、RMSE客观指标进行评价，实验结果表明，本文改进的BSF-HDDWT法具有更强的噪声抑制能力，同时可以更好地保留信号的细节特征。良好的消噪效果和保留信号细节的能力，为后续脑电信号的特征提取与和模式识别打下坚实的基础。

参考文献：

［1］马玉良，许明珍，佘青山，等.基于自适应阈值的脑电信号去噪方法［J］.传感技术学报，2014，27（10）：1368-1372.

［2］Mporas I，Tsirka V，Zacharaki E I，et al.Seizure Detection Using EEG and ECG Signals for Computer-Based Monitoring，Analysis and Management of Epileptic Patients［J］.Expert Systems with Applications，2015，42（6）：3227-3233.

［3］王兵，王魁，梁晓霖，等.脑电信号中工频干扰去除的综合研究［J］.传感技术学报，2010，23（1）：87-92.

［4］吴明权，李海峰，马琳.单通道脑电信号中眼电干扰的自动分离方法［J］.电子与信息学报，2015，37（2）：367-372.

［5］潘泉，孟晋丽，张磊，等.小波滤波方法及应用［J］.电子与信息学报，2007，29（1）：236-242.

［6］罗志增，周镇定，周瑛，等.双树复小波特征在运动想象脑电识别中的应用［J］.传感技术学报，2014，27（5）：575-580.

［7］吴光文，王昌明，包建东，等.基于自适应阈值函数的小波阈值去噪方法［J］.电子与信息学报，2014，36（6）：1340-1347.

［8］Donoho D L.De-Noising by Soft-Thresholding［J］.IEEE Transac⁃tions on Information Theory，1995，41（3）：613-627.

［9］Zhang X，Feng X，Wang W，et al.Image Denoising via 2D Diction⁃ary Learning and Adaptive Hard Thresholding［J］.Pattern Recog⁃nition Letters，2013，34（16）：2110-2117.

［10］Simoncelli E P.Bayesian Denoising of Visual Images in the Wave⁃let Domain［M］.New York：Springer，1999：291-308.

［11］Cho D，Bui T D.Multivariate Statistical Modeling for Image De⁃noising Using Wavelet Transforms［J］.Signal Processing：Image Communication，2005，20（1）：77-89.

［12］Sendur L，Selesnick I W.Bivariate Shrinkage Functions for Wave⁃let-Based Denoising Exploiting Interscale Dependency［J］.IEEE Transactions on Signal Processing，2002，50（11）：2744-2756.

［13］Min D，Jiuwen Z，Yide M.Image Denoising via Bivariate Shrink⁃age Function Based on a New Structure of Dual Contourlet Trans⁃form［J］.Signal Processing，2015，109：25-37.

［14］Sethunadh R，Thomas T.Spatially Adaptive Image Denoising Us⁃ing Inter-Scale Dependence in Directionlet Domain［J］.IET Im⁃age Processing，2014，9（2）：142-152.

［15］李向军，姜玉莉.一种改进的双变量收缩模型图像去噪［J］.现代电子技术，2014（8）：132-134，137.

［16］Qin Y，Tang B，Wang J.Higher-Density Dyadic Wavelet Trans⁃form and Its Application［J］.Mechanical Systems and Signal Pro⁃cessing，2010，24（3）：823-834.

［17］Selesnick I W.A Higher Density Discrete Wavelet Transform［J］.IEEE Transactions on Signal Processing，2006，54（8）：3039-3048.

陈顺飞（1992-），男，浙江丽水人，硕士研究生，研究方向为生物医学信息检测、信息融合与信息处理，142060149@hdu.edu.cn；

周镇定（1989-），男，硕士研究生，研究方向为生物医学信息检测、信息融合与信息处理。

罗志增（1965-），男，浙江慈溪人，博士生导师，1998年在浙江大学获博士学位，主要从事智能机器人技术、传感器及多信息融合、生物医学信息检测与利用等领域的研究，luo@hdu.edu.cn；

EEG Denoising Method Based on Improved Bivariate Shrinkage Function*

CHEN Shunfei，LUO Zhizeng*，ZHOU Zhending
（Robot and Intelligent Control Research Institute，Hangzhou Dianzi University，Hangzhou 310018，China）

Abstract：Traditional wavelet denoising methods have an assumption that the wavelet coefficients are independent in global thresholding.Then traditional bivariate shrinkage function model has a deficient in not considering the highest scale wavelet coefficients with no parents coefficients.To tackle these defects，in this paper an EEG denois⁃ing method is proposed based on high density discrete wavelet transform using bivariate shrinkage function.In this method，the children coefficients will achieve local and adaptive shrinking treatment using bivariate function.Then because of the appearance when the parents coefficients tend to zero，the threshold function approximate to the soft threshold，the soft threshold is used for denoising in the highest scale wavelet coefficients.The results show that this improved method is better than the soft threshold，the hard threshold method and bivariate shrinkage method from two perspectives，the actual signal processing effect and the objective quantitative indicators.

Key words：EEG signal；signal denoising；high density discrete wavelet transform；bivariate shrinkage function

doi：EEACC：523010.3969/j.issn.1004-1699.2016.02.016

收稿日期：2015-08-27修改日期：2015-10-14

中图分类号：TP391

文献标识码：A

文章编号：1004-1699（2016）02-0242-06

项目来源：国际自然科学基金项目（61172134，61201302）