基于模糊递归和最优硬阈值的局部投影降噪算法

王 东, 崔忠马, 陈文东, 舒 勤,*

(1. 四川大学电气工程学院, 四川 成都 610065; 2. 北京遥感设备研究所, 北京 100084)

0 引 言

实际测得的信号易受噪声干扰,从而破坏信号的结构特征,因此降噪是后续各种信号处理的基础[1-2]。传统信号降噪方法基于线性系统的假设,认为信号和噪声的频谱不完全重合。但是，由混沌系统产生的信号往往具有宽频谱、内在伪随机等混沌特性,信号和噪声频谱重合,使得传统线性滤波方法失效[3]。

Cawley等[4]和Sauer等[5]提出的基于混沌动力学理论的局部投影降噪算法为具有宽频谱特性的混沌信号降噪提供了一种新方法。该方法首先根据Takens定理重构信号,得到一个与原始动力系统微分同胚的相空间[6]。期望信号的吸引子被限制在一个低维流形上,而噪声的吸引子分散在流形周围。局部投影降噪算法根据信号与噪声在相空间中局部动力学特性与局部几何特性的不同，区分信号子空间和噪声子空间,利用几何投影去除噪声分量,再将相空间反重构为时间信号,从而达到降噪的目的。在实际情况中，许多信号都有混沌特性,目前该算法成功应用的场景包括人类语音信号处理[7]、振动信号分析[8]，还包括生物信号处理,如脑电信号处理[9]、心电信号处理[10]、呼吸声音信号处理[11],以及激光数据处理[12]、故障检测[13]等。

邻域选取和子空间划分是局部投影算法的两个研究重点。邻域选取的方法很多,但都存在不足。Cawley等[4]和Sauer等[5]直接指定邻点数目,这种方式易受人为因素影响。Kantz等[12]利用递归图估计邻域半径,效果优于原始方法,但计算复杂。冯飞龙等[14]利用小波分解估计初始邻域半径,再进行邻点搜索,徐礼胜等[15]采用经验模态分解(empirical mode decomposition, EMD)估计邻域半径[15],但这些方法都需要预先估计噪声水平。Przybya等[16]和Kotas等[17]利用K-means聚类算法确定邻域,但聚类数难以确定。对于子空间划分,现有方法直接指定子空间的维数,或根据特征值的大小对子空间进行划分。Chelidze等[18]通过重构信号的短时轨迹,利用平滑正交分解识别子空间。但考虑到混沌吸引子本身具有分数维的性质[19],系统局部邻域的动力学特性以及每个邻域内的噪声分量占据的空间也不尽相同,且在实际情况中可能并不知道原始动力系统的相空间维数,所以每个局部邻域应该进行不同的子空间划分。

针对以上问题,本文提出基于模糊递归图与最优硬阈值准则的局部投影降噪算法。首先，根据模糊递归图对邻域进行选择。为避免计算邻域协方差矩阵,直接将邻域矩阵进行奇异值(singular value decomposition, SVD)分解；然后，根据最优硬阈值对局部邻域的信号子空间和噪声子空间进行划分,避免人为因素的影响；最后，针对高斯白噪声,采用本文所提方法分别对仿真Lorenz信号与实测含噪心电图(electrocardiogram, ECG)信号进行仿真研究,并与其他局部投影算法以及其他ECG信号降噪方法进行对比,仿真结果验证了本文方法的有效性。

1 局部投影算法原理

1.1 相空间重构

根据Takens嵌入定理,对于无限长、无噪的d维混沌吸引子的标量时间信号{x(t)},在拓扑不变的意义下可以找到一个m维的嵌入相空间。其中,m≥2d+1。而对于有限长、含噪的信号，可采用坐标延迟重构[20]。

设x1,x2,…,xN为一长度为N的含高斯白噪声的单变量时间信号,选定一个时间延迟τ和嵌入维数m,构造如下的相空间矢量:

Xi=[xi,xi+τ,…,xi+(m-1)τ]T,i=1,2,…,M

(1)

式中:M=N-(m-1)τ。τ由平均互信息量法[21]确定,即

(2)

式中：yi=xi+τ。取I(τ)第一个极小值点对应的时间作为时间延迟。m由Cao[22]所提的方法确定,即

(3)

1.2 局部投影算法

基于相空间重构局部投影降噪算法,利用信号与噪声在相空间中轨线的动力学特性与几何特性的不同,保留信号分量,抑制噪声分量,最大程度地恢复原始信号的吸引子流形。

