级联自适应局部投影降噪方法

徐礼胜, 崔慧颖, 吴俊鼎, 王仲怡

(1. 东北大学 医学与生物信息工程学院, 辽宁 沈阳 110169; 2. 沈阳东软智能医疗科技研究院有限公司, 辽宁 沈阳 110167)

信号在采集过程中易受不同程度的噪声干扰,这些噪声破坏信号结构,降低了识别信号特征的准确性[1];因此,噪声抑制是信号研究的基础[2].当信号频谱与噪声频谱存在明显差别时,线性滤波器效果显著[3-4].而产生自混沌系统的信号通常具有宽频谱的特性,信号频谱与噪声频谱重叠时,传统线性滤波方法会导致信号失真.

基于非线性动力学理论的局部投影算法(local projection, LP)为宽频谱信号的降噪提供了一种新的选择.期望信号与噪声在相空间的不同维度子空间中具有不同的结构,局部投影算法则据此区分两者.依据Takens嵌入定理[5]重构信号得到的嵌入空间即为相空间,相空间中的光滑非线性流形反映系统状态,因此系统产生的期望信号在相空间中的轨迹位于或非常靠近光滑非线性流形[6],而噪声的相空间轨迹因噪声并不具备任何来源约束,因此分散在流形四周,甚至可能位于流形上.LP算法在重构相空间中为每个相点划分邻域,对邻域内的非线性流形进行局部线性逼近.依据特征值分解理论区分信号空间和噪声空间,通过向噪声子空间投影滤去噪声成分.该算法目前已被成功应用于呼吸音信号滤波[7]、核磁共振激光数据降噪[8]、脑机接口脑电图增强[9]、语音信号分析[10-11]等.

自适应邻域选取是局部投影算法的研究重点之一.Przybya等[12]使用K-means聚类算法形成邻域,取得了较好的降噪效果.Chelidze[13]重构了信号的短时轨迹链束,并利用平滑正交分解识别局部子空间,降噪效果优于基于传统正交分解的降噪方法.经验模态分解(empirical mode decomposition,EMD)作为一种完全由数据驱动的方法,不需要任何先验知识即可自适应地将复杂数据分解为有限个本征模态函数,适用于非线性时间序列的处理[14-16].因此,本文提出基于经验模态分解的邻域选取策略,根据自相关函数确定分解得到的噪声分量,再通过评估噪声分量在原始信号中的占比确定邻域半径.此外,目前有许多结合小波理论对信号进行局部投影降噪处理的研究[17-19],Wan等[20]通过将局部投影算法应用于脑电图的小波变换,去除了噪声干扰.小波阈值去噪方法理论成熟,降噪效果明显,但其在恢复信号非线性流形结构方面性能不佳.LP算法通过将噪声分量投影到低维子空间中以减少其与非线性流形的偏差,能够最大程度地恢复信号的真实流形结构.综上,本文提出局部投影算法与小波阈值去噪级联的降噪方法.

1 方 法

1.1 局部投影降噪算法

假设x(t)(t=1,2,…,l)是动力系统输出的一组观测时间序列.Takens嵌入定理已经证明,当嵌入维数m>2D2(D2为系统的关联维数)时,从x(t)中每隔时间τ取出一个观察值作为时间序列的状态测量值,即可重构出与原动力系统具有相同拓扑特性的m维相空间.重构得到的信号为

X(t)=[x(t),x(t+τ),…,x(t+(m-1)τ)]T.

(1)

其中,t=1,2,…,l-(m-1)τ.

理论上,当嵌入维数大于或等于m时,系统的动力学特性即可得到充分展现[21].但是,更大的嵌入维数也意味着更大的计算成本.因此,选取一个大于2D2的最小嵌入维数是必要的.基于嵌入窗思想的C-C算法[22]以计算量小、结果准确的特点成为目前最常用的计算嵌入维数和时间延迟的方法.

