机械早期及复合故障的集成噪声重构经验模式分解方法

摘要

本发明公开了一种机械早期及复合故障的集成噪声重构经验模式分解方法，它首先对信号中的噪声分量进行估算和重采样，并重构若干输入信号。利用噪声极值点的变化来改善信号不连续问题，为避免模式混淆的信号分解提供前提。然后，对这些重构输入信号进行经验模式分解(EMD)得到相应的基本模式分量(IMF)。最后，利用随机序列统计均值为零的原理，将所有对应的IMF做总体平均运算，消除信号中噪声分量的影响，从而达到信号降噪与特征增强的目的。本发明非常适合于信噪比相对较低的机械早期及复合故障信号处理及特征提取，识别精度高，运算实时性好，简单易行，便于工程实践中使用。

说明书

技术领域

本发明涉及大型机械设备（航空发动机、汽轮发电机组、连铸连轧机组、风电装备等）故障诊断技术，具体涉及一种机械设备的早期故障和复合故障的特征提取与识别方法。 

背景技术

机械设备一旦出现事故，将带来巨大的经济损失和人员伤亡。在工程实践中，我们期望对故障的发生与发展能够做到防微杜渐，不期望亡羊补牢般地处理事故。因此，研究开发有效的早期故障和复合故障诊断技术具有重要的工程应用价值。早期故障包含两方面含义，其一是指处于早期阶段的微弱故障或潜在故障，具有症状不明显、特征信息微弱、信噪比低等特点；其二是从物理意义上讲，某一故障是另一故障的早期阶段，并随着时间推移进一步诱发复合故障。复合故障是指关键设备同时发生两个或多个互相关联、交叉影响的故障，其动态信号一般表现为多个故障特征的相互耦合与调制。早期故障和复合故障的特征提取、故障定位与定量比较困难，已成为国内外本领域研究的难点和热点。 

近年来兴起的第二代小波变换、多小波变换、双树复小波变换、谱峭度等方法在早期故障特征提取上展现出广阔的应用前景。但是这些方法中如何构造和选择与故障特征波形相匹配、且具有优良性质的合适基函数成为故障诊断的关键所在。同时，这些方法在进行抑制噪声的同时，可能会“过扼杀”信号中的重要微弱故障特征信息，造成故障特征信息的丢失。 

经验模式分解(Empirical Mode Decomposition，EMD)方法可根据信号本身具有的时间尺度特征自适应地将其分解为一系列本征模式分量(intrinsic mode function,IMF)的线性和，已在设备早期故障诊断领域得到应用。但传统EMD方法对信号进行分解时，时常会出现由于某一或几个时间尺度的IMF不连续，产生模式混淆的现象，从而影响EMD方法进一步应用。为克服模式混淆问题，集成经验模式分解（Ensemble Empirical Mode Decomposition,EEMD）方法应运而生。该方法从噪声利用的角度，对信号添加白噪声来改善传统EMD方法中的模式混淆现象。然而EEMD应用效果依赖于添加的白噪声大小和总体平 均次数，目前尚缺乏可靠、通用的噪声协助准则。同时，对于信噪比相对较低的机械故障信号，EEMD中的平均过程仅能抵消所添加的白噪声，不能消除信号中的原有噪声从而提取信号中早期及复合故障特征。 

发明内容

本发明的目的在于提供一种应用于机械早期及复合故障特征提取的集成噪声重构经验模式分解(ensemble noise-reconstructed EMD)的方法，该方法不仅可继承EMD方法在自适应分解非平稳信号上的优势，同时它利用原信号中固有噪声分量可实现信号降噪并改善EMD中的模式混淆问题。 

为达到以上目的，本发明是采取如下技术方案予以实现的： 

一种机械早期及复合故障的集成噪声重构经验模式分解方法，其特征在于，包括下述步骤： 

1）对所采集到的原始信号x(t)中的噪声分量进行估算： 

a.首先对x(t)进行EMD分解，得到各IMF{c_k(t),k＝1,···,n}和残差分量r_n(t)，其中，k代表IMF个数； 

b.计算各个{c_k(t)}的能量{E_k,k＝1,···n}，设定第一个IMF分量c₁(t)为纯白噪声，即c₁(t)中的纯白噪声能量