对于相点Xn,将其在邻域内线性化展开:

(4)

(5)

(6)

第n个相点Xn的邻域加权矩阵为

(7)

对式(7)进行SVD分解,得到:

(8)

式中:Vn=[a1,a2,…,aK]由信号子空间与噪声子空间构成。将SVD分解得到的奇异值按从小到大进行排列，有σ1≥σ2≥…≥σK-1≥σK,其中大的奇异值对应信号子空间,小的奇异值对应噪声子空间。根据Vn得到由噪声引起的分量,减去第n个相点在噪声子空间中的投影,即得到修正后的相点:

(9)

为了进一步提高降噪效果,避免局部线性化产生较大误差,本文采用Moore等[23]提出的质心修正方法。修正后的质心为

(10)

2 改进的局部投影算法

2.1 邻域选择

在局部投影降噪算法中,邻域的选择十分重要。邻域选择得过小,会损失有效信息,受噪声干扰严重;邻域选择得过大,会使得线性逼近的效果不好。

递归图[24]可以用来提取时间信号中的相关信息。但是，递归图的缺点是,动力系统的递归模式可视化对相似阈值的选取十分敏感[25]。Pham等[26]在递归图基础上提出了模糊递归图。相比于传统递归图,模糊递归图不需要选取相似阈值,且模糊递归图以灰度图的形式展示,能够为模式分析提供更丰富的信息。本文采用模糊递归图来对邻域进行选择。

(11)

式中：μ(Xi,Xj)∈[0,1]表示Xi和Xj之间的一种模糊相似性的度量,其具有如下性质。

(1) 自反性:

μ(Xi,Xi)=1

(12)

(2) 对称性:

μ(Xi,Vj)=μ(Vj,Xi)

(13)

(3) 传递性:

μ(Xi,Xj)=max[min{μ(Xi,Vk),μ(Xj,Vk)}]

(14)

式中：i=1,2,…,M;k=1,2,…，c；c为聚类数,1

μ(Xi,Xj)的值根据模糊C均值(fuzzy C-Mean, FCM)聚类算法[27]计算得到。FCM算法通过最小化如下模糊目标函数实现:

(15)

式中：ω∈[1,+∞)为模糊度参数;U=[μij](i=1,2,…,M;j=1,2,…,c)是划分矩阵;V=(V1,V2,…,Vc)是聚类中心矩阵;Vj表示第j个聚类中心;d(Xi,Vj)表示某一范数,本文采用欧式范数。上述模糊目标函数满足

(16)

为了得到最优的U和V,通过迭代过程数值求解目标函数的最小值,迭代过程如下:

(17)

(18)

当满足‖Ut-Ut+1‖<ε时,停止迭代。其中,t为迭代次数,ε为给定的精度水平。

模糊递归图具有对称性,可以看作是相空间状态矢量之间的一种模糊关系。模糊递归图用灰度图表示,灰度值代表了状态矢量对之间的模糊关系,这与递归图是相互兼容的。灰度值越小,则表示这两个状态矢量越相似。一个灰度值为0的像素点代表了两个状态矢量完全相似,即代表动力系统中一个100%的递归事件[28]。

2.2 信号与噪声子空间划分

在进行投影修正之前,需要根据SVD分解得到的奇异值对信号子空间与噪声子空间进行划分。本文利用Gavish等[29]提出的最优硬阈值准则对信号子空间与噪声子空间进行划分,不需要预先估计噪声水平。最优硬阈值γ为

γ=ω(β)σmed

(19)

式中:σmed为奇异值的中位数;ω(β)=λ(β)/μβ,λ(β)为

(20)

对于前述的邻域加权矩阵Zn∈Rm×n;当m=n时,β=1;当m>n时,β=n/m;当m

(21)

可以通过数值方法求解得到式(21)积分方程中μβ的值,从而求得阈值γ。

根据最优硬阈值准则,将大于等于阈值γ的奇异值所对应的奇异向量所形成的子空间作为信号子空间,小于阈值γ的奇异值所对应的奇异向量所张成的子空间作为噪声子空间,由此实现了子空间的划分。

2.3 改进算法的基本步骤

基于以上分析,本文所提的局部投影降噪算法步骤如算法1所示。

