基于自适应Morlet小波参数字典设计的微弱故障检测方法研究

2021-04-28邓飞跃强亚文郝如江马怀祥

振动与冲击 2021年8期

邓飞跃，强亚文，郝如江，马怀祥，高飞

(1.石家庄铁道大学省部共建交通工程结构力学行为与系统安全国家重点实验室，石家庄 050043；2. 石家庄铁道大学机械工程学院，石家庄 050043)

齿轮箱被广泛应用于各种旋转机械设备中，是机械系统中传递动力的重要组成部件。但恶劣的工况往往造成齿轮产生多种故障，甚至引起重大的人员与经济损失。及时、准确地检测出齿轮故障，在工程实际中具有非常积极的意义。

基于信号处理的齿轮故障诊断方法，近年来不断涌现并得到了得到广泛应用。古莹奎等[1]提出一种结合完整集成经验模态分解，信号质量指数与奇异值分解的齿轮故障最优特征提取方法；Zhang等[2]应用一种新的自适应局部振荡特征分解方法，可有效识别齿轮瞬时故障特征频率；佘博等[3]将一种新的增量式监督局部切空间排列算法与支持向量机相结合，用于齿轮箱多故障状态识别。其中，稀疏表示法作为一种新的工具也得到了快速发展。它基于过完备字典，可以将故障信号表示为若干个原子的线性组合，在寻找最优稀疏系数的基础上，可准确提取信号固有的特征结构[4]。但是，如何设计一个有效的稀疏表示字典，一直是一个广泛研究的热点问题。

目前，学习字典与参数字典是原子字典设计的两种主要方法[5]。学习字典通过学习目标信号的特征信息能够自适应地更新字典，主要方法有：K奇异值分解法(K-sigular value decomposition, K-SVD)[6]、最优方向法(method of optimal directions, MOD)[7]、位移不变稀疏编码算法(shift invariant sparse coding algorithm, SISCA)[8]。但学习字典在应用中主要存在两个不足：①只有学习大量样本信号，学习字典才能有效揭示信号中内在固有特征；②由于信号中故障特征微弱，学习得到的字典原子中更多的是强背景噪声干扰成分。参数字典基于固定小波基函数，可充分利用小波波形与冲击波形相似的特点，有效提取信号中微弱的故障冲击特征[9]。Fan等[10]基于Morlet小波，采用相关滤波算法(correlation filtering algorithm，CFA)构造了用于齿轮故障诊断参数字典；Sun等[11]结合Laplace小波与CFA，提出了用于轴承故障诊断的参数脉冲小波字典；He等[12]利用CFA设计了一个包含稳态调制与冲击调制成分的参数字典，用于齿轮箱复合故障诊断。虽然小波基函数类型的选取多种多样，但均通过CFA来确定小波参数值。不可忽视的是，传统CFA方法在实际应用存在着抗噪性差、计算量大等缺陷，限制了参数字典在稀疏表示中的应用。

鉴于上述问题，本文提出了一种新的基于Morlet小波自适应参数字典设计方法。该方法基于局部分割和全局分析思想，采用相关系数与峭度指标综合评价小波函数与目标信号的局部匹配度与全局匹配度，利用鲸鱼优化算法(whale optimization algorithm，WOA)自适应确定小波字典参数。最后，结合正交匹配追踪(orthogonal matching pursuit， OMP)算法对振动信号进行稀疏分解，通过包络谱分析提取齿轮微弱的故障特征。

1 自适应Morlet小波参数字典设计

1.1 稀疏模型及正交匹配追踪

传感器拾取的故障振动信号y表示为

y=h+e

(1)

式中:h为故障特征信号;e为背景噪声。假定原子字典为D，根据信号稀疏模型，y进一步表示为

y=Dx+e

(2)

式中，x={x(1),x(2),…,x(m)}T为稀疏系数集合。信号最稀疏表示即x中大部分系数为零或者接近零，即求解‖x‖0最小l0范数

(3)

针对这个NP-hard问题，我们选择匹配追踪(matching pursuit, MP)等贪婪算法，将求解最小l0范数转化为求解最小l1范数问题：

(4)

OMP是MP的改进算法，通过迭代计算残余信号与字典D原子间的内积，不断对信号进行稀疏分解。假定此时D={dγ}γ∈Γ(Γ={γi,i=1,2,…})是Hilbert空间中一个过完备字典，则故障特征信号h可以分解为

h=〈R0h,dγ0〉dγ0+R1h

(5)

式中:R0h为初始故障特征信号h;R1h为首次迭代计算后的残余信号;dγ0为首次迭代中选择的字典原子。〈R0h,dγ0〉为R0h与dγ0的内积，dγ0和R1h彼此正交。为使R1h尽可能小，需选用与信号内积最大的原子作为最佳匹配原子。M次迭代后，h分解为

(6)

OMP算法中，当前原子需要与之前获取的原子进行施密特正交化处理，正交化过程为

(7)

式中，u0=dγ0。最终信号OMP分解为所选字典原子与相应系数的乘积

