CN104182617A - 基于内禀波形匹配的端点效应抑制方法 - Google Patents

Abstract

该方法公开了一种基于内禀波形匹配的端点效应抑制方法，属于信号分析、故障诊断领域，特别是基于波形匹配的端点效应抑制方法。该发明为对波形两端分别进行延拓，首先对接收波形进行消去中心趋势，然后选取末端的一段波形，与消除中心趋势的波形中其它段波形相匹配；如果匹配成功，则选取匹配成功的该段原始波形对接收波形进行延拓，若匹配不成功则用末端的一段波形进行延拓。从而具有准确性高、稳定性好的效果。

Description

基于内禀波形匹配的端点效应抑制方法

技术领域

本发明属于信号分析、故障诊断领域，特别是基于波形匹配的端点效应抑制方法。

背景技术

在IMF波形两端出现畸变或发散的现象叫做端点效应，又叫端点飞翼如图1所示。端点效应可能会导致以下问题：(1)端点效应会随着筛分过程不断向信号内部扩散，后续筛分出来的IMF端点效应会更严重；(2)端点效应使分解结果出现虚假模态；(3)端点效应使得端点附近的数据失去分析价值，舍去端点部分数据会降低分析资料的完整性。实际应用中，如果信号的数据序列比较长，通常的做法是不断抛弃信号两端的数据，使得受关注的信号段不受端点效应的影响。但是信号很短时，这种方法将完全行不通，此时就必须对端点效应进行抑制。

邵晨曦等提出了一种新的匹配算法，《一种自适应的EMD端点延拓方法》[J].电子学报,2007,35(10):1944-1948.邵晨曦,王剑,范金锋,et al，在计算匹配度时采用相对参考点相同长度的两段波形进行匹配。信号的延拓包括左右两端，下面以左端的延拓为例进行说明。设原始信号为s(t)。首先确定s(t)最左端的两个相邻极值点，不妨设其分别为极大值点和极小值点，分别记为P₀和P₁，从起始点到P₁的这段波形记为ω₀，设其长度为l；然后设E_max为s(t)的极大值点集合，以E_max-{P₀}中的每一个极大值点P_i作为参考点，计算该段相同长度的波形ω_i和ω₀的匹配度m(ω₀,ω_i,P_i)；接着记m(ω₀,ω_i,P_i)＝min{m(ω₀,ω_i,P_i),i＝1,2,…}，若m(ω₀,ω_i,P_i)＜α·l，其中α为一常数，则取ω_i0左侧包含了一个极大值和极小值的子波，作为原始s(t)左端的延拓，延拓完毕；最后直接指定端点处的极大和极小值：取原始信号最左端的两个相邻极大值点的均值作为左端点的极大值，取信号最左端的两个相邻极小值点的均值作为左端点的极小值，完毕。

邵晨曦提出的算法采用相对参考点的恒定长度来确定匹配子波的起始点，这可能导致错误的匹配结果。如图2所示，根据匹配误差计算公式这两段子波是相匹配的，但特征子波的右端点是极值，匹配子波的右端点却不是极值，显然这时将该匹配子波作为最佳匹配波是欠妥当的。

发明内容

本发明要解决的技术问题是研究设计一种准确性高、稳定性好的基于内禀波形匹配的端点效应抑制方法。

本发明的技术方案是基于内禀波形匹配的端点效应抑制方法，该方法包括对接收波形左右两端的延拓计算，首先对左端波形的延拓计算过程包括：

步骤1、接收到信号，设该信号波形中所有极大值点序列构成的集合为M＝{M₁,M₂,…,M_i,M_i+1,…,M_m}，设所有极大值点的相邻极小值点序列构成的集合为N＝{N₁,N₂,…,N_i,N_i+1,…,N_m}；设信号的起始点为S₁(x_s1,y_s1)，第一个极值点为极大值点，则第一段子波为s₁为以M₁为起点向右包含a个极大值点和a个极小值点的波段；则第i段子波起始点为S_i(x_si,y_si)，第i段子波为s_i为以M_i为起点向右包含a个极大值点和a个极小值点的波段；设极大值点的坐标为M_i＝(x_mi,y_mi)，极小值点的坐标为N_i＝(x_ni,y_ni)；根据相似波形的时间尺度比例，计算第i段子波的起始点S_i的坐标x_si，由于得到的x_si点可能不在采样点上，这时采用四舍五入取整；

