专利摘要

本发明公开了一种基于经验模式分解滤波的包络分析方法，该方法首先利用经验模式分解方法对原始信号进行分解，然后利用数据的重排和替代操作排除分解结果中的噪声分量和趋势项，接着再采用谱峭度方法对第一次滤波后的信号进行分析，得到最优滤波器的中心频率和带宽，然后利用该滤波器对第一次滤波后的信号再进行第二次滤波，然后采用有理样条迭代平滑包络分析方法对第二次滤波后的信号进行包络分析，最后根据包络谱确定旋转机械的故障类型。本发明适合于处理复杂的旋转机械故障信号，能够准确地判定出旋转机械的故障类型，具有良好的抗噪性和鲁棒性，便于工程应用。

权利要求

1.一种基于经验模式分解滤波的包络分析方法，其特征在于，包括以下步骤：

步骤1：利用加速度传感器以采样频率fs测取旋转机械的振动信号x（k）, k=1, 2, …,N，N为采样信号的长度；

步骤2：采用经验模式分解（Empirical Mode Decomposition, EMD）算法将振动信号x（k）分解成n个分量和一个趋势项之和，即，其中，ci（k）代表由EMD算法得到的第i个分量，rn（k）代表由EMD算法得到的趋势项；

步骤3：对ci（k）执行重排操作和替代操作，经重排操作得到的数据用cishuffle（k）表示，替代操作后得到数据用ciFTran（k）表示；

步骤4：对ci（k）、cishuffle（k）和ciFTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，ci（k）的广义Hurst指数曲线用Hi（q）表示；

步骤5：如果Hi（q） 与Hishuffle（q）或Hi（q） 与HiFTran（q）之间的相对误差小于5%，或者Hi（q） 、Hishuffle（q） 和HiFTran（q）三者都不随q而变化，则抛弃对应的ci（k）分量；

步骤6：对剩余的ci（k）分量求和，将该和记为信号经重排和替代滤波后的结果xf1（k）；

步骤7：对xf1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f0和带宽B；

步骤8： 根据中心频率f0和带宽B对xf1（k）进行带通滤波，得到xf2（k）；

步骤9：对信号xf2（k）执行有理样条迭代平滑包络分析，得到信号包络eov（k）；

步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

2.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于，所述步骤2中经验模式分解（EMD）算法包括以下步骤：

1） 第一个筛选过程：分别找出振动信号x（k）的上、下局部极值点，采用三次样条曲线分别拟合上、下局部极值点，得到振动信号x（k）的局部最大值包络线和局部最小值包络线，再将此两条包络线的相应各点的值取平均，得到一条平均曲线m1；

再求振动信号x（k）与此平均曲线m1的差，即h10=x(k)-m1，至此第一个筛选过程结束；

2）第二个筛选过程：h10重新被当作振动信号，重复上述步骤1），可得到h11= h10-m11，这里参数m11代表h10的均值曲线，重复这一过程j次，直到0.2<SD<0.3时筛选过程停止，这里， 此时，h1j= h1(j-1)-m1j，这时认为h1j是一个内秉模态函数（Intrinsic Mode Function, IMF），定义第1个IMF为c1=h1j；

3） 从振动信号中减去c1，可得r1=x(k)-c1，再将r1当作新的数据，并重复上述两步操作，这样得到第2个IMF；

4）重复步骤3）操作可得到一系列IMF，如果rn已经变成一条单调曲线，则筛选过程停止，最终将振动信号分解为如下形式： 。

3.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于，所述步骤3中数据重排操作包括以下步骤：

随机打乱分量ci（k）的排列顺序。

4.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于：所述步骤3中数据替代操作包括以下步骤：

1） 对分量ci（k）执行离散傅里叶变换，获得分量ci（k）的相位；

2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量ci（k）的原始相位；

3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据ciIFFT（k），求取数据ciIFFT（k）的实部。

5.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于：所述步骤4中MFDFA方法包括以下步骤：

1）构造xy(k)的轮廓Y(i)，k=1,2,…,N：

,

，

xy(k)代表权利要求1所述步骤4中的ci（k）或cishuffle（k）或ciFTran（k）；

2）将信号轮廓Y(i)分成不重叠的Ns段长度为s的数据，由于数据长度N通常不能整除s，所以会剩余一段数据不能利用；

为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2Ns段数据；

3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：