(8)

1.2 相关滤波法

原子小波波形与信号中故障冲击成分波形越相似，字典原子越能匹配信号中的冲击故障特征。Morlet小波是一个双边对称小波，其波形与齿轮故障产生的瞬态冲击成分较为相似，因此本文选用Morlet小波构造原子参数字典。Morlet小波表达式为

(9)

式中:γ={f,ζ,τ}为Morlet小波参数集合;f∈R+为频率;ζ∈[0,1)⊂R+为黏滞阻尼比;τ为时移参数。小波参数直接决定了小波的波形特征，是决定原子小波与信号故障冲击特征匹配的重要因素。假设离散参数f,ζ,τ分别属于集合F，Z，T

(10)

遍历所有的参数集合即可得到全部的Morlet小波集

Ψ={ψγ(t):γ∈F×Z×T}={ψ(f∈F,ζ∈Z,τ∈T)}

(11)

相关系数cγ能够定量表征Morlet小波Ψγ(t)与所分析信号y(t)之间的相似程度[13]

(12)

通过上述分析可知，CFA方法需要在三个参数集合范围内进行遍历迭代运算，计算量大且非常耗时。而且当信号中包含大量背景噪声时，相关函数的计算会受到噪声成分的影响，往往不能准确表征小波与信号故障冲击特征的匹配程度，因此CFA的抗噪性较差。

1.3 自适应优化的参数字典设计

在CFA中，所选小波支撑区间往往与所分析信号长度相等[14]。然而，Morlet小波中主要冲击成分只占全部波长的一小部分。因此，本文根据小波能量正态分布3δ准则重新确定Morlet小波支撑区间为

(13)

式中:ψfiζi(t)为小波参数fi,ζi构成的Morlet小波;L为分析信号长度;W为重新确定的小波支撑区间。根据新的小波支撑区间，对信号进行分割，构造空间矩阵为

(14)

式中,Y1,Y2,…,YL-W+1为各个子信号。计算ψfiζi(t)与每个子信号的相关函数得到Ci=[ci,1,ci,2,…,ci,L-W+1]，计算该序列的峭度为

(15)

式中：E(·)为数学期望；μ和σ分别为Ci的均值和标准差。遍历参数集合f⊂F,ζ⊂Z，得到F×Z个峭度，最大峭度值为

Kbest=max(K1,K2,…,KF×Z)

(16)

选取Kbest对应的小波参数f′,ζ′为最优小波基参数。在此基础上，取小波支撑区间等于所分析信号长度构建Morlet小波，通过逐点时移构成最优原子参数字典。cγ描述了单个小波与信号的局部匹配，而由cγ构成的序列Ci=[ci,1,ci,2,…,ci,L-W+1]描述了单个小波在逐点时移下与信号的整体匹配。该序列故障冲击特征越明显，所选小波与所分析信号的局部匹配与整体匹配程度越好。

WOA是一种新的自然启发式优化算法，它通过模拟鲸鱼“螺旋气泡网”狩猎行为实现参数过程自适应优化[15]，主要包括：寻找猎物、围捕猎物、泡网捕食等阶段，具体实现过程请参见文献[16]。本文所提的基于WOA的自适应参数字典设计方法具体步骤如下:

步骤1设置Morlet小波寻优参数范围，黏滞阻尼比范围为ζ∈[0,1)，频率范围为f∈(0,fs/2]，fs为所分析信号的采样频率；

步骤2参考文献[17]，初始化WOA中各项参数，其中鲸群数量A=50，最大迭代次数M=20，优化参数个数为2；

步骤3随机选取小波参数构建Morlet小波，根据本文所提方法，得到对应小波参数下的峭度值，并以此作为WOA的适应度函数；

步骤4执行WOA算法，不断迭代更新鲸鱼狩猎位置，进行bubble-net捕食操作；

步骤5通过多次迭代，判断满足收敛准则后，输出最优适应度函数对应的小波参数；

步骤6基于最优小波参数构建Morlet小波，逐点时移构成原子参数字典；

步骤7使用OMP方法对所分析信号进行稀疏分解，对分解后信号进行包络解调分析，提取故障特征。

流程图如图1所示。

2 仿真分析

本节通过仿真齿轮微弱故障信号，验证所提方法的有效性。传统故障信号模拟中，往往添加白噪声作为背景噪声，然而实际工程中白噪声并不存在，因此本文在仿真信号加入白噪声和有色噪声来更接近实际情况。仿真信号表达式为

(17)

式中:h(t)为齿轮故障信号;n(t)为添加信噪比为-5 dB的高斯白噪声;r(t)为有色噪声，均值为0标准差为1的高斯白噪声通过滤波后加入。仿真信号中频率f=20 Hz，阻尼比ζ=0.008，齿轮故障冲击周期T=0.25 s，冲击起始时刻τ0=2 s，信号采样频率fs=600 Hz，采样点数为8 192。