算法1 改进的局部投影降噪算法步骤输入:含噪信号x(t)输出:降噪信号x^(t)开始 相空间重构:利用式(2)计算平均互信息量并取第一个局部极小值对应的时间作为时间延迟τ;利用式(3)计算E1(m),取停止变化时对应的维数作为嵌入维数m;对x(t)进行相空间重构,得到式(1)表示的相空间X。 循环(i=1∶M) 1. 选定参考相点Xi。 2. 邻域选择:根据式(11)计算重构相空间的模糊递归图,得到第i个参考相点的邻域。 3. 计算邻域质心和邻域矩阵:由式(10)和式(7)分别计算第i个参考相点的邻域质心和邻域矩阵。 4. SVD分解:根据式(8)对邻域矩阵进行SVD分解,得到奇异值与右奇异向量。 5. 子空间划分:根据式(19)的最优硬阈值准则对信号子空间与噪声子空间进行划分。 6. 投影修正:根据式(9)进行投影修正。 结束 反重构:将相空间恢复为时间信号x^(t),恢复方式采用式(22)[30],以减小误差。结束

为减小反重构所产生的误差，进行如下操作：

(22)

为达到较好的效果，需将以上步骤重复迭代几次。

3 仿真实例与工程应用

3.1 Lorenz混沌信号仿真

本文首先采用Lorenz混沌系统信号进行仿真。Lorenz混沌系统由如下偏微分方程组描述：

(23)

当参数取σ=10, r=28, b=8/3时,系统呈现出混沌特性。采用四阶Runge-Kutta方法数值求解上述偏微分方程组,初始值取x0=10, y0=10, z0=10,积分步长为0.01。取X分量中的2 000个数据点进行仿真。对信号添加高斯白噪声,添加噪声后的信噪比为8 dB,然后使用本文所提的局部投影方法进行降噪处理。降噪效果如图1所示。

由图1可知,本文方法具有较好的降噪效果,在噪声较强的情况下也能够恢复信号。由图1可以看出，降噪后，信号比较光滑,看不到噪声的影响,同时基本保持了Lorenz混沌吸引子的几何形状。

3.2 参数对算法的影响

本小节讨论邻域选择和阈值这两个参数对算法的影响。首先是邻域选择对算法的影响。Lorenz信号添加噪声后的信噪比为8 dB。选取不同邻域对信号进行降噪,每个邻域仿真200次，然后取输出信噪比的平均值作为最终输出信噪比。由图2可知,当邻域选取合适时,降噪后信噪比达到最大。邻域选择过小或过大都会使得信噪比下降。这是因为邻域选择过小,信号受噪声影响明显;而邻域选择过大,对于相空间轨线的逐段线性逼近效果差。

图1 Lorenz混沌信号降噪前后时域波形图和相空间图Fig.1 Time-domain waveform and phase space diagram of Lorenz chaotic signal before and after denoising

图2 邻域选择对算法的影响Fig.2 Effect on the algorithm of neighborhood selection

仍然选取上述信号,在计算得到的阈值上叠加一个随机误差后再进行降噪处理。由图3可以看出,在有误差和无误差时，其降噪效果基本相同,这表明算法具有较好的鲁棒性。

图3 阈值误差对算法的影响Fig.3 Effect of threshold error on the algorithm

3.3 不同局部投影方法对比

为了进一步研究本文所提方法的降噪效果,将其他局部投影降噪算法与本文所提方法进行对比。选取降噪后的信噪比、均方误差、复杂度和耗时作为衡量降噪效果的评价指标。信噪比的计算公式为

(24)

(25)

首先,在信噪比为8 dB的情况下,计算混沌信号重要的特征量:李雅普诺夫指数、关联维数,然后再统计每个方法计算的耗时。每种方法设定相同参数,迭代8次,仿真200次,然后计算平均耗时。由表1可以看出,由本文方法降噪后的信号计算得到的李雅普诺夫指数和关联维数最接近原始Lorenz信号的对应值,这表明Lorenz混沌吸引子中的确定性结构得到了较好的保留。另外,由表1可以看出，本文方法较为耗时，这是因为本文方法在计算过程中涉及了模糊递归图的求解以及数值求解积分方程。但是本文方法在牺牲计算时间的情况下,具有更好的降噪效果。

表1 不同方法降噪后信号的李雅普诺夫指数、关联维数以及计算耗时(SNR=8 dB)

同样选定Lorenz混沌信号,添加0～20 dB的高斯白噪声,然后分别用这些方法进行降噪处理。为减小偶然误差,每种方法仿真200次,最后取200次的平均值作为最终的结果。

