基于局部投影与小波包分解的机电系统信号组合降噪方法与流程

本发明涉及流程工业生产系统监测数据信号降噪领域。，具体涉及一种基于局部投影与小波包分解的机电系统信号组合降噪方法。

背景技术：

流程工业生产系统生产设备繁多，且需要各种辅助系统，各结构单元之间不断地进行着物质、信息、能量的交换，系统内部关联耦合度高，是一个分布式的复杂机电系统。为监测和评估系统的实时状态，在生产系统中存在大量的监测点位，用dcs监测系统收集各监测点位的时间序列形成监测数据集。研究表明，流程工业系统监测数据具有混沌性，非线性等特点，而且系统工况波动剧烈，生产过程中存在着强烈的噪声和振动，对传感器的精度产生很大的影响。数据的不准确对后续基于数据的系统状态分析产生了较大的影响。而监测数据集的混沌性和强噪声等特点也给数据集降噪造成了很大困难。

针对混沌性数据的降噪问题，国内外专家做了很多相关的研究。farmer提出基于最大似然估计的方法有效地减小了混沌信号中的噪声；sun等使用了最小二乘支持向量机(ls-svm)类的方法进行了降噪；韩敏等对局部投影方法进行改进，提出了自适应的邻域选择方法，改善了局部投影方法的去噪效果。梁庆真等利用小波阈值方法与局部投影进行多次降噪，对心音信号起到了很好的降噪作用。其中基于局部投影的降噪方法，从系统嵌入学流形结构角度，进行局部线性多次逼近，从而有效抑制了混沌系统的噪声干扰。

技术实现要素：

本发明的目的在于提供基于局部投影与小波包分解的机电系统信号组合降噪方法，以克服现有技术的不足。

为达到上述目的，本发明采用如下技术方案：

基于局部投影与小波包分解的机电系统信号组合降噪方法，包括以下步骤：

步骤1)、根据takens嵌入定理，将时间序列经过重构得到一个m维相空间；

步骤2)、以m维相空间中的第n个相点即观测点y(n)作为参考点，计算每个点距离观测点y(n)的欧式距离，得到相点之间的欧式距离之后，根据距离寻找邻域相点满足dni＜ε；

步骤3)、根据得到的邻域相点，求取邻域相点的质心

步骤4)、在邻域中，计算构成的邻域矩阵an，构成m×m的协方差矩阵cn＝(an)^t(an)，计算cn的特征值s1≥s2≥…≥sm和对应的特征向量e1,e2,…,em最后得到q个最小特征值对应的特征向量并构成投影矩阵：

步骤5)、根据每一个邻域相点的邻域中心位置，对邻域相点进行修正，计算邻域相点的修正值；

步骤6)、将相空间进行反重构，恢复时间序列；

步骤7)、利用gp算法计算时间序列关联维度，根据关联维数的变化来确定最佳迭代次数，并输出降噪结果；若未达到最佳迭代次数，重新进行步骤1)至步骤6)，完成局部投影降噪；

步骤8)、在局部投影降噪结果基础之上，利用小波包方法进行降噪，抑制高频噪声的干扰。

进一步的，利用互信息法和cao方法计算时间序列最佳延迟时间和最佳嵌入维数，重构出时间序列的高维相空间。

进一步的，设一维含噪混沌时间序列由两部分组成：

y(n)＝x(n)+w(n),n＝1,2,...,n

式中：x(n)为无噪声混沌时间序列，w(n)为噪声时间序列，n为时间序列长度；根据takens嵌入定理，将y(n)经过重构得到一个m维相空间：

y(n)＝{y(n)-(m-1)τ,y(n)-(m-1)τ,...,y(n)}

式中:y(n)表示m维相空间中的第n个相点；m表示嵌入维数；τ表示延迟时间。

进一步的，以观测点y(n)为参考点，计算每个点距离该点的欧式距离dni：

式中，n,i＝1,2,…,n；得到相点之间的欧式距离之后，根据距离寻找邻域点满足dni＜ε。

进一步的，对相点邻域半径的确定，具体包括以下步骤：

(1)计算递归图rij；

rij＝θ(ε-||xi-xj||)i,j＝1,2,…,n

根据上式构造递归图，式中，ε为预先设置的邻域半径，θ是heaviside函数，当相点之间的距离小于ε时，令rij值为1，否则为0；

(2)计算递归图中受限对角线长度np(ε)：