原始故障信号及添加背景噪声后的信号,如图2所示。添加背景噪声后，信号中故障冲击成分被淹没，无法准确识别。采用本文所提方法，基于Morlet小波自适应识别小波参数，WOA输出的最大峭度值为6.165 2，对应结果见表1。构建小波参数字典，通过OMP信号进行稀疏分解，结果如图3(a)所示。信号波形中清晰的显示出周期性故障冲击成分，与原始仿真信号一致。然后，采用传统CFA进行分析对比，鉴于CFA方法运算非常耗时，小波参数搜寻范围及步长并没有设置很细致，分别为：f∈[0.1∶10∶fs/2]Hz，ζ∈[0.001∶0.000 5∶0.01]，t∈[0∶1/fs∶end]。结果如图4所示。最大相关系数为0.110 9，对应小波参数位置用虚线标出，具体结果见表1。从表1可知，本文所提方法得到的Morlet小波参数与初始设置非常接近，而CFA得到的参数结果与设置值相差较大。同时，由于所提方法参数遍历运算更为简洁，运算时间更少。CFA识别小波参数后，构建原子参数字典，OMP处理结果如图3(b)所示。图中并未提取出准确的周期性故障冲成分。相比传统CFA方法，所提方法识别原子小波参数更为准确，抗噪性更好，运算效率更高。

图2 齿轮故障仿真信号

图3 仿真信号分析结果

图4 仿真信号CFA识别小波参数结果

表1 仿真信号分析结果

3 试验分析

3.1 试验介绍

齿轮故障是在QPZZ-II旋转机械故障实验台上完成的。齿轮箱从动大齿轮出现断齿故障，主动轮与从动轮齿数分别为55和75，齿轮故障实验台结构及故障齿轮,如图5、图6所示。试验中，大齿轮与小齿轮的齿数比为75∶55，输入轴转速为870 r/min，采样频率为fs=5 120 Hz，采样点数为8 192点。通过计算可知，齿轮啮合频率为fm=797.5 Hz，故障齿轮的故障特征频率为fg=10.63 Hz。

图5 齿轮故障实验台

图6 齿轮断齿

3.2 实例介绍

齿轮故障信号如图7所示。时域波形中背景噪声较多，故障冲击成分并不规律，频谱中也只能提取出并不明显的一阶啮合频率。首先，采用常用的Hilbert包络解调对原信号进行分析，具体实现过程参考文献[18]，得到原信号Hilbert幅值为

(18)

图7 齿轮故障信号

图8 本文所提方法分析结果

表2 齿轮故障信号信号分析结果

在CFA方法中，参考Zi等的研究，为保证在低阻尼比时有较高分辨率，小波参数划分不均匀。小波参数搜寻范围及步长分别为：f∈[0.1∶20∶fs/2]，ζ∈{[10-6∶10-6∶0.001]∪[0.001∶0.001∶0.02]}，t∈[0∶1/fs∶end]。小波参数识别结果如图9所示。表2中是具体数值。通过对比两种方法结果，可以发现CFA与本文所提方法参数识别结果相差较大，运算时间也大幅增加。构建字典，经过OMP处理后的结果如图10所示。包络谱中虽然出现了齿轮故障频率，但在噪声频率干扰下非常微弱，无法准确识别。通过上述分析可知，小波参数选取对参数字典的设计至关重要，所提方法可以有效弥补传统CFA方法的不足。

图9 齿轮故障信号CFA识别小波参数结果

图10 基于CFA方法的分析结果

3.3 对比分析

首先采用文献[19]中小波降噪方法分析齿轮故障信号进行对比。小波函数选取“db4”小波，采用固定史坦无偏估计(rigrsure)准则确定阈值，小波硬阈值和软阈值降噪后的故障信号波形和频谱,如图11和图12所示。我们采用Hilbert包络解调方法对两结果分别进行分析，包络谱见图11(c)和图12(c)。相比齿轮故障原始故障信号波形，小波降噪后时域信号中噪声成分有所减少，频谱中部分频率成分也受到抑制，但包络谱中依然无法提取出与故障特征频率相关的频率成分。面对信号中强背景噪声的干扰，小波阈值降噪方法效果有限，无法检测出齿轮微弱的故障特征信息。

图11 齿轮故障信号小波硬阈值降噪结果

图12 齿轮故障信号小波软软值降噪结果

采用文献[20]所提K-SVD学习字典，与本文所提方法进行对比。K-SVD方法中相关参数设置为：字典原子长度和数量为n=100，N=4 065，线性组合原子数为L=10，迭代次数I=10。学习字典中前10个典型原子波形，如图13所示。从图13可知，各个原子波形中并未发现明显的瞬时故障冲击成分。同样，基于OMP方法对信号进行稀疏分解，得到的分析结果如图14所示。从图14可知，时域波形中故障冲击成分微弱且杂乱无章，包络谱中也没有提取出明显的齿轮故障特征频率及倍频成分。面对强背景噪声干扰下齿轮微弱故障特征提取，所提方法要优于K-SVD学习字典方法。