朱丹宸 ，张永祥，赵 磊，朱群伟

(海军工程大学 动力工程学院，武汉 430033)

滚动轴承在旋转机械中得到广泛使用，其工作状态直接影响着整个机器设备的运行稳定性，因此，对滚动轴承进行状态监测和故障诊断就具有重要意义。就目前而言，振动信号分析是轴承故障诊断的常用手段。滚动轴承故障一般包括内圈故障、外圈故障以及滚动体故障，这些故障产生的振动信号往往是非平稳、非线性的，同时周期性故障冲击一般都淹没在强背景噪声之中[1]，如何降低原始振动信号中的噪声成分，提取出轴承故障特征就成了主要研究内容。

TVD(Total Variation Denoising)是一种基于最优问题的降噪方法，最早被用于去除图像中的噪声成分[2]，其本质上是一种信号的稀疏表示方法，随后，TVD方法得到了广泛的研究并在一维信号处理中得到了运用。由于TVD是通过将代价函数最小化来获得最终输出，文献[3]提出利用Majorization-Minimization(MM)算法对TVD进行求解从而构成TV-MM算法，然而，该算法中，参数λ的选取对最终滤波结果有着较大影响。Zhang等[4]将变分模态分解以及TV-MM算法进行结合，利用加权峭度求解最优参数λ，并证明该方法优于一些经典的去噪算法。Yi等[5]提出了一种二阶TVD算法，将惩罚函数引入算法之中，通过排列熵和相关系数选取参数λ，提高了降噪效果。

近来，Alwodai等[6]利用MSB(Modulation Signal Bispectrum)分析电机中的电流信号，展现出比传统功率谱更好的频谱特性，实现了电机轴承的故障诊断。文献[7]同样利用电流信号，通过MSB分析，实现了齿轮磨损状况的监测。Rehab等[8]利用MSB从包络信号中提取出故障特征，体现出比传统包络分析更强的抗噪与故障监测能力。文献[9]提出了一种更为稳定的MSB解调方法，实现了滚动轴承的故障特征提取，且分析效果优于传统的Kurtogram算法。

受到强背景噪声的影响和信号自身特点的制约，TVD能够消除部分噪声，但同时也有一些残留噪声影响故障特征提取效果。本文提出将TVD与MSB相结合，利用MSB分析TVD处理后的滤波信号，进一步抑制噪声的干扰，提取出轴承故障特征。首先，二阶TVD被用于处理采样得到的振动信号，利用包络谱相关峭度选取参数λ，获得较优滤波结果；然后利用MSB分析滤波后的信号，通过故障特征频率成分占比p选取最佳的5个载波频率切片进行平均得到复合切片谱，提取出轴承故障特征；最后，通过分析复合切片谱，判断故障类型。仿真和实验分析表明，该方法能够有效抑制随机噪声，提高故障特征提取效果。

1 二阶全变分去噪

1.1 基本原理

全变分去噪可以看成是一个数值优化过程，包含二次数据保真项和凸正则化项，常用的全变分过程的基础是通过一阶或者二阶差分实现对原始信号的稀疏表示。本文选用二阶差分定义信号的全变分。假设一维信号x(n)，(0≤n≤N-1)，定义信号x(n)的二阶全变分为

(1)

式中：D2x为二阶差分，其中D2可以用矩阵形式表示为

(2)

‖·‖p为p阶范数，且‖x‖p可以表示为

(3)

当p=1以及p=2时，信号的一阶与二阶范数可分别表示为

‖x‖1=|x1|+|x2|+…+|xN|

(4)

(5)

假设信号y(n)由有效成分x(n)与白噪声w(n)组成，即

y(n)=x(n)+w(n)

(6)

由此，二阶全变分去噪模型可以通过一个选优过程表示

(7)

式中：F(x)为目标函数；λ>0为正则化参数。

在获得目标函数之后，本文利用MM算法求解式(7)，结果可以表示为

(8)

利用式(8)就可以实现信号的全变分去噪，但通过观察式(7)和式(8)可知，正则化参数λ对降噪过程会产生一定的影响，当λ→0时，降噪后的结果x将无限趋近于原始信号y，当λ→∞时，信号的保真度将会下降。因此，如何选取最为合适的参数λ是接下来将要研究的内容。

1.2 最优参数的选取方法

