CN110084208A - 一种自适应降噪并避免阶次混叠的计算阶次跟踪方法 - Google Patents

Abstract

本发明公开了一种自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，根据信号转速信息及预计最高分析阶次定义裕量频率；对信号进行VMD预分解，将中心频率低于裕量频率的模态保留，将高于裕量频率的抛弃，以滤除信号中的高频噪声以及非分析阶次带宽内的高阶成分；计算重构信号的排列熵PE：PE用来表示时间序列的随机程度，PE的值越大，时间序列也就越随机；采用差分进化算法对VMD参数进行优化，得到参数后自适应生成重构信号；计算重采样阶次，对获得的重构信号进行计算阶次跟踪，得到重采样信号后进行FFT以得到信号的阶次谱。本发明用于处理原始振动信号，以自适应降低采集的振动信号中的噪声干扰，凸显故障信息。

Description

一种自适应降噪并避免阶次混叠的计算阶次跟踪方法

技术领域

本发明涉及一种自适应降噪并避免阶次混叠的计算阶次跟踪方法，具体涉及一种基于经差分进化优化的变分模态分解和计算阶次跟踪方法，属于非平稳信号的处理技术领域。

背景技术

旋转机械变工况运行时的振动信号不仅包含机械设备零部件自身的振动信息，也包含转速与负载信息，这些信息融合在一起，使得振动信号呈现十分复杂的非平稳特征，使得信号产生严重的幅值和频率调制。如果对这些非平稳振动信号直接采用传统的频谱分析技术如傅里叶分析将产生严重的频率模糊现象，导致对相应问题的误判或漏判。

阶次跟踪是解决这一问题的一种最常用的方法，目前已经发展出了多种形式。计算阶次跟踪(Computed Order Tracking，COT)是阶次跟踪技术的一种，其比传统阶次跟踪更为简化，而且精度比较高，其通过等角度采样将时域非平稳信号转化为角域平稳信号，能将因速度变化产生的影响从振动信号中去除，克服传统频谱分析方法难以准确提取变工况条件下特征信息的问题，此外，其对信号中与转速无关的分量也有一定的抑制作用。如果原始信号中具有较大的噪声，容易对重要的特征信息造成干扰，影响实际判断。而且其在角域重采样时需要指定重采样阶次，如果设置的重采样阶次比较小将会导致模态混叠，增加重采样阶次可以避免这一现象，但是会影响算法效率，更重要的是，实际中往往不知道要分析的最高阶次，即便是增加了采样阶次，模态混叠也不能完全避免。

变分模态分解(Variational Mode Decomposition，VMD)可以通过迭代搜寻变分模型的最优解确定每个模态的中心频率和带宽，自适应地实现信号各模态的有效分离。而且其非递归、变分分解模式，能有效克服模态混叠和端点效应等缺点。该算法具有良好的噪声鲁棒性，其在降低信号噪声、凸显故障特征信息等方面均具有良好的作用，很适合用于处理非平稳信号。

发明内容

本发明所要解决的技术问题是：变工况下运用阶次跟踪技术时存在的阶次混叠问题。

为了解决上述技术问题，本发明通过以下技术方案实现：

一种自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，包括以下步骤：

步骤1)：根据信号转速信息及预计最高分析阶次定义裕量频率：

假设某段时间内参考轴的最高转速为n_max，最低转速为n_min，振动信号的采样频率为f_s，所要分析的信号的最高阶次为O_max，那么临界频率f_c的计算公式为：

那么，在阶次带宽O_max内，信号频率f的变化如式2所示：

f(O≤Omax)≤f_c 式2；

考虑到裕量m，定义裕量频率f_m为：

f_m＝f_c+m 式3；

步骤2)：对信号进行VMD预分解，将中心频率低于裕量频率的模态保留，将高于裕量频率的抛弃，以滤除信号中的高频噪声以及非分析阶次带宽内的高阶成分：