步骤2、用曲线l_m1连接s₁的所有极大值点，用曲线l_n1连接s₁的所有极小值点，计算出l_m1、l_n1中心线l₁，则有用曲线l_mi连接s_i的所有极大值点，用曲线l_ni连接s_i的所有极小值点，计算出l_mi、l_ni的中心线l_i，则有

将s₁减去l₁得到一个消去中心趋势的特征子波，记为将s_i减去l_i得到一个消去中心趋势的匹配子波

步骤3、求取特征子波与各匹配子波的匹配度，找出其中最大的匹配度并将其与事先设定的阈值比较，若大于或等于阈值则匹配成功，转到步骤4，若匹配不成功则转到步骤5；

步骤4、在输入的信号波段中将最佳匹配子波的起始点S_i左边的包含k个极大值和k个极小值的波段作为延拓样本，消去该样本的中心趋势，结果记为S^α，将S^α加上步骤2中的l₁，如果l₁长度不够，则将l₁延长，得到延拓波形S；将波形S延拓到信号端点S₁左边；

步骤5、如果ρ＜β，则表示匹配不成功，此时直接用s₁进行镜像延拓，如果数据长度不够，则重复延拓该段波形，即完成了数据左端的内禀波形延拓；

步骤6、按照步骤1到步骤5的相同方法对输入波形右端进行延拓。

所述步骤1中a取值为2、3或4。

所述步骤1中第i段子波的起始点S_i的坐标x_si计算公式为：

$x_{\mathrm{si}}=\mathrm{round}\left(x_{\mathrm{mi}}\frac{(x_{m1}-x_{s1})(x_{m(i+1)}-x_{\mathrm{mi}})}{x_{m2}-x_{m1}}\right)$

其中，x_si是点S_i的坐标，round(·)表示四舍五入取整。

所述步骤3的具体步骤为：

步骤3.1、设特征子波长度为T₁，匹配子波长度为T_i，设较短子波的长度为T，即T＝min{T₁,T_i}；

步骤3.2、从较长子波中截取前T个数据计算归一化相关系数，并乘以长度比例因子作为匹配度，当T₁＞T_i时，有：

${\mathrm{\ρ}}_{1i}=\frac{{\left(S_1^{\mathrm{\α}}\right)}^T\·S_i^{\mathrm{\α}}}{{\left|\right|S_1^{\mathrm{\α}}\left|\right|}_2\·{\left|\right|S_i^{\mathrm{\α}}\left|\right|}_2}\×|\frac{T_i}{T_1}-1|=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{S}_1^{\mathrm{\α}}\left(t\right)S_i^{\mathrm{\α}}\left(t\right)}{{\left[\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{\left|S_1^{\mathrm{\α}}\left(t\right)\right|}^2\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{\left|S_i^{\mathrm{\α}}\left(t\right)\right|}^2\right]}^{1/2}}\×|\frac{T_i}{T_1}-1|$

式中，ρ_1i表示特征子波和第i段匹配子波的匹配度，||·||₂是求二范数运算；

步骤3.3、依次计算出与各的匹配度，ρ₁₂，ρ₁₃，…，ρ_1(m-2)，求出最大相关系数ρ＝ρ_1i＝max(ρ₁₂,ρ₁₃,…,ρ_1(m-2))，如果有多个相同值为最大值，则取离左端点最远的一个，给定判断阈值β，若ρ≥β，则匹配成功，将这段波形作为最佳匹配子波，转步骤4；否则转步骤5。

所述步骤3.3中β＝0.95。

所述步骤4中k＝5。

本发明的有益效果为

(1)如果EMD算法将端点值作为极值处理，数据序列延拓可以避免端点处插值点不正确而导致的错误；如果原方法在端点部分用外插值拟合包络线，数据序列延拓后，端点处的插值将由外插值变成内插值，避免了外插值的发散现象；