在重构的相空间中,原始动力系统产生的混沌吸引子仅存在于维数为m0的低维子空间,m0维子空间的正交补空间叫做零子空间,其维数为Q=m-m0.信号无噪声的情况下,零子空间中无元素;当信号含噪声时,噪声因其自身的无方向性存在于任何子空间.显然,此时零子空间中的元素全部由噪声产生.因此,把邻域内的所有相点投影到m0维零子空间中,即可将低维动力学成分从高维空间中分离出来,使混沌吸引子流形回归到真实形态.

1.2 基于经验模态分解的邻域选取策略

邻域半径是保证局部投影算法降噪质量的重要参数[23].半径太小会导致部分相点丢失,使邻域内的期望信号被噪声淹没;半径太大则会包含过多无用相点,导致局部线性逼近的结果不准确.为避免凭经验选取邻域半径的盲目性,本文利用经验模态分解方法估计观测时间序列的噪声水平,据此选取半径进行邻域的划分.

经验模态分解方法认为任何时间序列都是由多个具有不同特征尺度的本征模态函数(intrinsic mode function,IMF)组成的,即

(2)

式中:x(t)为原始信号;imfi(t)(i=1,2,…,n)为分解得到的一组IMF;rn(t)为残余分量.IMF分量的频率从高到低代表了x(t)的不同成分,而噪声通常表现出高频特性.确定代表噪声的IMF分量后,即可评估信号的噪声水平,从而选取合适的邻域范围.

自相关函数是反映信号与其自身在不同时间尺度上相似程度的统计量,随机噪声在不同时间尺度上的随机性使其与自身表现为弱相关性,而一般信号则表现为较强相关性,因此可用自相关函数区分随机噪声和一般信号.

(3)

其中,M是重构噪声分量的绝对中值偏差，

(4)

定义邻域半径为

(5)

其中,σdata是实测时间序列的标准差,即根据估得的噪声标准差与原时间序列的标准差之比确定邻域半径.

1.3 级联自适应局部投影降噪算法

局部投影降噪算法对原始信号的非线性流形有明显的修复作用,但对高频噪声的去除效果不如频域降噪方法[24].如图1所示,脉搏信号在经过两次局部投影降噪处理后仍残留部分噪声.

图1 含噪脉搏信号和两次局部投影降噪后信号

传统局部投影方法通过迭代的方式去除残留噪声,实际上,2～3次迭代之后信号的信噪比提升量会逐渐收敛.如图2所示,持续迭代使信号平滑的代价是幅值的持续衰减,这种丢失信息的降噪方法显然是不可取的.

本文方法只进行1～2次局部投影处理,这样可以尽可能多地保留信号非线性特征,并且信噪比得到了大幅提升,而小波阈值方法正适用于中等信噪比的信号降噪处理.虽然小波阈值方法能够较好地抑制高频噪声,但恢复混沌吸引子流形结构的能力有限[25].因此,本文综合两种方法的优点,1～2次局部投影后再通过小波阈值方法去除残留的少量噪声成分即可达到既保留信号特征又平滑噪声的目的,本文的降噪方法流程如图3所示,具体步骤如下:

1) 相空间重构.观测时间序列在m维重构相空间中可表示为

X=[X(1),X(2),…,X(N)],N=l-(m-1)τ.

(6)

2)邻域划分.对相空间中的任意相点X(i),定义其邻域为

Ui={X(j)‖X(j)-X(i)‖

(7)

其中:X(j)为满足条件的邻点;‖·‖表示欧几里得距离;r为邻域半径,由前述基于经验模态分解方法得到.

3)局部投影.定义协方差矩阵：

(8)

图2 原始信号与三次局部投影降噪后的信号对比

计算C(i)的特征向量和对应的特征值.理论上,特征向量对应的方向可以衡量信号与噪声的不同方向.较大特征值对应的特征向量所在的方向代表信号子空间;反之,较小的Q个特征值λ1,λ2,…,λQ对应的特征向量eq(q=1,…,Q)所在的方向代表了由噪声组成的零子空间[8].引入权矩阵R对坐标进行修正以避免边界效应,定义r11=rmm=1 000,其余元素全为1[26].各个相点在m0维子空间中的投影即为噪声分量,