假设振动信号x(t)被分解成K个模态，K个模态的中心频率ω_k的集合表示为：

{ω_k}＝{ω₁,ω₂,…,ω_k} 式4；

其中，k＝1,2,…,K；

中心频率ω_k满足式5的模态将会被保留下来：

ω_k≤f_m 式5；

其中，k＝1,2,…,i；

重构信号y(t)的表达式如下：

其中，u_j为保留模态；

步骤3)：计算重构信号的排列熵PE：PE用来表示时间序列的随机程度，PE的值越大，时间序列也就越随机；实际的振动信号往往呈现出一定的规律性，而噪声是杂乱无章的，所以PE可以用于描述信号中的噪声含量，噪声越严重，PE的值就越大。

对于重构信号的时间序列{Y(i),i＝1,2,…,N}，对其进行相空间重构，可以得到矩阵：

其中，j＝1，2，…，G。G表示重构分量的个数；d表示嵌入维数；τ表示时延。

矩阵X中的每一行可看作一个重构分量，共有G＝N-(d-1)·τ个重构分量。将重构矩阵X中的每一行{y(j)y(j+τ)…y(j+(d-1)·τ)}按照升序重新排列，j₁,j₂,…,j_d表示重构分量中各个元素所在列的索引，即

y(i+(j₁-1)τ)≤y(i+(j₂-1)τ)≤…≤y(i+(j_d-1)τ) 式8；

如果重构分量中存在相等的值，即有

y(i+(j_p-1)τ)＝y(i+(j_q-1)τ) 式9；

此时就按照j_p和j_q值的大小来排序，也就是当j_p≤j_q时，有

y(i+(j_p-1)τ)≤y(i+(j_q-1)τ) 式10；

因此，对于任意一个时间序列重构所得的矩阵X中每一行都可以得到一组符号序列

S(l)＝{j₁,j₂,…,j_d} 式11；

式中，l＝1,2,…,G且G≤d！，d维相空间映射不同的符号序列{j₁,j₂,…,j_d}总共有d！种，符号序列S(l)是其中的一种排列。计算每一种符号序列出现的概率P₁,P₂,…,P_G，此时，时间序列Y(i)的G种不同符号序列的排列熵可以按照Shannon熵的形式定义为

通过ln(m！)对式13进行标准化，可得：

H_P(m)＝PE(m)/ln(m！) 式13；

步骤4)：采用差分进化算法对VMD参数进行优化，得到参数后按照步骤1)和步骤2)自适应生成重构信号：

采用差分进化算法对VMD进行改进，其中，差分进化算法的适应度函数选取为重构信号的排列熵PE，为了避免DE在自适应寻优的过程中，盲目追求PE的最小值而漏掉信号中的重要模态，裕量频率f_m与在组成重构信号的所有模态中最大的中心频率的差值不得超过偏差e，如果该差值超过偏差e，该重构信号所计算的PE的值在原基础上加上常数C＝1，如果该差值小于偏差e，那么PE的值仍保留为原信号的值；

为保证算法在初期种群具有良好的多样性并在后期能较快地收敛到最优解，同时具有较好的局部搜索能力，差分进化算法的缩放因子F和交叉因子CR的计算公式为：

由于分解个数与惩罚因子均为整数，且惩罚因子步长往往为数十倍才有明显变化，因此本发明在传统的DE/rand/1的差分进化策略的基础上进行了改进，以便算法运行时能够以合适的部长搜寻个体数目，提高算法的运行效率。改进后的VMD分解参数模态数目K以及惩罚因子α的差分策略如下：

v(1,m)＝round(x(1,r₁)+F×(x(1,r₁)-x(1,r₁))) 式16；

v(2,m)＝round(x(1,r₁)/10+F×(x(1,r₁)-x(1,r₁))/10)×10 式17；

其中，round表示返回一个数值，该数值是按照指定的小数位数进行四舍五入运算的结果；