4P`F)18O2626ZX3J9SV283R

Y$[`H}@AV[I7X%IJ6A_JQWC；

yv(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；本例中，m=1；

4） 计算第q阶波动函数的平均值：

0R(8CJGG2TXCM$2NT2@((TH；

5）如果xy(k)存在自相似特征，则第q阶波动函数的平均值Fq(s)和时间尺度s之间存在幂律关系：

Fq(s)~sH(q)；

当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：

S42M8HF%L8]B11YUNN21RY7；

6）对步骤5）中的公式两边取对数可得ln[Fq(s)]=H(q)ln(s)+c，c为常数，由此获得直线的斜率H(q)。

6.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于：所述步骤7中的谱峭度方法包括以下步骤：

1）构造一个截止频率为fc=0.125+ε的低通滤波器h(n)；ε>0；

2）基于h(n)构造通频带为[0;0.25]的准低通滤波器h0(n)和通频带为[0.25;0.5]的准高通滤波器h1(n)，

；

3）信号 经h0(n)、h1(n)滤波并降采样后分解成低频部分和高频部分，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2k个频带，其中表示滤波器树中第k层上的第i个滤波器的输出信号，i=0,…,2k-1，0≤k≤K-1；

4）分解树中第k层上的第i个滤波器的中心频率fki和带宽Bk分别为

fki=(i+2-1)2-k-1

Bk=2-k-1；

5）计算每一个滤波器结果 i=0,…,2k-1的峭度；

6）将所有的谱峭度汇总，得到信号总的谱峭度。

7.根据权利要求1所述的一种基于经验模式分解滤波的包络分析方法，其特征在于，所述步骤9中的有理样条迭代平滑包络分析方法包括以下步骤：

1）计算信号z(k)的绝对值∣z(k)∣的局部极值；在第1次迭代中，z(k)代表权利要求1所述步骤9中xf2（k）；

2）采用有理样条曲线拟合局部极值点得到包络线eov1（k）；

3）对z(k)进行归一化处理得到；

4）第2次迭代：把z1(k)重新作为新数据，重复执行上述步骤1）~3），得到；

5）第i次迭代：把zi-1(k) 重新作为新数据，重复执行上述步骤1）~3)，得到；

6) 如果第n次迭代得到的zn(k)的幅值小于或等于1，则迭代过程停止，最后得到信号z(k)的包络为。

说明书

技术领域

本发明涉及旋转机械状态监测与故障诊断领域，具体涉及一种基于经验模式分解滤波的包络分析方法。

背景技术

包络分析技术广泛应用于齿轮和滚动轴承的故障诊断中。现有的包络分析技术有下面三个缺陷：①现有的包络分析技术或者是直接对原始信号进行分析，或者是仅对原始信号进行简单的滤波后再进行分析，因此现有的方法容易受到噪声、趋势及其它成分的干扰，从而导致现有技术的分析精度较低；②现有的包络分析技术是以Hilbert变换为基础，而Hilbert变换要求被分析的信号必须是单分量的窄带信号，否则信号的频率调制部分将要污染信号的幅值包络分析结果，但是目前待分析的信号都不严格满足单分量且窄带的条件，这样就会导致现有技术因精度不高而容易出现误判问题；③由传统方法得到的包络谱存在端点效应。

发明内容

本发明要解决的问题是针对以上不足，提出一种基于经验模式分解滤波的包络分析方法，采用本发明的包络分析方法后，具有分析结果准确度和精确度高，并能准确地检测出旋转机械故障类型的优点。

为解决以上技术问题，本发明采取的技术方案如下：一种基于经验模式分解滤波的包络分析方法，其特征在于，包括以下步骤：

步骤1：利用加速度传感器以采样频率fs测取旋转机械的振动信号x（k）, （k=1,2, …,N），N为采样信号的长度；

步骤2：采用经验模式分解（Empirical Mode Decomposition, EMD）算法将信号x（k）分解成n个分量和一个趋势项之和，即 ，其中，ci（k）代表由EMD算法得到的第i个分量，rn（k）代表由EMD算法得到的趋势项；

步骤3：对ci（k）执行重排操作和替代操作，经重排操作得到的数据用cishuffle（k）表示，替代操作后得到数据用ciFTran（k）表示；

步骤4：对ci（k）、cishuffle（k）和ciFTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，ci（k）的广义Hurst指数曲线用Hi（q）表示；

步骤5：如果Hi（q） 与Hishuffle（q）或Hi（q） 与HiFTran（q）之间的相对误差小于5%，或者Hi（q） 、Hishuffle（q） 和HiFTran（q）三者都不随q而变化，则抛弃对应的ci（k）分量；

步骤6：对剩余的ci（k）分量求和，将该和记为信号经重排和替代滤波后的结果xf1（k）；

步骤7：对xf1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f0和带宽B；

步骤8： 根据中心频率f0和带宽B对xf1（k）进行带通滤波，得到xf2（k）；

步骤9：对信号xf2（k）执行有理样条迭代平滑包络分析，得到信号包络eov（k）；

步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

一种优化方案，所述步骤2中经验模式分解（EMD）算法包括以下步骤：

1） 第一个筛选过程：分别找出数据x（k）的上、下局部极值点，采用三次样条曲线分别拟合上、下局部极值点，得到信号x（k）的局部最大值包络线和局部最小值包络线，再将此两条包络线的相应各点的值取平均，得到一条平均曲线m1；

再求原始信号x（k）与此平均曲线m1的差，即h10=x(k)-m1，至此第一个筛选过程结束；

2）第二个筛选过程：h10重新被当作原始信号，重复上述步骤1），可得到h11= h10-m11，这里参数m11代表h10的均值曲线，重复这一过程j次，直到0.2<SD<0.3时筛选过程停止，这里 ， 此时，h1j= h1(j-1)-m1j，这时可以认为h1j是一个内秉模态函数（Intrinsic Mode Function, IMF），定义第1个IMF为c1=h1j；

3） 从原信号中减去c1，可得r1=x(k)-c1，再将r1当作新的数据，并重复上述两步操作，这样可以得到第2个IMF；

4）重复步骤3）操作可得到一系列IMF，如果rn已经变成一条单调曲线，则筛选过程停止，最终将原始信号分解为如下形式： 。

进一步地，所述步骤3中数据重排操作包括以下步骤：

随机打乱分量ci（k）的排列顺序。

进一步地，所述步骤3中数据替代操作包括以下步骤：

1） 对分量ci（k）执行离散傅里叶变换，获得分量ci（k）的相位；

2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量ci（k）的原始相位；

3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据ciIFFT（k），求取数据ciIFFT（k）的实部。

进一步地，所述步骤4中MFDFA方法包括以下步骤：

1）构造x(k)(k=1,2,…,N)的轮廓Y(i)：

；

x(k)代表权利要求1所述步骤4中的ci（k）或cishuffle（k）或ciFTran（k）；

2）将信号轮廓Y(i)分成不重叠的Ns段长度为s的数据，从数据的反方向以相同的长度分段，得到2Ns段数据；

3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：

；

yv(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；

4） 计算第q阶波动函数的平均值：

；

5）如果x(k)存在自相似特征，则第q阶波动函数的平均值Fq(s)和时间尺度s之间存在幂律关系：

Fq(s)~sH(q)；

当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：

；

6）对步骤5）中的公式两边取对数可得ln[Fq(s)]=H(q)ln(s)+c (c为常数)，由此可以获得直线的斜率H(q)。

进一步地，所述步骤7中的谱峭度方法包括以下步骤：

1）构造一个截止频率为fc=0.125+ε的低通滤波器h(n)；ε>0；

2）基于h(n)构造通频带为[0;0.25]的准低通滤波器h0(n)和通频带为[0.25;0.5]的准高通滤波器h1(n)，

；

3）信号 经h0(n)、h1(n)滤波并降采样后分解成低频部分 和高频部分 ，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2k个频带，其中 表示滤波器树中第k层上的第i个滤波器的输出信号，i=0,…,2k-1，0≤k≤K-1；

4）分解树中第k层上的第i个滤波器的中心频率fki和带宽Bk分别为

fki=(i+2-1)2-k-1

Bk=2-k-1；

5）计算每一个滤波器结果 ( i=0,…,2k-1)的峭度 ；

6）将所有的谱峭度汇总，得到信号总的谱峭度。

进一步地，所述步骤9中的有理样条迭代平滑包络分析方法包括以下步骤：

1）计算信号z(k)的绝对值∣z(k)∣的局部极值；在第1次迭代中，z(k)代表权利要求1所述步骤9中xf2（k）；

2）采用有理样条曲线拟合局部极值点得到包络线eov1（k）；

3）对z(k)进行归一化处理得到 ；

4）第2次迭代：把z1(k)重新作为新数据，重复执行上述步骤1）~3），得到 ；

5）第i次迭代：把zi-1(k) 重新作为新数据，重复执行上述步骤1）~3)，得到 ；

6) 如果第n次迭代得到的zn(k)的幅值小于或等于1，则迭代过程停止，最后得到信号z(k)的包络为 。

本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：

1)利用EMD对原始信号进行分解，然后利用数据的重排和替代操作排除其中的噪声和趋势分量，仅仅保留信号分量中的有用成分，从而避免了噪声和趋势分量对包络分析结果的影响，分析结果准确度和精确度高。