轴承故障信号在时域范围内包含周期性的冲击特征，所以相关峭度(Correlated Kurtosis,CK)作为一种衡量信号中冲击特性的指标，由于它既保留了峭度的特性，也包含了相关函数的特性，常被用于轴承故障诊断之中，其计算公式为[10]

(9)

式中：xn为信号序列；N为信号的长度；T为感兴趣信号的长度；M为偏移的周期个数。同时，轴承故障冲击特性在频域范围内会表现为故障特征频率及其倍频。由此，本文利用频域相关峭度自适应选取全变分去噪中的参数λ。频域相关峭度可以通过原信号的包络谱定义为

(10)

式中：E(x)n为信号序列xn的包络谱；此时T选取为故障特征频率。当信号的包络谱中包含明显的故障特征频率成分时，频域相关峭度较大，相反，当故障特征频率成分没有在信号的包络谱中得到明显体现时，频域相关峭度较小。

2 调制信号双谱

2.1 基本原理

MSB作为一种考虑边频带的双谱分析方法能对调制信号进行有效分析，该方法可以较好地抑制随机噪声和非周期成分的干扰，清晰反映出信号中的调制成分。MSB在频域范围内，利用原始信号x(t)的离散傅里叶变换形式X(f)可以表示为

BMS(fc,fx)=
E〈X(fc+fx)X(fc-fx)X*(fc)X*(fc)〉

(11)

式中：E〈·〉为数学期望；fc为载波频率；fx为调制频率。对MSB进行归一化可得

(12)

2.2 最优切片谱选取

为了进行有效的故障特征提取，选择合适的分析频段是十分重要的，在利用MSB进行分析时，也即是选取合适的fc切片位置，而一般来说，故障特征较为明显的切片位置有多个，选取多个切片可以综合各切片处的故障特征信息，同时也可以减小存在于单个切片中的干扰成分。由此，本文提出利用前三阶故障特征频率占比p来选取最佳的5个fc切片，并对其进行平均得到复合切片谱，从而提取出轴承故障特征。

设f1，f2，f3分别为轴承故障的特征频率及其二倍频和三倍频，则在fc切片位置，特征频率占比p可以定义为

(13)

(14)

3 故障特征提取流程

基于TVD与MSB的滚动轴承故障特征提取步骤如下，流程如图1所示：

步骤1获取振动信号，在0.1～5内，间隔为0.1选取λ并计算相应的TVD处理结果；

步骤2以最大包络谱相关峭度为准则确定参数λ并计算最优去噪结果；

步骤3对步骤2得到的去噪信号进行MSB分析；

步骤4根据特征频率成分占比p选取MSB中的5个最优载波频率切片；

步骤5对步骤4中得到的5个切片谱进行平均，得到复合切片谱，提取出轴承故障特征，判断故障形式。

4 轴承故障仿真信号分析

为了验证本文前一部分所述故障特征提取方法是正确而又有效的，先利用仿真信号进行研究分析，以滚动轴承内圈故障为例构造仿真信号

(15)

图1 滚动轴承故障特征提取流程图Fig.1 Flow chart of fault feature extraction for rolling element bearings

式中：轴的转频fr为42 Hz；Ai为以1/fr为周期的幅值调制；n(t)为随机白噪声；r为系统的阻尼系数且设定为0.05；两个相邻冲击的间隔T为1/185 s；ti为第i各周期内，由滚动体滑移引起的延迟，ti=0.01T；B(t)为谐波分量；fm=100；φA=π/2；φw=0；A0，B0，C0分别为2，0.5，0.5；系统自然频率fn为3 200 Hz；采样频率fs为32 768 Hz；仿真时间为1 s。

为了分析TVD的去噪效果，将随机白噪声分别设为0，-5 db与-10 db，分别得到仿真信号1、仿真信号2和仿真信号3，其时域波形如图2所示。由于受到噪声的干扰，周期性冲击成分在三组仿真信号中都不能得到体现，且所添加白噪声越强，干扰越明显。

利用TVD处理仿真信号，根据图3(a)、图3(c)和图3(e)给出的三组λ值变化曲线，选取的最优参数λ分别为0.7，1.3和2.1，得到的降噪结果如图3(b)、图3(d)和图3(f)所示。对比图3和图2，虽然经过TVD处理后的时域波形中能够看出噪声的减小，但周期性冲击成分仍然难以得到体现。