步骤5)：计算重采样阶次，对获得的重构信号进行计算阶次跟踪，得到重采样信号后进行FFT以得到信号的阶次谱：

自适应得到的重构信号中，满足式13的频率成分都已经被滤除掉了：

f≥f_m 式18；

因此，重构信号中包含的最高阶次为：

在角域中存在着与Shannon采样定理相似的角域采样定理：角域重采样阶次应至少为信号中最高分析阶次的2倍，即：

O_s≥2O_max 式20；

联立式19和式20，角域重采样阶次满足式21：

最后确定好角域重采样阶次，对重构信号进行等角度重采样，然后对得到的角域重采样信号进行FFT即可得到重构信号的阶次谱。

优选地，所述步骤1)中的裕量m根据信号的实际情况以及经验选取，其值一般不大于100Hz。

优选地，所述步骤3)中嵌入维数d的取值范围为5-8，优选为6；时延τ设置为1。

优选地，所述步骤4)中根据经验可选取偏差e为裕量m的2-3倍。

本发明将变分模态分解、差分进化算法以及排列熵进行融合，提出了DE-VMD算法，用于处理原始振动信号，以自适应降低采集的振动信号中的噪声干扰，凸显故障信息，同时克服变工况下运用阶次跟踪技术时存在的阶次混叠问题。

附图说明

图1为DE-VMD的算法流程图；

图2为信号x₁的时域图；

图3为转速曲线；

图4为噪声信号的时域图；

图5为信号x的时域图；

图6为重构信号的阶次谱；

图7为原始信号的阶次谱。

具体实施方式

为使本发明更明显易懂，兹以优选实施例，并配合附图作详细说明如下。

本发明提供的一种自适应降噪并避免阶次混叠的计算阶次跟踪方法，包括以下步骤：

步骤1)：根据信号转速信息及预计最高分析阶次定义裕量频率：

假设某段时间内参考轴的最高转速为n_max，最低转速为n_min，振动信号的采样频率为f_s，所要分析的信号的最高阶次为O_max，那么临界频率f_c的计算公式为：

那么，在阶次带宽O_max内，信号频率f的变化如式2所示：

f(O≤Omax)≤f_c 式2；

考虑到裕量m，定义裕量频率f_m为：

f_m＝f_c+m 式3；

步骤2)：对信号进行VMD预分解，将中心频率低于裕量频率的模态保留，将高于裕量频率的抛弃，以滤除信号中的高频噪声以及非分析阶次带宽内的高阶成分：

假设振动信号x(t)被分解成K个模态，K个模态的中心频率ω_k的集合表示为：

{ω_k}＝{ω₁,ω₂,…,ω_k} 式4；

其中，k＝1,2,…,K；

中心频率ω_k满足式5的模态将会被保留下来：

ω_k≤f_m 式5；

其中，k＝1,2,…,i；

重构信号y(t)的表达式如下：

其中，u_j为保留模态；

步骤3)：计算重构信号的排列熵PE：PE用来表示时间序列的随机程度，PE的值越大，时间序列也就越随机；实际的振动信号往往呈现出一定的规律性，而噪声是杂乱无章的，所以PE可以用于描述信号中的噪声含量，噪声越严重，PE的值就越大。

对于重构信号的时间序列{Y(i),i＝1,2,…,N}，对其进行相空间重构，可以得到矩阵：

其中，j＝1，2，…，G。G表示重构分量的个数；d表示嵌入维数；τ表示时延。

矩阵X中的每一行可看作一个重构分量，共有G＝N-(d-1)·τ个重构分量。将重构矩阵X中的每一行{y(j)y(j+τ)…y(j+(d-1)·τ)}按照升序重新排列，j₁,j₂,…,j_d表示重构分量中各个元素所在列的索引，即

