基于多指标优化TQWT和TEO的轴承声发射故障诊断

2023-02-28陈俊潼周凤星严保康汪峰

机床与液压 2023年3期

陈俊潼，周凤星，严保康，汪峰

(武汉科技大学信息科学与工程学院，湖北武汉 430081)

0 前言

滚动轴承是被广泛应用在旋转机械中的支承元件，其性能和工作状态对机械设备的正常运转至关重要[1]。因此针对滚动轴承故障诊断，尤其是早期故障诊断的研究具有重大现实意义。

声发射(Acoustic Emission，AE)信号检测技术是一种新型无损检测技术，具有响应频带宽、灵敏度高的特点，可以有效提取早期微弱故障信号。但是实际采集的轴承声发射信号带有大量噪声和调制成分，有必要对其进行处理，提高信噪比并解调，以便提取故障特征信息。

对轴承故障声发射信号这种非平稳信号，小波变换降噪是一种有效的方法。孙守保等[2]通过AR模型和最小熵解卷积(Minimum Entropy Deconvolution,MED)滤波结合复Morlet小波变换成功提取了声发射信号的双冲击特征。张瑞等人[3]使用经验小波变换(Empirical Wavelet Transform，EWT)构建小波滤波器组，实现对声发射信号的消噪和特征提取。徐嗣嘉等[4]将小波包能量谱与共振解调法结合，通过分析子频带能量提取了声发射信号特征。以上研究证明了小波变换在轴承故障声发射信号处理的有效性。但是经典小波变换对信号的分解效果依赖于小波基的选取，如果小波基与故障特征不匹配会导致故障特征成分的损失，并在子频带中引入无关噪声。可调Q因子小波变换(TunableQfactor Wave-let Transform,TQWT)是在小波变换基础上改良的时频分析方法，相比于经典小波变换，TQWT可对品质因子Q、冗余度r、分解层数J三个参数进行调节，实现小波基与待提取故障特征的最佳匹配，从成分复杂的信号中提取出微弱的故障特征分量。目前TQWT方法主要应用在振动信号上，关于声发射信号TQWT研究还很少。孔运等人[5]设计自适应TQWT滤波器算法，结合冲击特征指标实现了信号降噪和特征提取。余发军和周凤星[6]依据谱峭度指标用TQWT对信号进行分解，有效提取了早期故障冲击成分。张龙等人[7]利用最小熵解卷积突出信号冲击特征，再用TQWT对其进行分解重构，得到轴承故障特征信号。

Teager能量算子能够跟踪信号的瞬时能量，并对信号解调，与Hilbert解调法相比，能量算子解调法能抑制端点效应，近年来越来越多的研究者将其运用在轴承故障诊断中。杨斌等人[8]将自相关变换与能量算子解调法结合，提取出故障特征。郭璐等人[9]通过快速傅里叶变换得到Teager能量谱图，输入卷积神经网络，得到故障诊断模型。

基于以上分析，本文作者提出一种多指标优化TQWT与Teager能量算子相结合的故障诊断方法，通过改进的TQWT方法对原始声发射信号进行自适应分解和频带选取，对选出的子频带降噪后重构信号,再使用Teager能量算子对重构信号解调，对解调后的信号进行频谱分析得到故障特征频率，实现故障诊断。将该方法分别应用于仿真故障信号与实测信号，验证了该方法的有效性。

1 理论介绍

1.1 TQWT可调Q因子小波变换基本理论

TQWT是在离散二进小波变换基础上改进的可以灵活调节Q因子的离散小波变换。与经典小波变换不同，TQWT直接在频域上利用小波基的频域响应对被处理信号进行滤波。在频域上设计双通道过采样滤波器组，将信号分解为高频部分和低频部分，经过尺度变换后，高频部分作为本层输出，低频部分作为下层的输入继续分解。这样通过高、低通滤波器组的层层分解将信号分解为不同子频带,TQWT分解的信号具有可重构性。具体分解过程如图1所示，重构即分解的逆向过程。

图1 TQWT分解过程Fig.1 TQWT decomposition process

图1中每一层的输出W为TQWT分解的小波系数序列，可单支重构为子频带时间序列。定义第j层子频带的中心频率为fc,带宽为Bw[10]，计算公式分别为

(1)

式中:fs为采样频率;α、β分别为滤波器的高通变换尺度和低通变换尺度，定义为

(2)

式中:Q为小波的品质因子，Q越大，小波时域振荡次数越多；r是冗余系数，表现为分解所得子频带的重合度；J是最大分解层数，影响低频段分解性能，分解完成后会获得J+1个子频带。三参数(Q,r,J)共同定义了TQWT的性质和特性。最大分解层数J可以由Q、r确定，表达式为

(3)

式中：⎣·」为向正无穷取整；N为原信号长度。

1.2 Teager能量算子

20世纪80年代，TEAGER等提出一种非线性算子TEO(Teager Energy Operator)，又被称为Teager能量算子，最初被应用在声音信号增强建模上。连续信号x(t)的能量算子ψc定义如下:

(4)

ψ[x(n)]=x2(n-1)-x(n)x(n-2)

(5)