c.根据纯白噪声下IMF能量公式，分别计算各个{c_k(t)}得到各IMF在置信区间95%的纯白噪声能量和置信区间99%下的纯白噪声能量$>\{\hat{E}{n}_{k,99\%},k=\mathrm{2,3}...n\};$>

d.将各IMF的能量{E_k}与相应的纯白噪声能量对比：若 $>\hat{E}{n}_{k,95\%}\le E_k\le \hat{E}{n}_{k,99\%},$>或$>|\hat{E}{n}_{k,95\%}-E_k|\ge \alpha ,$>或$>|\hat{E}{n}_{k,99\%}-E_k|\ge \alpha ,$>其中α∈[1,2]，为一个给定容差值，则将c_k(t)作为待处理的含噪IMF；否则，c_k(t)作为估算的故障特征分量； 

e.将以上待处理的含噪IMF组成{c_l(t),l＝1,···}； 

f.对待处理的含噪IMF{c_l(t),l＝1,···}进行如下反硬阈值处理得到提纯噪声分量

其中，为c_l(t)阈值处理后的第i个样本，且阈值门限其中，N为c_l(t)的长度，C为常数且C∈[0.3,0.5]； 

g.合并所有提纯噪声分量得到信号x(t)中估算噪声分量

2）通过标准均匀分布下的随机重排采样点方式对估算噪声分量进行重采样得到然后重构相应的输入信号

3）对进行EMD分解得到IMF分量{c_j,k(t),k＝1,···,n}及残差分量{r_j,n(t)}； 

4）重复步骤2）和步骤3）直至满足以下运算停止准则： 

a.假设进行了m次EMD运算，则得到平均运算后的IMF分量和残差分量如下： 

$>{\bar{c}}_k\left(t\right)=\frac{\underset{j=1}{\overset{m}{\Sigma }}{c}_{j,k}\left(t\right)}{m},{\bar{r}}_n\left(t\right)=\frac{1}{m}\underset{j=1}{\overset{m}{\Sigma }}{r}_{j,k}\left(t\right)$>

并重构得到降噪信号$>\bar{s}\left(t\right)=\underset{k=1}{\overset{n}{\Sigma }}{\bar{c}}_k\left(t\right)+{\bar{r}}_n\left(t\right);$>

b.计算得到平均运算后的噪声分量及其能量

c.根据IMF噪声估算方法获得估算噪声分量的能量

d.计算相对噪声误差比

若得到的相对噪声误差比例er小于给定的误差允许值ε，则停止运算，否则重复运算直至满足该误差，其中ε∈[0.01,0.1]； 

5）当满足停止准则时，对步骤3）中计算得到的所有IMF分量进行平均运算处理，得到最终平均IMF分量和平均残差分量： 

$>{\tilde{c}}_k\left(t\right)=\frac{\underset{j=1}{\overset{r}{\Sigma }}{c}_{j,k}\left(t\right)}{r},{\tilde{r}}_n=\frac{\underset{j=1}{\overset{r}{\Sigma }}{r}_{j,n}\left(t\right)}{r}$>

从和中提取早期及复合故障特征。 

上述方法中，所述对信号x(t)进行EMD分解包含以下步骤： 

1）把x(t)作为待处理信号，确定其所有局部极值点，然后将所有极大值点和所有极小值点分别用三次样条曲线连接起来，得到x(t)的上、下包络线，使x(t)的所有数据点都处于这两条包络线之间； 

2）从待处理信号x(t)中减去其上、下包络线均值m(t)，得到h₁(t)＝x(t)-m(t)； 

3）检测h₁(t)是否满足基本模式分量的两个条件，如果不满足，则把h₁(t)作为待处理信号，重复上述操作，直至h₁(t)是一个基本模式分量c₁(t)＝h₁(t)； 

4）分解出第一个基本模式分量c₁(t)之后，从x(t)中减去c₁(t)，得到剩余值序列r₁(t)＝x(t)-c₁(t)； 