y(i+(j₁-1)τ)≤y(i+(j₂-1)τ)≤…≤y(i+(j_d-1)τ) 式8；

如果重构分量中存在相等的值，即有

y(i+(j_p-1)τ)＝y(i+(j_q-1)τ) 式9；

此时就按照j_p和j_q值的大小来排序，也就是当j_p≤j_q时，有

y(i+(j_p-1)τ)≤y(i+(j_q-1)τ) 式10；

因此，对于任意一个时间序列重构所得的矩阵X中每一行都可以得到一组符号序列

S(l)＝{j₁,j₂,…,j_d} 式11；

式中，l＝1,2,…,G且G≤d！，d维相空间映射不同的符号序列{j₁,j₂,…,j_d}总共有d！种，符号序列S(l)是其中的一种排列。计算每一种符号序列出现的概率P₁,P₂,…,P_G，此时，时间序列Y(i)的G种不同符号序列的排列熵可以按照Shannon熵的形式定义为

通过ln(m！)对式13进行标准化，可得：

H_P(m)＝PE(m)/ln(m！) 式13；

步骤4)：采用差分进化算法对VMD参数进行优化，得到参数后按照步骤1)和步骤2)自适应生成重构信号：

采用差分进化算法对VMD进行改进，其中，差分进化算法的适应度函数选取为重构信号的排列熵PE，为了避免DE在自适应寻优的过程中，盲目追求PE的最小值而漏掉信号中的重要模态，裕量频率f_m与在组成重构信号的所有模态中最大的中心频率的差值不得超过偏差e，如果该差值超过偏差e，该重构信号所计算的PE的值在原基础上加上常数C＝1，如果该差值小于偏差e，那么PE的值仍保留为原信号的值；

为保证算法在初期种群具有良好的多样性并在后期能较快地收敛到最优解，同时具有较好的局部搜索能力，差分进化算法的缩放因子F和交叉因子CR的计算公式为：

由于分解个数与惩罚因子均为整数，且惩罚因子步长往往为数十倍才有明显变化，因此本发明在传统的DE/rand/1的差分进化策略的基础上进行了改进，以便算法运行时能够以合适的部长搜寻个体数目，提高算法的运行效率。改进后的VMD分解参数模态数目K以及惩罚因子α的差分策略如下：

v(1,m)＝round(x(1,r₁)+F×(x(1,r₁)-x(1,r₁))) 式16；

v(2,m)＝round(x(1,r₁)/10+F×(x(1,r₁)-x(1,r₁))/10)×10 式17；

其中，round表示返回一个数值，该数值是按照指定的小数位数进行四舍五入运算的结果；

步骤5)：计算重采样阶次，对获得的重构信号进行计算阶次跟踪，得到重采样信号后进行FFT以得到信号的阶次谱：

自适应得到的重构信号中，满足式13的频率成分都已经被滤除掉了：

f≥f_m 式18；

因此，重构信号中包含的最高阶次为：

在角域中存在着与Shannon采样定理相似的角域采样定理：角域重采样阶次应至少为信号中最高分析阶次的2倍，即：

O_s≥2O_max 式20；

联立式19和式20，角域重采样阶次满足式21：

最后确定好角域重采样阶次，对重构信号进行等角度重采样，然后对得到的角域重采样信号进行FFT即可得到重构信号的阶次谱。

实施例

本实施例采用了如式22所示的调幅调频信号x₁(t)模拟齿轮局部故障时的理想振动信号，如图2所示；采用了式25所示的一倍调制频率f_r(t)模拟速度的变化情况如图3所示；为模拟噪声污染，加入噪声分量η(t)如图4所示，同时为模拟未知高阶成分干扰，分别加入高阶分量x₂(t)和x₃(t)，从而得到要分析的信号x(t)如图4所示。式22中信号的载波频率被f_r(t)所调制，载波频率为调制频率f_r(t)的6倍。设置采样频率为2048，分析的时间间隔为[0,1]。则：