不同降噪方法的降噪后输出信噪比和均方误差的结果对比分别如图4和图5所示。基于小波的局部投影与基于EMD的局部投影的降噪效果接近,因为这两种方法均是采用预先估计噪声水平方式确定邻域的。平滑子空间局部投影通过平滑正交分解识别子空间,效果优于标准局部投影方法。递归局部投影采用递归图计算最优邻域半径,K均值局部投影利用K均值聚类方法选取邻域,能够实现更好的去噪效果。而本文方法同时考虑了邻域选取、子空间划分以及质心修正问题,相比上述方法,本文所提方法具有最高的输出信噪比和最低的均方误差,这说明方法较好地保留了信号时域波形的形状,具有较好的降噪性能。

图4 不同局部投影降噪方法的输出信噪比Fig.4 Output signal to noise ratio of different local projective denoising methods

图5 不同局部投影降噪方法的均方误差Fig.5 Mean square error of different local projective denoising methods

综上所述,本文所提的基于模糊递归图和最优硬阈值准则的局部投影方法能够有效滤除混沌信号中的噪声。

3.4 实测ECG信号仿真

ECG可以记录心脏的生理活动,在医疗临床领域应用广泛。但是，ECG信号在仪器测量过程中容易被噪声污染,从而影响各种生理特征的检测、提取和识别,容易造成误诊,因而ECG信号降噪对实际临床诊断具有重要的意义和价值。目前，常见的ECG信号去噪方法有奇异谱分析法、EMD方法、小波阈值方法等[31]。

ECG信号属于非平稳非线性信号,研究表明ECG信号具有复杂的动力学特性,表现出一定的混沌特性,因此非线性动力系统的方法作为一种研究心脏病的有效工具,被广泛应用于ECG信号的研究当中。

本文选取MIT-BIH(Massachusettes Institute of Technology-Boston’s Beth Israel Hospital)噪声压力测试数据库中编号为118e12的实测数据,其采样频率为f=360 Hz,采集到的信号含有基线漂移、噪声等影响。

首先，利用最小二乘拟合去掉信号内的基线漂移；然后，再使用本文方法、奇异谱分析方法、EMD方法和小波阈值方法分别对ECG信号进行处理及对比。处理结果如图6所示。

由图6对比去噪前后的相图,可以发现噪声大大减少，本文方法去噪后的相图表现出了较为清晰的吸引子结构,说明原来含有随机噪声的ECG信号在经过降噪处理后，其确定性成分得以加强,展现出了ECG信号原来的特征。对比本文方法,奇异谱分析(singular spectrum analysis, SSA)方法、EMD方法和小波阈值(wavelet threshold, WT)方法虽然也有较好的去噪效果,但时域波形还是不够光滑,而且相空间轨线细节部分呈现出杂乱堆叠的现象。

图6 ECG信号降噪前后波形图和相空间图Fig.6 Waveform and phase space diagram of ECG signal before and after denoising

为进一步说明降噪效果,计算如图7所示的原始信号和降噪信号在嵌入维数不同时的伪最近邻点比例。当伪最近邻点比例等于0或者小于某个值而保持不变,可认为不含伪最近邻点;当时间信号受到噪声影响时,混沌吸引子轨线受到影响,伪最近邻点的比例也会受到影响。由图7可知,降噪前伪最近邻点的比例在嵌入维数为4时降到较低,但没有降至0。降噪后,当嵌入维数为7时，伪最近邻点的比例就降至0。这进一步说明了应用本文所提方法对实测ECG信号可进行有效的降噪处理。

图7 ECG降噪前后的虚假最近邻点比例Fig.7 Proportion of false nearest neighbor points before and after ECG noise reduction

4 结 论

本文提出的基于模糊递归图和最优硬阈值的局部投影算法,主要针对原始局部投影算法中的邻域选择问题和子空间划分问题。在相空间重构的基础上,首先使用模糊递归图对邻域进行选择,同时为避免求取原始方法中的协方差矩阵,本文直接使用SVD分解,然后采用最优硬阈值准则对信号子空间和噪声子空间进行划分。通过对Lorenz混沌信号进行仿真,并与其他局部投影算法相比较,表明该方法具有一定的优越性。最后，将其应用于实测含噪的ECG信号,有效地降低了原始信号中的噪声,使信号的特征更加明显。同时，将本文所提算法与其他类型的ECG信号降噪方法进行了对比,证明了本文方法的有效性,为后续的信号处理奠定了基础。