计算递归图对角线长度hi，

定义一个门限值t，使其值等于hi的均值与三倍标准差之和。由于相邻的对角线不能同时很长，定义n*(ε)递归图中受限对角线长度，令np(ε)值为满足条件的对角线的个数；

(3)计算领域偏置系数β(ε)：

定义递归图中对角线的平均长度，其表达式为:

n*(ε)＝∑ijrij/(2n-1)

定义β(ε)如下式，

将ε从小到大进行搜索，当β(ε)取得极小值时，其对应的ε作为最佳邻域半径εopt。

进一步的，采用二阶方法计算邻域质心，通过将邻域中心向外移的方法使其与投影的超平面近似相切：

具体的，将邻域相点的平均值作为邻域中心，其公式如下：

式中，k为邻域点数，为邻域平局值，即为该领域的质心；

利用二阶多项式进一步地抑制噪声：

式中，为最终得到的对一阶邻域质心修正后的二阶邻域质心，为一阶邻域质心，μk是第k个领域，为第k个领域μk内各自的邻域质心。

进一步的，在邻域y(n)中，计算构成的邻域矩阵an，构成m×m的协方差矩阵cn＝(an)^t(an)，计算cn的特征值s1≥s2≥…≥sm和对应的特征向量e1,e2,…,em最后得到q个最小特征值对应的特征向量并构成投影矩阵：

q(n)＝∑q(req)(req^t)

r为q维对角阵，第1个和第m个元素为1,其余为1000。

进一步的，取降噪过程中的关联维数变化变缓处所对应的迭代次数作为最佳迭代次数，采用gp算法计算时间序列的关联维数。

与现有技术相比，本发明具有以下有益的技术效果：

本发明基于局部投影与小波包分解的机电系统信号组合降噪方法，首先利用互信息法和cao方法计算时间序列最佳延迟时间和最佳嵌入维数，重构出时间序列的高维相空间，由于强噪声的干扰作用，动力学系统的流形结构遭到严重破坏，因此先利用递归图方法计算合适的邻域半径，采用局部投影方法对相空间流形结构进行多次修复，并根据关联维数的变化来判别最佳的修复效果，修复后的相空间流形结构虽然得到改善，但由于高频噪声的干扰，其流形曲线并不是特别光滑，因此在局部投影降噪的基础上采用小波包降噪方法从频域进行降噪，抑制高频噪声干扰，得到较为光滑的流形曲线，取得更好的降噪效果。本方法同时兼顾了动力学系统流形结构的修复和频域的降噪，对于信噪比为5-25db的中高强度噪声有很好的降噪效果，本方法对于中高强度噪声的抑制优势，可以得到较好的发挥。

进一步的，采用二阶方法计算邻域质心，通过将邻域中心向外移的方法使其与投影的超平面近似相切，对邻域质心进行了估计，实现噪声的有效抑制。

进一步的，对邻域质心进行了估计，本质上是通过将邻域中心向外移的方法使其与投影的超平面近似相切。

附图说明

图1为关联维数与降噪后的snr值的归一化联合曲线示意图。

图2为原始lorenz序列的相空间示意图。

图3为不同方法降噪后相空间对比图；图3(a)、图3(b)、图3(c)、图3(d)分别为加噪声信号、小波包方法、局部投影方法以及本方法降噪后的信号相空间图。

图4为不同降噪方法对a_rki7650监测点位降噪效果。

图5为不同降噪方法对a_rki7650监测点位降噪后的相空间对比图，图5a、图5b、图5c、图5d分别为原始信号、局部投影降噪信号、小波包降噪信号与本文方法降噪信号。

图6为不同方法对5个监测点位降噪后的最大lyapunov指数对比图。

图7为不同方法对5个监测点位降噪后的snr值对比图。

图8为本发明方法流程图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

本发明一种基于局部投影与小波包分解的机电系统信号组合降噪方法，该方法先通过互信息法和cao方法计算时间序列最佳延迟时间和最佳嵌入维数，重构出时间序列的高维相空间，由于强噪声的干扰作用，动力学系统的流形结构遭到严重破坏，因此先利用递归图方法计算合适的邻域半径，采用局部投影方法对相空间流形结构进行多次修复，并根据关联维数的变化来判别最佳的修复效果，修复后的相空间流形结构虽然得到改善，但由于高频噪声的干扰，其流形曲线并不是特别光滑，因此本发明在局部投影降噪的基础上采用小波包降噪方法从频域进行降噪，抑制高频噪声干扰，得到较为光滑的流形曲线，取得更好的降噪效果。