x₁＝{1+cos[2π×(4t³+t²+10t+12)]}·cos[2π×6(4t³+t²+10t+12)] 式22；

x₂(t)＝cos[2π×14(4t³+t²+10t+12)] 式23；

x₃(t)＝cos[2π×20(4t³+t²+10t+12)] 式24；

f_r(t)＝12t²+2t+10 式25；

x(t)＝x₁(t)+x₂(t)+η(t) 式26；

对式22所示的调幅调频信号进行积化和差分析，得到的表达式如式27所示，可以看出式22中信号可由3个信号分量组成，其中第6阶为主要阶次，5阶、7阶分量为边频带成分，且幅值均为6阶分量的1/2。

本实施例的具体实施步骤如下所述：

步骤1：根据转速曲线可知，信号的最高转速为n_max＝1440，最低转速为n_min＝600。本例所要分析的最高阶次为7。那么可计算其临界频率f_c＝7×1440/60＝168Hz。

设置裕量m＝12，那么可计算的裕量频率f_m＝180Hz。

步骤2：设置种群数目NP＝24，最大迭代次数G_m＝30，根据相关经验，DE参数分别设置为：F_max＝0.9，F_min＝0.4，CR_max＝0.6，CR_min＝0.1，K_h＝15，K_l＝2，α_h＝6000，α_l＝100。根据设定参数，执行DE-VMD算法，最终得到的重构信号如图5所示。

步骤3：信号经过重构以后，包含的最高阶次为O_r＝180×60/600＝18。由采样定理取重采样阶次O_s＝36，对获得的重构信号进行角域重采样，得到重采样信号后进行FFT以得到信号的阶次谱，如图6所示。

本发明同时给出了未采用本发明方法，直接应用采样定理确定重采样阶次的计算阶次跟踪方法，该方法将采样阶次取最大分析阶次的两倍以上，即O_s≥14，取O_s＝18，计算的阶次谱如图7所示。从图7可以发现，由于强噪声干扰，第5阶成分已经难以准确判断，而且在第4阶和第2阶出现了明显的分量，说明发生了阶次混叠，这是因为采样阶次为18时，原信号中的第14阶和第20阶通过第9阶混叠到了第4阶和第2阶。而在图6中，第20阶高阶成分已经被滤除，没有出现在阶次谱中，除了主要阶次及其边频带，未知的14阶成分也被提取出来了，没有发生阶次混叠现象，并且信号中噪声干扰被滤除掉，各阶次特征十分明显。因此本发明提出的方法能有滤除信号中的噪声，凸显特征信息，而且该算法能够提取出实际信号中隐含的阶次成分，避免阶次混叠现象。

Claims (4)

1.一种自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，包括以下步骤：

步骤1)：根据信号转速信息及预计最高分析阶次定义裕量频率：

假设某段时间内参考轴的最高转速为n_max，最低转速为n_min，振动信号的采样频率为f_s，所要分析的信号的最高阶次为O_max，那么临界频率f_c的计算公式为：

那么，在阶次带宽O_max内，信号频率f的变化如式2所示：

f(O≤Omax)≤f_c 式2；

考虑到裕量m，定义裕量频率f_m为：

f_m＝f_c+m 式3；

步骤2)：对信号进行VMD预分解，将中心频率低于裕量频率的模态保留，将高于裕量频率的抛弃，以滤除信号中的高频噪声以及非分析阶次带宽内的高阶成分：

假设振动信号x(t)被分解成K个模态，K个模态的中心频率ω_k的集合表示为：

{ω_k}＝{ω₁,ω₂,…,ω_k} 式4；

其中，k＝1,2,…,K；

中心频率ω_k满足式5的模态将会被保留下来：

ω_k≤f_m 式5；

其中，k＝1,2,…,i；

重构信号y(t)的表达式如下：

其中，u_j为保留模态；