(9)

其中,为避免各个相点首尾坐标畸变带来的误差,修正后的相点可表示为

(10)

4) 更新所有相点,对修正后的所有相点进行反重构.

5) 对反重构得到的信号进行小波阈值降噪处理,抑制高频噪声.

图3 级联降噪方法框图

2 结果与讨论

2.1 Lorenz系统时间序列实验

本节以Lorenz系统产生的混沌时间序列检验所提算法的有效性,Lorenz系统的动力学方程为

(11)

选取σ=10,R=28,b=8/3,设定积分步长h=0.003,此时系统表现为混沌状态，采用四阶龙格-库塔法求解方程.取x分量上的4 000个数据点作为实验信号,以添加噪声后信噪比为20 dB的信号为例,首先对仿真信号进行经验模态分解,并计算邻域半径用于局部投影降噪.图4分别为原始信号、含噪信号、小波阈值方法降噪信号和本文方法降噪信号的相空间重构图.图4b显示,信噪比为20 dB时吸引子流形已经完全被噪声破坏,无法分辨.图4c展示了小波阈值去噪方法的处理结果,高频噪声虽然得到了抑制但吸引子流形的结构较为混乱.当嵌入维数m=11,时间延迟τ=1,邻域半径r=0.592 8,并选用“db4”小波基时,本文方法的降噪结果如图4d所示,信号的相空间轨迹得到了基本恢复,并且流形结构更加有序,显著改善了单独使用小波阈值方法的不足.特别地,对于吸引子流形中轨迹密集的位置,本文方法的处理结果具有更清晰的细节特征,更加接近原始信号,这也说明利用本文方法选取的邻域可以提高逼近非线性流形的准确率.

为定量评估算法性能,向Lorenz信号中添加不同信噪比的噪声,分别计算小波阈值降噪与本文方法降噪后信号的信噪比(signal to noise ratio,SNR)和均方误差(mean square error,MSE).通常来说,SNR越大MSE越小的降噪后的信号越接近原始信号.结果如表1所示,相较于小波阈值去噪方法,本文方法在中高强度的噪声范围内能获得更高的SNR和更低的MSE,说明该方法在滤除噪声的同时能够最大程度地接近原始信号.

表1 不同噪声强度下Lorenz信号降噪后的信噪比和均方误差

2.2 脉搏信号实验

分别采集桡动脉、肱动脉、颈动脉三处的脉搏波.为对比评估本文算法在脉搏信号降噪中的性能,将三处实测信号降噪后得到的信号作为干净信号,采用SNR和MSE作为评价指标,分别使用小波阈值去噪方法、本文方法对添加不同强度噪声的三处脉搏信号进行降噪处理.图5结果显示本文方法整体优于小波阈值去噪方法,三处脉搏信号的实验结果均显示本文方法能够获得更大的信噪比和更小的均方误差.噪声能量较大时,信号几乎完全被噪声覆盖,两种方法性能差别较小,都不能很好地恢复信号.噪声能量相对较小时,信号结构受噪声影响较小,局部投影方法恢复信号结构的性能发挥程度有限,而小波阈值去噪方法对此类噪声敏感性更强.因此,随着SNR的提高,两种方法在降噪效果上的差别逐渐变小.中等噪声强度下,本文方法的优势更加明显.此时,含噪信号经过局部投影处理,信号基本结构得以恢复,为小波阈值去噪方法平滑噪声的处理提供了前期基础.相比单独使用小波阈值去噪方法,本文方法明显提高了信噪比,降低了均方误差.其中,本文方法在颈动脉脉搏信号的降噪处理中优势更加明显.颈动脉脉搏波波形面积大,主波下降速度慢,潮波、重搏波位置更高,信号波幅变化相对较小,主波持续时间更长,与桡动脉、肱动脉脉搏波形态差别较大.“db4”小波基无法很好地匹配该波形,致使小波阈值方法性能有所下降.本文方法在不同部位的脉搏信号表现出稳定的降噪效果.