5）把r₁(t)作为新的“原始”信号重复上述操作，依次可得第二、第三直至第n个基本模式分量，记为c₁(t),c₂(t),···,c_n(t)，这个处理过程在满足预先设定的停止准则后即可停止，最后剩下原始信号的余项r_n(t)。 

本发明继承了EMD方法在自适应分解非平稳信号上的优势，并具有以下优点： 

1）与传统基于噪声抑制的降噪方法不同，本发明利用了原信号中固有噪声分量的相互抵消作用进行降噪，可避免“过扼杀”重要微弱故障特征信息的不足，从而有效识别故障信号中的早期及复合故障特征。 

2）与传统EEMD方法不同，本发明不必再额外添加白噪声序列，仅利用原有噪声的重采样便可有效改善和消除模式混淆现象，避免了人为添加白噪声的干扰。 

3）本发明识别精度高，运算实时性好，简单易行，非常适合于信噪比相对较低的机械早期及复合故障信号处理及特征提取。 

附图说明

以下结合附图及具体实施方式对本发明作进一步的详细描述。 

图1为实施本发明方法的仿真信号图。 

图2为本发明方法的流程框图。 

图3为图2信号中各IMF的{E_k},和

图4为图2信号中待处理的含噪IMF的反硬阈值处理。其中，（a）图为各IMF的反硬阈值处理（横线表示阈值门限）；（b）图为c₅(t)反硬阈值处理的放大图（阈值门限内的浅色信号代表反硬阈值处理后信号）。 

图5为图2信号中估计噪声和实际噪声的对比图。 

图6为图2信号中估计噪声和实际噪声的概率密度分布。 

图7为图2信号中采用本发明的分析结果。其中，（a）图为采用本发明分析得到的IMF分量；（b）图为IMF分量c₂的放大图；（c）图为IMF分量c₄的放大图。 

图8为图2信号中采用EMD的分析结果。其中，（a）图为采用EMD分析得到的IMF分量；（b）图为IMF分量c₂的放大图；（c）图为IMF分量c₄的放大图。 

图9为图2信号中采用EEMD的分析结果。其中，（a）为采用EEMD分 析得到的IMF分量；（b）图为IMF分量c₂的放大图；（c）图为IMF分量c₄的放大图。 

图10为实例1齿轮早期故障中振动信号[（a）图]及其频谱[（b）图]。 

图11为采用本发明对图10故障的分析结果。其中，（a）图为采用本发明分析得到的IMF分量；（b）图为IMF分量c₂的放大图；（c）图为IMF分量c₄的放大图。 

图12为实例2齿轮复合故障中振动信号[（a）图]及其频谱[（b）图]。 

图13为采用本发明对图12故障的分析结果。齿轮复合故障中采用本发明的分析结果。其中，（a）图为采用本发明分析得到的IMF分量；（b）图为IMF分量c₂的放大图；（c）图为IMF分量c₄的放大图。 

具体实施方式

本发明具体方法包含以下内容： 

对所采集到的原始信号x(t)中的噪声分量进行估算： 

首先对x(t)进行EMD分解，得到各IMF{c_k(t),k＝1,···,n}和残差分量r_n(t)，其中，k代表IMF个数；计算各个{c_k(t)}的能量{E_k,k＝1,···n}，设定第一个IMF分量c₁(t)为纯白噪声，即c₁(t)中的纯白噪声能量

根据纯白噪声下IMF能量公式，分别计算各个{c_k(t)}得到各IMF在置信区间95%的纯白噪声能量和置信区间99%下的纯白噪声能量 将各IMF的能量{E_k}与相应的纯白噪声能量对比：若$>\hat{E}{n}_{k,95\%}\le E_k\le \hat{E}{n}_{k,99\%},$>或$>|\hat{E}{n}_{k,95\%}-E_k|\ge \alpha ,$>或$>|\hat{E}{n}_{k,99\%}-E_k|\ge \alpha ,$>其中α∈[1,2]，为一个给定容差值，则将c_k(t)作为待处理的含噪IMF；否则，c_k(t)作为估算的故障特征分量； 

将以上待处理的含噪IMF组成{c_l(t),l＝1,···}； 

对待处理的含噪IMF{c_l(t),l＝1,···}进行如下反硬阈值处理得到提纯噪声分量