步骤3)：计算重构信号的排列熵PE：对于重构信号的时间序列{Y(i),i＝1,2,…,N}，对其进行相空间重构，可以得到矩阵：

其中，j＝1，2，…，G。G表示重构分量的个数；d表示嵌入维数；τ表示时延。

矩阵X中的每一行可看作一个重构分量，共有G＝N-(d-1)·τ个重构分量。将重构矩阵X中的每一行{y(j) y(j+τ) … y(j+(d-1)·τ)}按照升序重新排列，j₁,j₂,…,j_d表示重构分量中各个元素所在列的索引，即

y(i+(j₁-1)τ)≤y(i+(j₂-1)τ)≤…≤y(i+(j_d-1)τ) 式8；

如果重构分量中存在相等的值，即有

y(i+(j_p-1)τ)＝y(i+(j_q-1)τ) 式9；

此时就按照j_p和j_q值的大小来排序，也就是当j_p≤j_q时，有

y(i+(j_p-1)τ)≤y(i+(j_q-1)τ) 式10；

因此，对于任意一个时间序列重构所得的矩阵X中每一行都可以得到一组符号序列

S(l)＝{j₁,j₂,…,j_d} 式11；

式中，l＝1,2,…,G且G≤d！，d维相空间映射不同的符号序列{j₁,j₂,…,j_d}总共有d！种，符号序列S(l)是其中的一种排列。计算每一种符号序列出现的概率P₁,P₂,…,P_G，此时，时间序列Y(i)的G种不同符号序列的排列熵可以按照Shannon熵的形式定义为

通过ln(m！)对式13进行标准化，可得：

H_P(m)＝PE(m)/ln(m！) 式13；

步骤4)：采用差分进化算法对VMD参数进行优化，得到参数后按照步骤1和步骤2自适应生成重构信号：

采用差分进化算法对VMD进行改进，其中，差分进化算法的适应度函数选取为重构信号的排列熵PE，为了避免DE在自适应寻优的过程中，盲目追求PE的最小值而漏掉信号中的重要模态，裕量频率f_m与在组成重构信号的所有模态中最大的中心频率的差值不得超过偏差e，如果该差值超过偏差e，该重构信号所计算的PE的值在原基础上加上常数C＝1，如果该差值小于偏差e，那么PE的值仍保留为原信号的值；

差分进化算法的缩放因子F和交叉因子CR的计算公式为：

VMD分解参数模态数目K以及惩罚因子α的差分策略如下：

v(1,m)＝round(x(1,r₁)+F×(x(1,r₁)-x(1,r₁))) 式16；

v(2,m)＝round(x(1,r₁)/10+F×(x(1,r₁)-x(1,r₁))/10)×10 式17；

其中，round表示返回一个数值，该数值是按照指定的小数位数进行四舍五入运算的结果；

步骤5)：计算重采样阶次，对获得的重构信号进行计算阶次跟踪，得到重采样信号后进行FFT以得到信号的阶次谱：

自适应得到的重构信号中，满足式13的频率成分都已经被滤除掉了：

f≥f_m 式18；

因此，重构信号中包含的最高阶次为：

在角域中存在着与Shannon采样定理相似的角域采样定理：角域重采样阶次应至少为信号中最高分析阶次的2倍，即：

O_s≥2O_max 式20；

联立式19和式20，角域重采样阶次满足式21：

最后确定好角域重采样阶次，对重构信号进行等角度重采样，然后对得到的角域重采样信号进行FFT即可得到重构信号的阶次谱。

2.如权利要求1所述的自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，所述步骤1)中的裕量m不大于100Hz。

3.如权利要求1所述的自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，所述步骤3)中嵌入维数d的取值范围为5-8；时延τ设置为1。

4.如权利要求1所述的自适应降噪并避免阶次混叠的计算阶次跟踪方法，其特征在于，所述步骤4)中偏差e为裕量m的2-3倍。