2)利用有理样条迭代平滑包络分析方法将信号包络与频率调制部分完全分离，能够避免频率调制部分对信号包络分析结果的影响，从而提高包络分析的精度。

3) 能够准确地检测出旋转机械的故障类型。

4)得到的包络谱能够避免端点效应。

下面结合附图和实施例对本发明做进一步说明。

附图说明

附图1为本发明实施例中包络分析方法的流程图；

附图2为本发明实施例中采用低通和高通滤波器对信号进行初步分解的示意图；

附图3为本发明实施例中采用树状滤波器结构快速计算谱峭度的示意图；

附图4为本发明实施例中具有内圈故障的滚动轴承振动信号；

附图5为本发明实施例中采用传统包络分析方法对内圈故障滚动轴承振动信号的分析结果；

附图6为本发明实施例中采用本发明包络分析方法对内圈故障滚动轴承振动信号的分析结果；

附图7为本发明实施例中具有外圈故障的滚动轴承振动信号；

附图8为本发明实施例中采用传统包络分析方法对外圈故障滚动轴承振动信号的分析结果；

附图9为本发明实施例中采用本发明包络分析方法对外圈故障滚动轴承振动信号的分析结果。

具体实施方式

实施例，如图1、图2、图3所示，一种基于经验模式分解滤波的包络分析方法，包括以下步骤：