(2)数据序列延拓给出了端点外插值信息，消除了在端点处插值的方向错误；

(3)该方法延长了数据序列，分解时端点效应将释放到延拓出来的数据段上，原信号的端点效应将会减弱；

(4)波形延拓法能够最大限度的保留原信号的内在特征，有利于得到精度更高的分解结果。

附图说明

图1为端点效应现象示意图；

图2为邵晨曦算法匹配错误情形示意图；

图3为求取内禀波过程示意图；

图4为消去趋势项后的内禀波形图；

图5为匹配成功时的波形延拓示意图；

图6为各分量信号和原始信号示意图；

图7为采用不同端点效应抑制方法的抑制结果图；

图8为不进行端点抑制的分解结果示意图；

图9为采用本算法的分解结果示意图；

图10为具有内圈缺陷的滚动轴承的故障数据。

图7中(a)为镜像延拓法分解结果图；(b)为对称极值延拓法分解结果图；(c)为特征波法分解结果图；(d)为本征波匹配法分解结果图；(e)为邵晨曦等提出的匹配算法分解结果图；(f)为本算法分解结果图。

具体实施方式

本具体实施方式采用不同端点效应抑制方法对同一信号进行端点效应抑制，并用综合评价指标来评价不同方法的抑制效果。一般情况下，信号分析的要求是在兼顾分解效率的前提下，取得较好的分解精度；同时，抑制效果的好坏，主要体现在分解精度上。因此本发明在评价抑制效果时，将下式中各分量指标的系数比例设定为a₁:a₂:a₃:a₄＝1:3:2:4，则综合评价指标Index为：

$\begin{array}{c}\mathrm{Index}=0.1\×{\mathrm{CE}}^{*}+0.3\×{\mathrm{EL}}^{*}+0.2\×{\mathrm{\ρ}}_u^{*}+0.4\×{\mathrm{\Δ}}^{*}\\ =0.1\×{\left(\frac{t_{\mathrm{algorithm}}-t_{s\tan \mathrm{dard}}}{t_{s\tan \mathrm{dard}}}\right)}^{*}+0.3\×{|\frac{\underset{i=1}{\overset{n+1}{\mathrm{\Σ}}}{E}_i^{\′}}{\underset{i=1}{\overset{n+1}{\mathrm{\Σ}}}E}-1|}^{*}+0.2\×{(1-\frac{\underset{i=1}{\overset{n+1}{\mathrm{\Σ}}}\left|{\mathrm{\ρ}}_i\right|}{n+1})}^{*}+\\ 0.4\×{\left(\frac{\underset{i=1}{\overset{n+1}{\mathrm{\Σ}}}{\mathrm{\Δ}}_i}{n+1}\right)}^{*}\end{array}$

本发明具体实施方式采用的仿真信号由两个正弦波和一个调幅波叠加而成，解析表达式为：