如图8所示，基于局部投影与小波包分解的机电系统信号组合降噪方法，包括以下步骤：

步骤1)、根据takens嵌入定理，将时间序列经过重构得到一个m维相空间；

具体的，利用互信息法和cao方法计算时间序列最佳延迟时间和最佳嵌入维数，重构出时间序列的高维相空间：

设一维含噪混沌时间序列由两部分组成：

y(n)＝x(n)+w(n),n＝1,2,...,n

式中：x(n)为无噪声混沌时间序列，w(n)为噪声时间序列，n为时间序列长度；根据takens嵌入定理，将y(n)经过重构可以得到一个m维相空间。

y(n)＝{y(n)-(m-1)τ,y(n)-(m-1)τ,...,y(n)}

式中:y(n)表示m维相空间中的第n个相点；m表示嵌入维数；τ表示延迟时间。

步骤2)、以m维相空间中的第n个相点即观测点y(n)作为参考点，计算每个点距离观测点y(n)的欧式距离，得到相点之间的欧式距离之后，根据距离寻找邻域相点满足dni＜ε；

具体的，以观测点y(n)为参考点，计算每个点距离该点的欧式距离dni：

式中，n,i＝1,2,…,n；得到相点之间的欧式距离之后，根据距离寻找邻域相点满足dni＜ε；

对相点邻域半径的确定，具体包括以下步骤：

(1)计算递归图rij；

rij＝θ(ε-||xi-xj||)i,j＝1,2,…,n

根据上式构造递归图，式中，ε为预先设置的邻域半径，θ是heaviside函数，当相点之间的距离小于ε时，令rij值为1，否则为0；

(2)计算递归图中受限对角线的长度np(ε)

计算递归图中对角线(不包括主对角线)长度hi，

hi＝∑k-j＝irjk

定义一个阈值t，使其值等于hi的均值与其标准差之和。由于相邻的对角线不能同时很长，因此定义np(ε)为递归图中受限对角线的长度，且满足条件

(3)计算偏置系数δ(ε)

定义递归图中对角线的平均长度为n*(ε)，其表达式为

n*(ε)＝∑ijrij/n

一个较好的ε应使n*(ε)与np(ε)应尽可能的接近，因此，定义偏置系数β(ε)如下式，

将ε从小到大进行搜索，当δ(ε)取得极小值时，其对应的ε作为最佳邻域半径εopt。

步骤3)、根据得到的邻域相点求取邻域相点的质心；

邻域质心的计算采用了二阶方法，通过将邻域中心向外移的方法使其与投影的超平面近似相切，对邻域质心进行了估计实现噪声抑制。

具体的，最初cawley将邻域相点的平均值作为邻域中心，其公式如下：

式中，k为邻域点数，为邻域平局值，即为该领域的质心。这种方法由于局部线性化会产生较大的误差，各个邻域的中心与曲线都不是相切的，而向曲线内部发生了一定程度的偏移，对降噪效果影响较大；后来sauer为了进一步地抑制噪声而利用二阶多项式：

式中，为最终得到的对一阶邻域质心修正后的二阶邻域质心，为一阶邻域质心，μk是第k个领域，为第k个领域μk内各自的邻域质心。该方法通过将邻域中心向外移的方法使其与投影的超平面近似相切，获得邻域质心；

步骤4)、在邻域y(n)中，计算构成的邻域矩阵an，构成m×m的协方差矩阵cn＝(an)^t(an)，计算cn的特征值s1≥s2≥…≥sm和对应的特征向量e1,e2,…,em最后得到q个最小特征值对应的特征向量并构成投影矩阵：

具体的，在邻域y(n)中，计算构成的邻域矩阵an，构成m×m的协方差矩阵cn＝(an)^t(an)，计算cn的特征值s1≥s2≥…≥sm和对应的特征向量e1,e2,…,em最后得到q个最小特征值对应的特征向量并构成投影矩阵：

q(n)＝∑q(req)(req^t)

r为q维对角阵，第1个和第m个元素为1,其余为1000；其目的是减少边界误差带来的影响；

步骤5)、根据每一个相点的邻域中心位置，对相点进行修正，计算相点的修正值：

具体的，根据以下公式：

从而得到每个相点的修正值；

步骤6)、将相空间进行反重构，恢复时间序列：