对于某个具有时变幅值和时变相位的声发射传感器响应信号h(t)[11]：

h(t)=exp(-2πξfnt)sin(2πfnt)

(6)

式中：fn为传感器固有频率；ξ是阻尼系数。使用Teager能量算子对其处理得：

ψc[h(t)]=[exp(-2πξfnt)]2·(2πfn)2

(7)

由公式(7)可知:Teager能量算子能很好地对原始信号进行解调，相比信号能量定义为幅值的平方，Teager能量算子增加了瞬时频率平方的乘积[12]。对于轴承故障声发射信号这种高频冲击信号，Teager能量算子可以有效增强故障特性。

另一方面，Teager能量算子对噪声十分敏感，因此常将其与其他方法结合以减弱噪声的影响[13]。

2 多指标优化的TQWT方法

为了从轴承早期故障声发射信号提取出微弱故障特征，提出一种多指标优化的TQWT方法。该方法可以对原始信号进行自适应分解，并筛选出包含故障信息的子频带进行降噪和重构，有利于提高后续能量算子解调效果。

2.1 参数范围确定

根据公式(3)可知，最大分解层数J可由Q、r确定；根据TQWT基本性质可知，Q的取值不小于1，同时为了防止Q过大导致小波共振属性过高产生奇异信号，将Q的取值范围设为[1,20]；综合考虑计算量和实际需要，将r的取值设为3[14-15]。

2.2 参数自适应选择

TQWT参数的选取可以视为设计中心频率和带宽不同的滤波器组，通过这些滤波器得到一系列子频带。研究表明声发射信号的波形熵可以作为故障特征指标，子频带波形信息熵越小，说明该频带中含有越多故障成分[16]。峭度是常用的描述信号冲击成分指标，峭度越大信号中的冲击性成分越多。对于长度为N的序列{x(n)},峭度的公式为

(8)

式中：σ为序列标准差；u为平均值。结合峭度指标与波形熵指标，使子频带峭度最大化的同时使波形熵尽可能小。对于选定的参数(Q,r)，分解得到J+1个子频带，其中第j层子频带对应改进峭度波形熵比Cj定义为

(9)

式中：xij为第j层子频带的第i个元素；N为子频带序列长度；Kj表示第j层子频带峭度。取所有子频带峭度波形熵比中的最大值C为参数(Q,r)定义下的故障特征指标。在给定范围内遍历参数的取值，C取到最大值时，对应参数为最佳分解参数。

2.3 子频带选择和降噪

分解完成后得到J+1个子频带，由于声发射信号频带宽，具有多频率解调特性，为了获得相对完整的故障信息并去除噪声，将3种故障特征融合为一个新的故障特征指标对子频带进行筛选。为了使前后优化方向一致，选用峭度作为故障指标之一，峭度虽然对冲击敏感但是容易受到局部极值影响，峰度系数可以反映信号整体波形的极端程度，不易受局部极点影响。峰度系数M表示为

(10)

稀疏值S可以表示信号的故障稀疏程度，为

(11)

分别求每个子频带信号的峭度K、峰度M、稀疏度S，得到3个特征指标序列{K(i)}、{M(i)}、{S(i)},使用这3个序列构成多指标的状态矩阵：

X=(xij)3×(J+1)=

(12)

式中：xij表示矩阵中第i行第j列元素。分别对矩阵中同一类指标求信息熵：

(13)

式中：Ei存在一个最大值ln(J+1)。信息熵小的指标会给出更有价值的故障信息,因此在特征融合过程中给信息熵小的指标更大权重[17-18]。为了消除不同指标间的差异，将各子带的3个指标采用min-max标准化后根据信息熵的大小对其赋予不同的权重，计算得到每个子频带的融合特征参数。

(14)

(15)

(16)

式中：Wi表示各特征指标的权重；Fj为各子频带的融合特征指标。求出Fj的平均值，将Fj值大于平均值2倍的子频带挑选出来，对这些分量采用无偏似然估计软阈值去噪，之后使用去噪后的分量重构信号。设序列为x(t)，序列长度为n，阈值σ的具体计算方式为

(17)

(18)

3 故障诊断总流程

提出的多指标优化TQWT与Teager能量算子结合的故障诊断方法具体步骤如下：

(1)使用峭度波形熵比优化的自适应TQWT算法对故障信号进行分解,得到多个子频带。

(2)根据融合特征指标选出包含故障特征成分的子频带，对其降噪后重构信号。

(3)使用Teager能量算子对重构信号解调，得到解调谱，从解调谱中得到故障特征频率，实现特征提取和故障诊断。

4 仿真信号检验

轴承外圈故障声发射信号可描述为近似周期性双冲击响应，公式如下：

(19)

式中：x(t)表示滚动轴承故障产生的声发射信号；x1(t)表示第一个冲击响应；x2(t)表示第二个冲击响应；ε(t)表示噪声；h(t)由公式(6)确定，式中ξ取0.05，fn取3 kHz；采样速率设定为20 kHz，采样时间1 s，并对其加入高斯白噪声使得信噪比为-13 dB。原始仿真信号时域波形如图2(a)所示，加入噪声仿真信号如图2(b)所示。对加入噪声的仿真信号进行频谱分析，低频段如图2(c)所示，可以看出信号已经完全被噪声淹没，从频谱图上也无法得到明显故障信息。

