一种多频率成分奇异值分解降噪方法

一种多频率成分奇异值分解降噪方法

基于变异值的确定方法信号的消噪一直是信号诊断和故障诊断的非常重要步骤。奇异值分解(singularvaluedecomposition,简称SVD)作为一种数据处理方法已经被成功地运用到信号消噪处理中并且被证明是有效的。该方法有两个关键:a.如何确定分解后重构的有效秩阶次;b.如何确定重构矩阵的行列数。针对有效秩阶次的选择,常用的方法是试凑法和阈值法,均依赖于经验,缺乏依据。朱启兵提出了基于结构风险最小化原则的奇异值分解降噪方法。王维提出了基于非监督动态聚类算法来确定矩阵有效重构阶次的方法。康春玉采用主分量分析的方法,根据奇异值的大小来确定有效秩的阶次。孙鑫晖提出了通过奇异熵增量确定降噪阶次的方法。针对重构矩阵的行列数的选择问题,普遍采取的方法是根据具体信号选择不同行列数进行试凑,这种试凑法需要大量的计算,并且严重依赖使用者的信号分析经验。Kanjilal提出了通过奇异值比谱来确定行列数的方法,但是该方法在信号由多个周期分量组成或噪声较严重的情况下效果并不明显。赵学智提出通过分析所有分解分量信号所含信息量的变化趋势来确定合理的矩阵结构。上述方法在实际应用中取得了较好的效果,但也存在着一定的局限性。本文提出了一种根据噪声信号的快速傅里叶变换结果来决定有效秩阶次,以降噪信号的信噪比和均方差大小为依据确定重构矩阵结构的SVD方法。1svd通信的基本理论给定一个秩为r的m×n维矩阵X,则X的奇异值分解形式为X=U[Σ000]VΗX=U[Σ000]VH(1)角线元素矩阵U,V分别为m×m,n×n维正交矩阵;Σ为r×r维对角阵,其对角线元素为矩阵X的非零奇异值σi,且以非增顺序排列,即σ1≥σ2≥…≥σr;0为零元素矩阵。矩阵X的秩为r,从式(1)中除去X的零奇异值,得到X奇异值分解的精简形式X=r∑i=1σiuivΗi(2)相空间重构矩阵对于一个测得的信号y(i)(i=1,2,3,…,N),基于相空间重构理论,可以由其构造重构矩阵AA=[y(1)y(2)⋯y(Ν-L+1)y(2)y(3)⋯y(Ν-L+2)⋯⋯⋯⋯y(L)y(L+1)⋯y(Ν)](3)i估计和重构矩阵估计若源信号y(i)是由有用信号和噪声共同组成,则矩阵A也是由有用信号和噪声共同组成的矩阵,那么矩阵A的奇异值σi可以反映信号和噪声能量集中的情况。前p个较大的奇异值将主要反映有用信号,较小的奇异值则主要反映噪声,把这部分反映噪声的奇异值置零就可以去除信号中的噪声,利用式(2)进行重构矩阵的估计,将矩阵中相应的项相加,取平均值就可以还原出信号。利用上述方法对信号进行降噪,要达到较好的效果,关键是确定有效秩的阶次和重构矩阵的结构。2mse与snr的定义降噪效果一般用信号的均方误差(meansquareerror,简称MSE)和信噪比(signaltonoiseratio,简称SNR)来衡量:MSE越小,SNR越大,降噪效果越好。MSE与SNR的定义形式如下ΜSE=1ΝΝ∑k=1(x(k)-(ˆx(k))2(4)SΝR=10log[Ν∑k=1x2(k)Ν∑k=1(x(k)-(ˆx(k))2](5)其中:x(k)为无噪声信号的第k个数据点;ˆx(k)为含噪声信号的第k个数据点;N为信号长度。2.1源信号的主导个数由第i个非零奇异值σi重构得到重构信号分量yi=σiuivHi,分别对这些分量信号进行傅里叶变换后发现,分量的频率成分均为源信号的频率成分组成,而由较大的奇异值重构得到的分量信号其频率成分与源信号中主频率相对应。显然,当有用信号未被噪声完全淹没时,源信号中主频率是有用信号的频率。因此可以通过奇异值与有用信号频率之间的某种对应关系来确定有效秩的阶次,从而达到降噪的效果。给定一个混有噪声的信号s=s0+ξ(n),其中:s0=2cos(0.4πt)+sin(0.2πt)为有用信号;ξ(n)为强度为1的高斯白噪声。该信号快速傅里叶变换(FFT)的结果如图1所示。它包含两个主频率成分,f1=0.1Hz,f2=0.2Hz,这是有用信号的频率成分,其他均为噪声的频率成分。分别对有用信号和源信号s利用式(3)构造相应的重构矩阵A,并对其进行奇异值分解。随着矩阵A的行数变化,其奇异值的分布如图2所示。可以看到,当行数小于主频个数时,无论是有用信号还是源信号,其非零奇异值的个数等于重构矩阵行数。当行数大于主频个数时,有用信号的非零奇异值有4个,不随重构行数的变化而变化;源信号的前4个奇异值也呈现相同的变化规律,主要反映有用信号的信息,将这些奇异值称为大奇异值。可以看到,随着行数的增加,大奇异值的个数恒为4,是源信号中主频个数的2倍,而其他奇异值相对较小且分布比较集中,反映出噪声的特点。用前4个奇异值进行重构得到降噪信号s′,其波形与源信号以及有用信号波形的对比如图3所示。从图3(b)可以看出,s′与s0几乎重合,说明前4个奇异值很好地重构了有用信号,抑制了噪声。i取不同的值,用前i个奇异值进行信号重构,比较重构信号的SNR与MSE值发现,取前4个奇异值进行重构,降噪效果最好。因此可以确定该信号有效秩的阶次为4。对不同主频个数的信号进行处理和分析,得到相同的结论:随着重构矩阵行数的增加,大奇异值的个数是源信号中主频率成分个数的2倍,并且恒定不变。因此,有效秩的阶次与源信号的主频个数存在一个确定的关系,即p=2n,其中:p为有效秩的阶次;n为源信号快速傅里叶变换后主频的个数。利用这种倍数关系可以确定有效秩的阶次。2.2噪声重构下的噪声误差分析确定有效秩的阶次之后,还需要确定重构矩阵的结构,即矩阵的行列数。本文根据降噪效果的好坏来确定最佳的矩阵结构。选取不同行数L(L>4)对信号s(N=800)进行奇异值分解与重构,降噪信号的SNR与MSE随L变化如图4所示。可以看出,当L增加到一定程度时,降噪信号的SNR基本稳定,仅有小的波动,此时MSE的变化呈现出这种规律。由计算结果可知,当L=440时,降噪信号的MSE取得最大值,SNR取得最小值,此时降噪效果最好。对不同的数据长度N以及不同频率组成的源信号进行分析后发现,SNR与MSE分别取最大和最小值时的L值不一定相等,但最佳L值基本出现在L=N/2处的一个领域内,并且L在该邻域内取值时,降噪效果较好且差异较小,均能满足要求。因此,重构矩阵的结构可以根据N来确定,工程应用中不妨取L=N/2(当N不是偶数时,舍弃最后一个数据点,不影响最终结果)。基于上述分析,对于一个含噪声的测试信号,其降噪的基本步骤如下:(1)取信号数据长度的一半作为重构矩阵的行数,根据式(3)构造重构矩阵并进行奇异值分解;(2)对信号进行快速傅里叶变换,确定主频个数n,以2n作为有效秩的阶次;(3)用前2n个奇异值根据式(2)进行重构,得到重构矩阵A2n,将A2n中对应的元素相加后平均得到降噪后的信号。3噪声信号的svd噪声仿真分别用不同频率成分的信号对该方法进行验证,结果如图5所示。3个源信号分别为s1=sin(0.5πt)+ξ(n)s2=5sin(0.5πt)+7sin(0.012πt)cos(0.073πt+0.275π)+ξ(n)s3=cos(0.1πt)+sin(0.3πt)+cos(0.5πt)+sin(0.7πt)+ξ(n)其中:ξ(n)为强度分别为5,18和8dB的高斯白噪声。从快速傅里叶变换结果可以看到,s1,s2,s33个信号的主频数分别为1,3,4,有效秩的阶次分别为2,6,8。3个信号数据长度均为600,取行数L=300。信号降噪前后的SNR与MSE值如表1所示。图5中,降噪信号的波形接近于有用信号的波形。从表1的结果可以看出,降噪后信号的SNR有很大提高,MSE明显降低,说明本文给出的降噪方法是有效的。含白噪声信号的SVD降噪研究表明,该方法对含色噪声的信号有较好地降噪效果。对仿真信号s0=2cos(0.4πt)+sin(0.2πt)添加色噪声η(n)得到源信号s4=s0+η(n)。η(n)为由强度为12dB的高斯白噪声ξ1(n)和服从均匀分布、幅值为1mV的白噪声ξ2(n)通过下面的非线性组合得到η(n)=ξ1(n)+ξ2(n)1+ξ12(n)(6)用本文的降噪方法对s4进行降噪处理,结果如图6所示。降噪后,信号的SNR由3.4496dB提高到21.6103dB,而MSE由9.3354下降到0.1163,说明降噪效果非常明显。从图6可以看到,降噪信号与无噪声信号s0的波形吻合较好。通过仿真计算发现,随着噪声强度的增加,信号的信噪比降低,有用信号逐渐被噪声淹没。当有用信号与噪声信号的频率分布比较均衡时,在傅里叶变换结果中无法区分主频与噪声频率,大奇异值包含的有用信号信息和噪声信息也比较均衡,无法利用大奇异值来抑制噪声,SVD降噪法失效。4l=n/2(1)有效秩的阶