步骤7)、利用gp算法计算时间序列关联维度，根据关联维数的变化来确定最佳迭代次数，并输出降噪结果；

根据流程工业系统监测数据的混沌性特点，纯净的混沌信号的关联维数一般比较小，当受到高维噪声干扰时，含噪声信号的关联维数变大，因此，在局部投影方法迭代降噪的过程中，降噪后的时间序列关联维数也会迅速变小，取降噪过程中的关联维数变化变缓处所对应的迭代次数作为最优次数，在此采用gp算法计算时间序列的关联维数；

若未达到最佳迭代次数，重新进行步骤1)至步骤6)；

步骤8)、在局部投影降噪结果基础之上，利用小波包方法进行降噪，抑制高频噪声的干扰：

在局部投影降噪的基础上采用小波包降噪方法从频域进行降噪，抑制高频噪声干扰，得到具有较为光滑流形的曲线，提升传统局部投影方法的降噪效果。

lorenz序列的仿真分析：

为了验证本文降噪方法的有效性，采用lorenz仿真序列进行降噪实验。分别采用小波包方法、局部投影方法和本文降噪方法对不同噪声强度的仿真序列进行降噪。并对降噪结果进行讨论分析。

lorenz方程的表达式为

对仿真序列的参数如下：σ＝19,r＝45.92,b＝4,设定起始点为[-1,0,1],积分步长h＝0.01,并将前8000个过渡点去掉，取后1500个点作为研究对象。对仿真序列施加14db的噪声干扰，使得加噪后的含噪声序列snr值约为14，然后分别利用小波包方法、局部投影方法和本文方法进行降噪处理。

在局部投影降噪过程中，随着迭代次数的增加，时间序列的关联维度会迅速变小，在此绘制关联维数与降噪后的snr值的归一化联合曲线如图1所示。由图1可知，当迭代次数为2时，snr值达到极大值，而此时的关联维数也从迅速下降变为了缓慢下降，在实际的降噪过程中，取关联维数变化变缓和处为最佳迭代次数。

图2为原始lorenz序列的相空间，图3(a)、图3(b)、图3(c)、图3(d)分别为加噪声信号、小波包方法、局部投影方法以及本方法降噪后的信号相空间图；

由图2和图3对比可以发现，加噪声后信号的相空间处于非常混乱的状态。经过小波包降噪之后，信号的高频噪声得到了很好的抑制，相空间流形曲线因此变得光滑，但是小波包方法对相空间的整体流形结构并没有明显的修复效果。

经过局部投影方法降噪之后，信号的相空间的整体流形结构得到了很好的修复，但相空间中的流形曲线并不光滑，这是由于局部投影方法对于信号的整体趋势的恢复能力较强，而对于高频噪声的抑制能力较弱所导致的。

而经过本文方法进行降噪之后，相空间的整体流形结构和流形曲线的光滑程度均得到了较好的修复。

对比三种方法降噪前后的snr值可知,当含噪声信号初始信噪比为14db时，小波包降噪方法与局部投影降噪方法降噪效果相差无几，而本文方法同时从相空间的流形结构和时间序列的频域进行降噪，降噪效果明显优于其他两种方法。

由于不同的降噪方法适合于不同的工况，当噪声强度不同时，降噪方法的降噪效果也不尽相同。为探究三种方法在不同噪声强度下的降噪效果，在此分别将信号加以不同强度的噪声，使含噪初始信号的snr值变化区间为5-25db时，对比不同的噪声强度下三种方法的降噪效果如图4所示。

图4中横坐标snr值为仿真序列的原始信噪比，信噪比越大代表信号中噪声强度越小，纵坐标snr增益为降噪后信号的snr与初始snr之间的差值，snr增益越大，表明该方法的降噪效果越好。

由图4可知，当噪声强度介于5-25db时，小波包的降噪效果明显劣于局部投影方法，这是由于当系统存在很强噪声时，系统的动力学流形结构遭到了严重破坏，仅从频域进行降噪很难达到很好的降噪效果。而局部投影方法直接针对系统的流形结构进行修复，对噪声起到了较好的抑制作用。而本文方法在中高强度噪声领域的降噪效果明显优于其它两种方法。本文方法同时兼顾了动力学系统流形结构的修复和频域的降噪，对于中高强度噪声有很好的降噪效果。

工业压缩机组监测序列的降噪分析：