图2 仿真信号时频图

使用文中的方法对信号进行处理，根据前文确定的品质因子Q变化为[1,20], 冗余系数r取3,搜索步长为1，在范围内搜索使得C最大的Q值。C随Q变化情况如图3所示,在Q=8时C取到最大值1.116 2，分解得到83个子频带。根据融合特征指标下对子频带进行筛选，图4列出了所有子频带对应的融合特征值。

图3 C随Q变化情况(仿真数据)Fig.3 Case that C changes with Q(simulation signal)

图4 子频带筛选标准(仿真数据)Fig.4 Sub-bands screening criteria(simulation signal)

融合特征值的平均值为0.003 1，选出特征值大于平均值2倍的子频带，j=15、16的子频带满足条件。对选出的子频带降噪后重构信号如图5所示，使用Teager能量算子对重构后的信号进行解调得到能量算子序列如图6所示，对能量算子序列做频谱分析，得到图7所示低频段的频谱图，可以清楚看到故障特征频率及其谐波，与设计的外圈故障频率10 Hz符合，实现特征提取和故障诊断。

图5 去噪后的重构信号Fig.5 Reconstructed signal after denoising

图6 重构信号Teager能量算子序列Fig.6 Reconstructed signal Teager energy operator sequence

图7 低频段Teager能量算子频谱Fig.7 Teager energy operator spectrum in low-frequency stage

5 轴承故障诊断实验

5.1 实验说明

为了进一步证明方法的有效性，设计实验进行验证。实验轴承为N205深沟球轴承，轴承节径为38.5 mm，滚动体直径6.5 mm，滚动体个数为13，接触角为0°，在滚动体上设置了宽度约0.4 mm、深度约0.5 mm的裂纹缺陷。实验设备为QPZZ-II故障诊断模拟实验平台，试验台如图8所示。

图8 故障诊断模拟实验平台Fig.8 Fault diagnosis simulated experimental platform

实验时将声发射传感器置于轴承上方，以1 MHz的采样频率对声发射故障信号采样，每次采样131 072个点，电机转速为600 r/min，对应转频10 Hz。

5.2 滚动体故障信号分析

经过计算得滚动体故障频率fOR=57.5 Hz，为了模拟工业环境中的强噪声环境，不对原信号做硬件滤波处理，并利用软件向原信号加入5 dB白噪声，所得信号时域波形如图9(a)所示。可以看出信号存在冲击干扰而且故障成分被噪声淹没。对其进行频谱分析，低频段频谱如图9(b)所示，无法通过该频谱分辨故障特征。

图9 强噪声背景下滚动体故障信号时频图

采用文中的方法对其进行处理。首先使用自适应TQWT对信号进行分解，C随品质因子Q变化情况如图10所示。

图10 C随Q变化情况(实验数据)Fig.10 Case that C changes with Q(test signal)

最优分解参数Q=3时，Cmax=4.123 3，计算出J=49。分解得到50个子频带，子频带对应的融合特征参数如图11所示。

图11 子频带筛选标准(实验数据)Fig.11 Sub-bands screening criteria(test signal)

经计算融合特征参数平均值为0.024 9，大于两倍平均值的子频带有j=25、26、27、28、29、30、32。对选出的子频带进行无偏似然软阈值去噪后重构信号，重构信号如图12所示。

图12 重构信号Fig.12 Reconstructed signal

使用Teager能量算子对重构后的信号进行处理得到能量算子序列，并进行频谱分析，结果如图13所示。

图13 Teager能量算子序列(a)与频谱(b)Fig.13 Teager energy operator sequence(a) and spectrum(b)

对Teager能量算子频谱进行分析，发现第一个峰值对应的频率为10 Hz，与电机转频符合，之后的峰值出现在61 Hz及其倍频附近。考虑到电机运行过程中的转速波动与滚动体打滑等因素影响，提取出的峰值频率与计算所得滚动体fOR=57.5 Hz是符合的，成功提取出滚动体故障特征并实现故障诊断。

5.3 与现有方法对比分析

现有的改进TQWT方法往往只使用了一种特征指标对TQWT进行优化，如文献[5]使用了能量熵作为优化指标，文献[6-7]使用了峭度作为优化指标。使用单一指标往往存在局限性而且在工程中鲁棒性不强。而且现有的方法都是对振动信号进行处理，只选取出所有子频带中唯一的最优子频带，由于声发射信号具有多频带解调特性，仅选取单一频带会导致大量故障信息被遗漏。

以上文处理的滚动体故障信号为例，仅使用峭度指标对TQWT参数进行优化并选取最优子频带，所得最优分解参数Q=5，选出的解调频带仅有j=27子频带，重构所得信号如图14所示。可以明显看出图14中的故障冲击比图12中的微弱。对其进行能量算子解调后进行频谱分析，如图15所示。该频谱中的峰值与图13(b)相比不够明显，特征频率的倍频处峰值也较少。