图4对比了三组仿真信号在降噪前后的包络谱，分析仿真信号1，降噪前后的包络谱中都能明显识别冲击频率185 Hz以及其倍频370 Hz，555 Hz，轴承故障特征得到了有效提取，且图4(b)中，围绕故障特征频率及其倍频，间隔为转频的调制边频带更为明显，说明利用TVD能够进一步降低噪声提高故障特征提取效果；分析仿真信号2，原始信号的包络谱无法识别故障特征频率，通过TVD处理后，降噪信号的包络谱能够识别故障特征频率及其倍频成分，但白噪声对故障特征提取的干扰是较为明显的；分析仿真信号3，降噪前后信号的包络谱都不能识别故障特征频率及其倍频。综上所述，当信号中的背景噪声较强时，TVD的降噪效果并不理想，需要利用其他分析手段进一步加以处理。

图2 仿真信号的时域波形Fig.2 Time domain waveform of simulation signals

图3 TVD处理后的结果Fig.3 Time domain waveform after TVD processing

图4 原始仿真信号与去噪信号的包络谱Fig.4 Envelope spectrum of original simulation signal and denoised simulation signal

由于在强背景噪声下，仅利用TVD处理仿真信号不能实现较好的故障特征提取，接下来将利用MSB进一步分析图3(d)和图3(f)中的信号。对于图3(d)中的信号，通过比较载波频率各切片的p值，选取288 Hz，658 Hz，3 004 Hz，2 920 Hz和472 Hz处(以p值从大到小进行排序)的切片构成复合切片谱，结果如图5(a)所示；对于图3(f)中的信号，选取472 Hz，288 Hz，658 Hz，2 920 Hz和2 776 Hz处的切片构成复合切片谱，结果如图5(b)所示。

图5中能够明显识别故障特征频率及其倍频成分，白噪声的干扰得到了明显的抑制，且故障特征提取效果要明显优于图4(d)和图4(f)的结果。

图5 本文算法处理后的结果Fig.5 Results of the proposed method

为了进行对比，利用文献[11]提出的Fast Kurtogram算法分析仿真信号2，结果如图6所示。快速峭度图确定的最佳分析频带中心频率为16 213 Hz，带宽为341 Hz，其平方包络谱如图6(b)所示，从图6(b)可知，不能明显识别故障特征频率及其倍频成分，说明对于该仿真信号，快速峭度图方法失效。由此可见，本文的算法在抑制白噪声，提取弱故障冲击方面具有一定的优势。

图6 Fast Kurtogram算法分析仿真信号2的结果Fig.6 Analysis results of simulation signal 2 by Fast Kurtogram

5 实测信号分析

为了验证所提出的方法在滚动轴承内圈和外圈故障特征提取中的有效性，使用实验室滚动轴承故障模拟平台(见图7)进行轴承故障模拟。测点选取在机匣外侧以增加振动信号的复杂程度，测量方向为径向。

图7 实验平台Fig.7 The test rig

5.1 轴承内圈故障分析

轴承内圈故障信号测量时所用的轴承为6 010，该轴承的内径为50 mm，外径为80 mm，滚动体直径为9 mm，滚动体个数为13，接触角α=0°。测试时的转速为3 000 r/min，采样频率为32 768 Hz，通过计算可得故障特征频率的理论值为370.01 Hz。

为了提高分析效果，先对测量得到的振动信号进行标准化处理。图8(a)为标准化处理后，轴承内圈故障信号的时域波形，且从图8(b)的频谱中无法识别轴承故障特征频率及其倍频成分。

图8 轴承内圈故障信号Fig.8 The inner ring fault signal

利用本文的算法处理该内圈故障信号，由图9(a)中的λ值变化曲线，可以确定TVD处理时的参数λ为1.6，得到的降噪后时域波形如图9(b)所示。从图9(b)可知，轴承故障产生的周期性冲击得到了一定的体现，之后利用MSB分析降噪后的故障信号，比较各载波频率切片的p值，选取2 154 Hz，5 292 Hz，3 990 Hz，4 040 Hz和2 522 Hz(以p值从大到小进行排序)处的切片构成复合切片谱，得到的结果如图9(c)所示。