图4 Lorenz信号降噪前后相空间重构图

图5 两种方法在不同位置脉搏波上的降噪结果对比

图6展示了对实测脉搏信号分别使用小波阈值去噪方法、本文方法的相空间重构结果和时域图像,其中时间延迟τ=1、嵌入维数m=6、邻域半径r=1.369 7,选用“db4”小波基.原始信号如图6a所示,从时域图像中尚可观察到信号的部分周期特征,但相空间轨迹已无法分辨.对比图6b与图6c,小波阈值去噪方法虽然使信号的重构轨迹变得平滑,但结构杂乱,轨迹堆叠现象明显,时域图像中各周期波形形态差异较大;而本文方法在恢复信号相空间基本结构的基础上去除了残留的高频噪声,得到的信号曲线光滑,轨迹更加收敛,混沌吸引子呈现出有规律的几何结构,且信号时域波形规整有序.实测脉搏信号降噪结果说明本方法恢复信号原始特征的能力明显优于小波方法.

图6 实测脉搏信号降噪前后时域和相空间重构图像

2.3 心电信号实验

为验证改进方法在生理信号降噪方面的普适性,本节选择实测心电信号(electrocardiogram signal, ECG)作为验证对象.心电信号在采集过程中易受噪声干扰,噪声会对心电信号QRS波的检测产生不良影响,降低计算机自动诊断的准确性[27].心电信号中的噪声主要有基线漂移、工频干扰和肌电干扰[28].其中,肌电干扰的频谱范围较宽,与ECG频谱存在重叠,传统方法很难将此类噪声与心电信号完全分离.

图7a中原始心电信号受到噪声干扰.P波、Q波、S波和T波被覆盖,其对应的相空间轨迹杂乱重叠,无法分辨.小波方法对心电信号中的噪声有一定的抑制作用,图7b中的信号波形得到了大致恢复,但仍存在明显瑕疵,不同周期中的P波、Q波、S波和T波恢复程度差异较大,更换不同小波基所得结果并无明显差别,说明小波方法的结果对小波基的选取依赖性较大.利用本文方法对心电信号作降噪处理,取时间延迟τ=1、嵌入维数m=9、邻域半径r=0.257 1,选用“db4”小波基.由图7c可知,本文方法对心电信号中的噪声有明显抑制作用且降噪效果稳定可靠,并不受信号形态影响,降噪后的心电信号完全平滑,形态稳定,各个周期波峰清晰.相空间轨迹结构稳定,且幅值较小的波所对应的轨迹同样得到了基本恢复,规整有序.

图7 实测心电信号降噪前后时域和相空间重构图像

3 结 语

具有非线性、非平稳性特征的微弱信号易受噪声干扰,而传统线性滤波器在期望信号频带与噪声频带发生重叠时性能不佳.针对此类问题,本文提出了一种基于EMD的噪声自适应邻域选取策略和级联局部投影降噪方法.通过向Lorenz信号中添加不同程度高斯白噪声验证提出方法的有效性,不同降噪方法下含噪Lorenz信号去噪前后的SNR和MSE结果表明,本方法能够有效去除强高斯白噪声,并在恢复信号非线性动力学特性方面更有效.心血管解剖结构的分形特点及生物电传导系统的分形结构使与之密切相关的生理信号表现出混沌现象与非线性动力学行为.因此,本文量化分析了所提方法在脉搏信号和心电信号中的降噪效果.脉搏信号仿真实验的定量结果表明,与小波阈值降噪方法相比,本文方法抑制噪声的能力有明显提高.实测心电信号降噪结果也证实本文方法消除噪声干扰的能力显著优于小波方法.降噪前后的信号相空间重构结果表明本文方法能够在中、高强度噪声下恢复信号的非线性结构,保留信号特征.