流程工业系统工艺过程复杂，生产波动剧烈，因此dcs系统采集的数据集存在强噪声的特点。这在一定程度上影响了基于数据的系统状态分析的精度。由上述仿真分析可知，本文方法对于中高强度的噪声有更好的降噪效果。在此，选取压缩机组中5个变量进行降噪处理。变量信息如下表所示。

上述监测序列分别用局部投影方法，小波包方法和本文方法三种方法对5个监测变量进行降噪处理。为对比降噪结果，本文采用最大lyapunov指数和snr值作为评价指标。

流程工业系统监测数据具有混沌性，其最大lyapunov指数一般大于0，并且而由于具有高维动力学特性的噪声干扰，系统动力学流形轨迹覆盖整个轨迹区间的速度会加快，导致最大lyapunov指数变大，而经过降噪后的时间序列其最大lyapunov指数都会变小。因此，在此认为最大lyapunov指数越小，降噪效果越好。

原始信号为含噪声信号，将降噪后的信号作为不含噪声的原始信号，计算降噪后的snr值，snr值越小，则认为降噪程度越大，降噪效果越好。

不同方法对5个监测变量降噪后的最大lyapunov指数和snr值如图6、图7所示。对比不同方法的降噪结果可知，经本文方法降噪后，各监测变量的最大lyapunov指数和snr值均比其他两种方法更小，说明本文方法对实际监测数据具有更好的降噪效果。

为进一步对比三种降噪方法的降噪效果，本文列出三种不同降噪方法对a_rki7650(汽轮机键相)监测点位降噪后的相空间轨迹，如图5所示。

图5a、图5b、图5c、图5d分别为原始信号、局部投影降噪信号、小波包降噪信号与本文方法降噪信号；由图5可知，原始信号由于存在生产过程中的强噪声干扰，相空间轨迹杂乱无章，看不出明显的特征；经过局部投影降噪后，相空间有了比较明显的质心，结合上述仿真分析的结论，局部投影方法在一定程度上修复了系统的动力学流形结构。而经过小波包降噪处理后，时间序列中的高频噪声得到了很好的抑制，因此相空间轨迹变得十分光滑，但系统在动力学上的整体流形结构并没有发生改变。而经过本文方法降噪处理过后，在去除高频噪声干扰的同时，对系统的动力学流形结构也进行了很好的修复，降噪效果优于其他两种方法。

本文针对于流程工业系统监测数据的混沌性和强噪声的特点，同时从动力学系统流形结构角度和频域进行降噪。首先利用互信息法和cao方法计算时间序列最佳延迟时间和最佳嵌入维数，重构出时间序列的高维相空间。由于强噪声的干扰作用，动力学系统的流形结构遭到严重破坏，因此本发明先利用递归图方法计算合适的邻域半径，采用局部投影方法对相空间流形结构进行多次修复，并根据关联维数的变化来判别最佳的修复效果。修复后的相空间流形结构虽然得到改善，但由于高频噪声的干扰，其流形曲线并不是特别光滑。因此本发明在局部投影降噪的基础上采用小波包降噪方法从频域进行降噪，抑制高频噪声干扰，得到较为光滑的流形曲线，取得更好的降噪效果。

对比小波包降噪，局部投影方法降噪与本发明降噪方法的降噪效果如下表，并结合图1降噪前后相空间轨迹可知：当含噪信号的初始信噪比为14db时，小波包降噪方法与局部投影降噪方法降噪效果相差无几；而本文方法同时从相空间的流形结构和时间序列的频域进行降噪，降噪效果明显优于其他两种方法。

表1.不同降噪方法降噪前后snr值对比表

由于不同的降噪方法适合于不同的工况，当噪声强度不同时，降噪方法的降噪效果也不尽相同。

当噪声强度逐渐减小，噪声对系统流形结构的破坏程度越来越小，从频域进行降噪的小波包方法优势越来越明显。当初始snr值大于13db时，小波包方法降噪效果明显提升，大幅超过局部投影方法。

而本文方法在中高强度噪声领域的降噪效果明显优于其它两种方法。本文方法同时兼顾了动力学系统流形结构的修复和频域的降噪，对于信噪比为5-25db的中高强度噪声有很好的降噪效果。而当噪声强度非常小时，本方法降噪效果与小波包方法相当，甚至弱与小波包降噪方法。而流程工业系统监测序列中具有很强的噪声，本方法对于中高强度噪声的抑制优势，可以得到较好的发挥。