图9(c)中轴承的转频50 Hz及其倍频，轴承内圈故障特征频率为368 Hz(与理论值相接近)以及其倍频736 Hz，1 104 Hz能够得到明显识别，围绕内圈故障特征频率及其倍频，间隔为转频的调制边频带也能得到体现，噪声对分析的干扰得到明显抑制。通过与故障特征频率的理论值相对比，可以判断该轴承存在内圈故障。

图9 本文算法处理轴承内圈故障信号的结果Fig.9 Results of inner ring fault signal by proposed method

为了进行对比，与仿真分析类似，利用文献[11]提出的Fast Kurtogram算法分析该故障信号，结果如图10所示。快速峭度图确定的最佳分析频带中心频率为341 Hz，带宽为682 Hz，其平方包络谱如图10(b)所示。从图10(b)可知，不能明显识别故障特征频率及其倍频成分，快速峭度图方法失效。由此可见，本文的算法在提取弱故障冲击方面具有一定的优势。

图10 Fast Kurtogram算法处理轴承内圈故障信号的结果Fig.10 Results of inner ring fault signal by Fast Kurtogram

5.2 轴承外圈故障分析

轴承外圈故障信号测量时所用的轴承为7 010 C，该轴承的内径为50 mm，外径为80 mm，滚动体直径为8.7 mm，滚动体个数为19，接触角a=15°。测试时的转速为3 000 r/min，采样频率为32 768 Hz，通过计算可得故障特征频率的理论值为413.6 Hz。图11(a)为标准化处理后，轴承外圈故障信号的时域波形，受背景噪声的干扰，周期性故障冲击不能得到识别，同时，图11(b)的频谱中，轴承故障特征频率及其倍频成分也没有得到体现。

图11 轴承外圈故障信号Fig.11 The outer ring fault signal

利用本文的算法处理该外圈故障信号，由图12(a)中的变化曲线可得TVD处理时选取的参数λ为2.5，得到的降噪后时域波形如图12(b)所示。相比原始故障信号，从图12(b)可知，背景噪声得到了明显减少，之后利用MSB分析降噪后的故障信号，选取载波频率2 502 Hz，4 168 Hz，3 334 Hz，4 584 Hz和2 918 Hz处的切片构成复合切片谱，得到的结果如图12(c)所示。

图12(c)中轴承外圈故障特征频率416 Hz(与理论值相接近)以及其倍频832 Hz，1 248 Hz能够得到明显识别，噪声的干扰得到明显抑制。通过与故障特征频率的理论值相对比，可以判断该轴承存在外圈故障。

图12 本文算法处理轴承外圈故障信号的结果Fig.12 Results of outer ring fault signal by proposed method

进行对比，图13为Fast Kurtogram算法分析该故障信号的结果，快速峭度图确定的最佳分析频带中心频率为14 336 Hz，带宽为1 365 Hz，其平方包络谱如图13(b)所示，由图13(b)可知，外圈故障特征频率416 Hz及其三倍频1 248 Hz的幅值十分突出，但受到噪声干扰，故障特征频率的二倍频832 Hz能得到识别但不明显，分析效果欠佳，通过比较，体现了本文算法的有效性。

图13 Fast Kurtogram算法处理轴承外圈故障信号的结果Fig.13 Results of outer ring fault signal by Fast Kurtogram

6 结 论

通过仿真和实验室实测信号验证表明，将TVD和MSB结合进行轴承故障特征提取是可行的。本文得到的结论主要有：

(1)TVD能够有效突出信号中的冲击成分，减少噪声干扰，但受到信号本身特点的制约与强背景噪声的影响，TVD的效果并不是十分理想。通过MSB进行解调可以有效抑制信号中的随机噪声干扰。将MSB与TVD相结合，可以利用两种方法的优势，提高故障特征提取的有效性。相比于经典的快速峭度图方法，本文提出的算法在抑制背景噪声，提取弱故障冲击特征方面更为有效。

(2)针对TVD和MSB解调时的不确定性，利用频域相关峭度选取TVD中的参数λ，通过特征频率成分占比p选取最佳切片位置，确保了结果的准确性。

(3)本文所提出的方法也存在一定的缺陷，TVD作为一种基于信号稀疏表示的方法，其分析结果会受到信号自身特点的影响，因此，如何根据信号自身特点，选取更为合适的稀疏表示方法，获得更理想的降噪效果，值得进一步研究。