其中，为c_l(t)阈值处理后的第i个样本，且阈值门限其中，N为c_l(t)的长度，C为常数且C∈[0.3,0.5]； 

合并所有提纯噪声分量得到信号x(t)中估算噪声分量

通过标准均匀分布下的随机重排采样点方式对估算噪声分量进行重采样得到然后重构相应的输入信号

对进行EMD分解得到IMF分量{c_j,k(t),k＝1,···,n}及残差分量{r_j,n(t)}； 

重复噪声分量估算和EMD分解，直至满足以下运算停止准则： 

假设进行了m次EMD运算，则得到平均运算后的IMF分量和残差分量如下： 

$>{\bar{c}}_k\left(t\right)=\frac{\underset{j=1}{\overset{m}{\Sigma }}{c}_{j,k}\left(t\right)}{m},{\bar{r}}_n\left(t\right)=\frac{1}{m}\underset{j=1}{\overset{m}{\Sigma }}{r}_{j,k}\left(t\right)$>

并重构得到降噪信号$>\bar{s}\left(t\right)=\underset{k=1}{\overset{n}{\Sigma }}{\bar{c}}_k\left(t\right)+{\bar{r}}_n\left(t\right);$>

计算得到平均运算后的噪声分量及其能量根据IMF噪声估算方法获得估算噪声分量的能量

计算相对噪声误差比

若得到的相对噪声误差比例er小于给定的误差允许值ε，则停止运算，否则重复运算直至满足该误差，其中ε∈[0.01,0.1]； 

当满足停止准则时，对多次EMD分解得到的所有相应的IMF分量进行平均运算处理，得到最终平均IMF分量和平均残差分量： 

$>{\tilde{c}}_k\left(t\right)=\frac{\underset{j=1}{\overset{r}{\Sigma }}{c}_{j,k}\left(t\right)}{r},{\tilde{r}}_n=\frac{\underset{j=1}{\overset{r}{\Sigma }}{r}_{j,n}\left(t\right)}{r}$>

从和中提取早期及复合故障特征。 

对信号x(t)进行EMD分解包含以下步骤： 

把x(t)作为待处理信号，确定其所有局部极值点，然后将所有极大值点和所有极小值点分别用三次样条曲线连接起来，得到x(t)的上、下包络线，使x(t)的所有数据点都处于这两条包络线之间； 

从待处理信号x(t)中减去其上、下包络线均值m(t)，得到h₁(t)＝x(t)-m(t)； 

检测h₁(t)是否满足基本模式分量的两个条件，如果不满足，则把h₁(t)作为待处理信号，重复上述操作，直至h₁(t)是一个基本模式分量c₁(t)＝h₁(t)； 

分解出第一个基本模式分量c₁(t)之后，从x(t)中减去c₁(t)，得到剩余值序列r₁(t)＝x(t)-c₁(t)； 

把r₁(t)作为新的“原始”信号重复上述操作，依次可得第二、第三直至第n个 基本模式分量，记为c₁(t),c₂(t),···,c_n(t)，这个处理过程在满足预先设定的停止准则后即可停止，最后剩下原始信号的余项r_n(t)。 

为了验证本发明所述方法的正确性，给出两个具体实例进一步说明。 

实例1 

构造包含正弦分量（模拟旋转机械特征信号）、冲击分量（模拟局部故障特征信号）和噪声分量的混合信号x(t)为 

x(t)＝s₁(t)+s₂(t)+λn(t) 

s₁(t)＝sin(160πt),$>s_2\left(t\right)=\underset{k=1}{\overset{8}{\Sigma }}f(t-\mathrm{kT}),$>f(t)＝e^-200｜t｜sin(1200πt) 

其中,T＝0.125s,n(t)为均值为0方差为1的高斯白噪声，且λ＝0.35。对该仿真信号采用采样频率f_s＝5120Hz以及点数n＝5120进行仿真，构造的仿真信号如图1所示。 