步骤1：利用加速度传感器以采样频率fs测取旋转机械的振动信号x（k）, （k=1,2, …,N），N为采样信号的长度；

步骤2：采用经验模式分解（Empirical Mode Decomposition, EMD）算法将信号x（k）分解成n个分量和一个趋势项之和，即 ，其中，ci（k）代表由EMD算法得到的第i个分量，rn（k）代表由EMD算法得到的趋势项；

步骤3：对ci（k）执行重排操作和替代操作，经重排操作得到的数据用cishuffle（k）表示，替代操作后得到数据用ciFTran（k）表示；

步骤4：对ci（k）、cishuffle（k）和ciFTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，ci（k）的广义Hurst指数曲线用Hi（q）表示；

步骤5：如果Hi（q） 与Hishuffle（q）或Hi（q） 与HiFTran（q）之间的相对误差小于5%，或者Hi（q） 、Hishuffle（q） 和HiFTran（q）三者都不随q而变化，则抛弃对应的ci（k）分量；

步骤6：对剩余的ci（k）分量求和，将该和记为信号经重排和替代滤波后的结果xf1（k）；

步骤7：对xf1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f0和带宽B；

步骤8： 根据中心频率f0和带宽B对xf1（k）进行带通滤波，得到xf2（k）；

步骤9：对信号xf2（k）执行有理样条迭代平滑包络分析，得到信号包络eov（k）；

步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

步骤2中经验模式分解（EMD）算法包括以下步骤：

1） 第一个筛选过程：分别找出数据x（k）的上、下局部极值点，采用三次样条曲线分别拟合上、下局部极值点，得到信号x（k）的局部最大值包络线和局部最小值包络线，再将此两条包络线的相应各点的值取平均，得到一条平均曲线m1；

再求原始信号x（k）与此平均曲线m1的差，即h10=x(k)-m1，至此第一个筛选过程结束；

2）第二个筛选过程：h10重新被当作原始信号，重复上述步骤1），可得到h11= h10-m11，这里参数m11代表h10的均值曲线，重复这一过程j次，直到0.2<SD<0.3时筛选过程停止，这里 ， 此时，h1j= h1(j-1)-m1j，这时可以认为h1j是一个内秉模态函数（Intrinsic Mode Function, IMF），定义第1个IMF为c1=h1j；

3） 从原信号中减去c1，可得r1=x(k)-c1，再将r1当作新的数据，并重复上述两步操作，这样可以得到第2个IMF；

4）重复步骤3）操作可得到一系列IMF，如果rn已经变成一条单调曲线，则筛选过程停止，最终将原始信号分解为如下形式： 。

步骤3中数据重排操作包括以下步骤：

随机打乱分量ci（k）的排列顺序。

步骤3中数据替代操作包括以下步骤：

1） 对分量ci（k）执行离散傅里叶变换，获得分量ci（k）的相位；

2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量ci（k）的原始相位；

3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据ciIFFT（k），求取数据ciIFFT（k）的实部。

步骤4中MFDFA方法包括以下步骤：

1）构造x(k)(k=1,2,…,N)的轮廓Y(i)：

；

x(k)代表权利要求1所述步骤4中的ci（k）或cishuffle（k）或ciFTran（k）；

2）将信号轮廓Y(i)分成不重叠的Ns段长度为s的数据，从数据的反方向以相同的长度分段，得到2Ns段数据；

3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：

；

yv(i)为拟合的第v段数据的趋势，若拟合的多项式趋势为m阶，则记该去趋势过程为（MF-）DFAm；

4） 计算第q阶波动函数的平均值：

；

5）如果x(k)存在自相似特征，则第q阶波动函数的平均值Fq(s)和时间尺度s之间存在幂律关系：

Fq(s)~sH(q)；

当q=0时，步骤4）中的公式发散，这时H(0)通过下式所定义的对数平均过程来确定：

；

6）对步骤5）中的公式两边取对数可得ln[Fq(s)]=H(q)ln(s)+c (c为常数)，由此可以获得直线的斜率H(q)。

步骤7中的谱峭度方法包括以下步骤：

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。