$\left\{\begin{array}{c}{x}_1\left(t\right)=2\sin (2\mathrm{\π}\×25t)\\ x_2\left(t\right)=4\sin (2\mathrm{\π}\×10t)\sin (2\mathrm{\π}\×1/10t)\\ x_3\left(t\right)=\sin (2\mathrm{\π}\×5t)\\ x\left(t\right)=x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\end{array}\right.,t\∈\left[\mathrm{0,1}\right]$

仿真信号是一个非稳态信号，仿真信号及三个分量如图6所示。仿真信号不含调频分量，是因为调频波的特征尺度不断变化，分解时调频波的能量将泄漏到其他分量中，即出现模态混叠现象。模态混叠也会使分解误差增大，从而导致无法确定分解误差的诱发因素，因此本发明具体实施方式采用的仿真信号不含调频波分量。同时，由于端点效应对短信号影响较大，对该仿真信号在t∈[0,1]采样1000个点作为原始数据进行后续分析。

不抑制端点效应时，采用SD准则进行分解，设SD＝0.2320，分解结果各项基本指标如表1所示；采用G.Rilling准则进行分解，参数取默认值θ₁＝0.05，θ₂＝0.05，α＝0.5，分解结果各项基本指标如表1所示。从表1可知，采用G.Rilling准则时，分解结果的各项基本指标均优于采用SD准则，因此后续分解均采用G.Rilling准则。

采用不同端点效应抑制方法对仿真信号进行端点效应抑制，分解结果如图7所示；计算各项基本指标，结果如表2所示；对各基本指标进行一致化和归一化处理，并计算综合评价指标，结果表3所示。

对比表1和图7可知，各种抑制方法对端点效应都具有抑制作用，抑制端点效应后IMF对称性更好，能量泄漏也越少，分解结果更接近真实波形。从图7可知和表2可知：(1)不同端点效应抑制方法的抑制效果不同，并且本发明所提基本指标的评价结果与分解结果相一致，说明了本发明所提基本指标的正确性。(2)比较EL和IO指标，可知两个指标具有一定相关性，但并非严格正相关，这是由两个指标的评价侧重点不一样造成的。EL是从能量的角度度量虚假模态的多少，IO指标是从正交性的角度度量IMF的局部对称性。一般情况下，虚假分量越多，正交性越差，但也有例外，如图7(e)所示，虽然能量泄漏较多，但是IMF对称性很好，这是由于波形匹配不正确导致的。(3)镜像延拓法和对称极值延拓法的分解结果基本相同，各基本指标也基本相等，但对称极值延拓法的计算效率要高于镜像延拓法，因此对称极值延拓法的实用性更强。(4)本算法的能量泄漏、分解精度、正交性均优于其他算法，说明本算法的虚假分量更少，分解精度更高，IMF局部对称性更好，因此抑制效果更显著。

比较图7和表3中不同抑制方法的分解结果和对应的综合评价指标，可知抑制效果和综合评价指标的评价结果是一致的，从而证明了本发明所提综合评价指标的正确性。从表3可知，本算法的综合指标最好，在该综合指标评价结果中是最优的端点效应抑制方法。

下面通过对具有内圈缺陷的轴承故障信号的分解来验证本算法，该故障信号是一段具有内圈缺陷的滚动轴承的振动数据，采样频率为12000Hz，截取的故障数据长度为1024个数据点。如图10所示。

由于端点效应在短数据中表现更明显，因此从该故障数据中截取前512个数据点对各种端点效应抑制方法进行验证分析。当不抑制端点效应时，分解结果如图8所示；当采用本算法抑制端点效应后，分解结果如图9所示。比较图8和图9可知，采用本算法进行端点效应抑制后，EMD筛分得到的IMF数量更少，筛分速度更快。从图8可知，IMF1～IMF7的幅值不断减小，但IMF8、IMF9的幅值大于IMF7的幅值，这说明IMF8、IMF9都是虚假模态。从图9可知，IMF1～IMF5的幅值不断减小，因此抑制端点效应后的分解没有产生虚假模态。由此可知，本算法对端点效应的抑制是有效的，也说明了前面仿真结果的正确性。

表1  不进行端点效应抑制的指标

表2  不同抑制方法分解的基本指标

表3  不同抑制方法分解的基本指标(归一化和无量纲化后结果)

Claims (6)

1.基于内禀波形匹配的端点效应抑制方法，该方法包括对接收波形左右两端的延拓计算，首先对左端波形的延拓计算过程包括：

步骤1、接收到信号，设该信号波形中所有极大值点序列构成的集合为M＝{M₁,M₂,…,M_i,M_i+1,…,M_m}，设所有极大值点的相邻极小值点序列构成的集合为N＝{N₁,N₂,…,N_i,N_i+1,…,N_m}；设信号的起始点为S₁(x_s1,y_s1)，第一个极值点为极大值点，则第一段子波为s₁为以M₁为起点向右包含a个极大值点和a个极小值点的波段；则第i段子波起始点为S_i(x_si,y_si)，第i段子波为s_i为以M_i为起点向右包含a个极大值点和a个极小值点的波段；设极大值点的坐标为M_i＝(x_mi,y_mi)，极小值点的坐标为N_i＝(x_ni,y_ni)；根据相似波形的时间尺度比例，计算第i段子波的起始点S_i的坐标x_si，由于得到的x_si点可能不在采样点上，这时采用四舍五入取整；

步骤2、用曲线l_m1连接s₁的所有极大值点，用曲线l_n1连接s₁的所有极小值点，计算出l_m1、l_n1中心线l₁，则有用曲线l_mi连接s_i的所有极大值点，用曲线l_ni连接s_i的所有极小值点，计算出l_mi、l_ni的中心线l_i，则有

将s₁减去l₁得到一个消去中心趋势的特征子波，记为将s_i减去l_i得到一个消去中心趋势的匹配子波

步骤3、求取特征子波与各匹配子波的匹配度，找出其中最大的匹配度并将其与事先设定的阈值比较，若大于或等于阈值则匹配成功，转到步骤4，若匹配不成功则转到步骤5；

步骤4、在输入的信号波段中将最佳匹配子波的起始点S_i左边的包含k个极大值和k个极小值的波段作为延拓样本，消去该样本的中心趋势，结果记为S^α，将S^α加上步骤2中的l₁，如果l₁长度不够，则将l₁延长，得到延拓波形S；将波形S延拓到信号端点S₁左边；

步骤5、如果ρ＜β，则表示匹配不成功，此时直接用s₁进行镜像延拓，如果数据长度不够，则重复延拓该段波形，即完成了数据左端的内禀波形延拓；

步骤6、按照步骤1到步骤5的相同方法对输入波形右端进行延拓。

2.如权利要求1所述的基于内禀波形匹配的端点效应抑制方法，其特征在于步骤1中a取值为2、3或4。

3.如权利要求1所述的基于内禀波形匹配的端点效应抑制方法，其特征在于步骤1中第i段子波的起始点S_i的坐标x_si计算公式为：

$x_{\mathrm{si}}=\mathrm{round}\left(x_{\mathrm{mi}}\frac{(x_{m1}-x_{s1})(x_{m(i+1)}-x_{\mathrm{mi}})}{x_{m2}-x_{m1}}\right)$

其中，x_si是点S_i的坐标，round(·)表示四舍五入取整。

4.如权利要求1所述的基于内禀波形匹配的端点效应抑制方法，其特征在于步骤3的具体步骤为：

步骤3.1、设特征子波长度为T₁，匹配子波长度为T_i，设较短子波的长度为T，即T＝min{T₁,T_i}；

步骤3.2、从较长子波中截取前T个数据计算归一化相关系数，并乘以长度比例因子作为匹配度，当T₁＞T_i时，有：

${\mathrm{\ρ}}_{1i}=\frac{{\left(S_1^{\mathrm{\α}}\right)}^T\·S_i^{\mathrm{\α}}}{{\left|\right|S_1^{\mathrm{\α}}\left|\right|}_2\·{\left|\right|S_i^{\mathrm{\α}}\left|\right|}_2}\×|\frac{T_i}{T_1}-1|=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{S}_1^{\mathrm{\α}}\left(t\right)S_i^{\mathrm{\α}}\left(t\right)}{{\left[\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{\left|S_1^{\mathrm{\α}}\left(t\right)\right|}^2\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{\left|S_i^{\mathrm{\α}}\left(t\right)\right|}^2\right]}^{1/2}}\×|\frac{T_i}{T_1}-1|$

式中，ρ_1i表示特征子波和第i段匹配子波的匹配度，||·||₂是求二范数运算；

步骤3.3、依次计算出与各的匹配度，ρ₁₂，ρ₁₃，…，ρ_1(m-2)，求出最大相关系数ρ＝ρ_1i＝max(ρ₁₂,ρ₁₃,…,ρ_1(m-2))，如果有多个相同值为最大值，则取离左端点最远的一个，给定判断阈值β，若ρ≥β，则匹配成功，将这段波形作为最佳匹配子波，转步骤4；否则转步骤5。

5.如权利要求1所述的基于内禀波形匹配的端点效应抑制方法，其特征在于步骤3.3中β＝0.95。

6.如权利要求1所述的基于内禀波形匹配的端点效应抑制方法，其特征在于步骤4中k＝5。