对该仿真信号采用本发明进行特征提取，按照图2所示的流程进行分析。首先对仿真信号进行EMD分解，计算得到各IMF分量的{E_k,k＝1,···n}， $>\{\hat{E}{n}_{k,95\%},k=\mathrm{2,3}...n\}$>和$>\{\hat{E}{n}_{k,99\%},k=\mathrm{2,3}...n\},$>如图3所示。根据$>\hat{E}{n}_{k,95\%}\le E_k\le \hat{E}{n}_{k,99\%},$>或者或者（其中，参数α＝2），选择{c_l(t),l＝1,2，3,5,6,8,9}组成待处理的含噪IMF。对这些IMF进行反硬阈值处理（其中，参数C=0.5），如图4所示。图中，横线代表各IMF的阈值门限{T_l}，深色曲线代表原信号{c_l(t)}，阈值门限内的浅色曲线代表反硬阈值处理后的信号 将以上阈值处理后的信号合并得到信号估算噪声分量如图5所示。图中，实线表示仿真信号中所估计的噪声分量虚线表示仿真实际噪声分量λn(t)，可以看到，估计噪声分量和实际噪声分量的幅值大小基本一致。同时，图6为仿真信号中估计噪声分量（实线）和仿真实际噪声分量λn(t)（虚线）的概率密度分布，可以看到，两者的概论密度分布也非常相似，也证明了本发明中的噪声估计技术的有效性。对估算的噪声分量进行随机重排的重采样得到噪声分量然后重构相应的含噪信号并进行EMD分解。计算噪声重构及总体平均后的噪声分量和相对噪声误差比er。重复以上步骤直至er小于给定的误差允许值ε（参数ε=0.1）。最后平均所有IMF分量，输出的分析结果如图7所示。从图7中可以看出，IMF分量c₂和c₄分别对应于仿真分量s₂(t)和s₁(t)，且c₂中的周期性冲击序列非常明显。采用峭度指标ks来衡量本发明对于冲击序列的提取效果，并且采用s₁(t)和c₄的相 关系数e来评价正弦信号的完整性，也即评估改善模式混淆现象的效果。本发明中，c₂的ks=17.9129，且c₄的e=0.9699。 

对该仿真信号采用EMD和EEMD方法进行分析，分别如图8和图9所示。从图8中可以看到，周期性冲击序列在c₂中不易识别，且在c₄和c₅中发生了模式混淆现象（如图中双箭头所示）。从图9中可以看到，虽然EEMD可以改善模式混淆现象，但是c₂中提取的周期性冲击特征远不如本发明中的分析结果。同时，在EMD中峭度指标ks=5.1656且相关系数e=0.9203，而在EEMD中ks=5.9622且e=0.9600，这些指标都不如本发明的指标值。因此，本实施案例中冲击序列分析说明本发明可以很好地提高信号的信噪比，提取与识别微弱特征信息，并且本案例中的正弦分量分析说明本方法可以明显改善EMD方法中的模式混淆现象。 

实施案例2： 

在齿轮试验台上采用齿轮试验来验证本发明对早期故障和复合故障特征提取能力。其中，大齿轮齿数223，小齿轮齿数48，其中小齿轮有一个轮齿被切除用于计数。采样频率设定为25kHz，大齿轮转速为145r/min，小齿轮转速为659r/min。 

图10为齿轮早期裂纹故障的振动信号及其频谱。按照图2所示的流程，采用本发明对齿轮早期故障信号进行分析，输出的分析结果如图11所示。从图11中可以看出，在c₂和c₄中分别出现明显的周期性冲击特征。c₄中周期性冲击是由小齿轮中缺少的一个轮齿所引起。而c₂中周期性冲击间隔为0.4144s，接近于大齿轮的旋转频率，说明该大齿轮存在局部故障，如齿根裂纹或疲劳剥落，与齿轮实际故障状态相符。 

实施案例3： 

图12为上述齿轮试验中齿轮存在2个断齿时的振动信号及其频谱。同样按照图2所示的流程，采用本发明对齿轮复合故障信号进行分析，输出的分析结果如图13所示。其中，c₄中周期性冲击特征依旧对应于小齿轮缺少的轮齿所产生的振动特征。而c₂呈现为若干包含双周期性冲击特征信号，与大齿轮存在两个断齿故障相符。同时，两组周期性冲击的间隔为0.0483s，说明两个断齿在圆周上相隔26个轮齿，与实际相符。