基于DTCWPT和相关峭度图的滚动轴承故障诊断

Fault Diagnosis of Rolling Bearings Based on DTCWPT and Correlation Kurtogram

摘要: 本文提出一种基于双树复小波包变换和相关峭度图的滚动轴承故障诊断新方法，解决了在故障诊断领域广泛应用的以峭度为指标的传统快速峭度图和使用矩形窗对峭度作平均化处理的子带平均峭度图易受偶发性冲击干扰、特征提取准确性降低的问题。所提方法利用了DTCWPT的近似平移不变性、频谱泄露程度低、频带排列有序优点的同时在高频段和低频段细分频带并逐层分解分解信号，然后对滤波信号进行上采样和卷积操作重构为与原信号长度一致的重构信号，结合对周期性特征敏感、抑制偶发性冲击、抗噪能力强的相关峭度为子带选择指标，构造相关峭度图，从中选择最优滤波信号，最后对其进行平方包络谱验证滤波信号的故障成分提取情况，实现滚动轴承的故障诊断。与快速峭度图、子带平均峭度图、基于小波包变换的相关峭度图对比，所提方法在信号分解和指标选择两方面进行改进。通过仿真信号和实验数据进行比较，所提方法可有效避免偶发性脉冲的影响，在表现复杂背景下提取滚动轴承的周期性故障特征成分更加准确，所提方法的鲁棒性更好，工程应用价值更高。

Abstract: This paper proposes a new rolling bearing fault diagnosis method based on Dual-Tree Complex Wavelet Packet Transform (DTCWPT) and correlation kurtogram. It solves the problems that traditional fast kurtogram (using kurtosis as indicator) and subband average kurtogram (using rectangular window for kurtosis averaging)—both widely used in fault diagnosis—are easily disturbed by occasional shocks and have reduced feature extraction accuracy. The proposed method uses DTCWPT’s advantages (approximate shift invariance, low spectral leakage, ordered frequency bands) to subdivide and decompose signals in both high and low frequency ranges. It then upsamples and convolves the filtered signal to reconstruct one with the same length as the original. Using correlation kurtogram (sensitive to periodic features, suppresses occasional shocks, strong anti-noise ability) to select the optimal filtered signal, it finally verifies fault component extraction via squared envelope spectrum to achieve diagnosis. Compared with fast, subband averaged kurtogram, and wavelet packet decomposition correlated kurtogram, the method improves signal decomposition and indicator selection. Tests with simulated and experimental data show that it effectively avoids occasional pulses, extracts periodic fault features more accurately under complex backgrounds, and has better robustness and higher engineering value.

1. 引言

滚动轴承是机械设备中支撑系统和传递动力的核心部件，其工作状态直接影响设备加工过程的安全、稳定运行。在长时间的外部载荷和材料疲劳的作用下，滚动轴承不可避免地会出现剥落、磨损、裂纹等缺陷故障形式。出现局部故障时，会产生带有共振调制效应的具有特定频率的周期性脉冲。在加工过程中不可避免会因外部撞击或电磁干扰产生偶发性冲击，其强度通常大于故障成分，严重影响故障提取。在偶发性冲击的影响下精准地提取信号故障特征，是故障诊断领域的主要挑战。

为在复杂的背景噪声下确定故障特征成分所在频段，Antoni在频域峭度的基础上结合非静态过程的Wold-Gramér分解定义谱峭度，检测信号中瞬态脉冲，并提出了快速峭度图(Fast Kurtogram, FK)自适应在不同中心频率与带宽的窄带信号，选择最优频带[1]。该方法通过1/3-二叉树结构的多分辨滤波器组逐步将信号分解到频率和频率分辨率的二维平面上生成不同中心频率和带宽的窄带信号，然后计算各滤波信号复包络的峭度值得到峭度图，最后从中选取最优频带进行故障特征的识别。FK在轴承故障诊断领域中得到广泛应用，但是随着现代机械设备的组成结构、动力传递系统逐渐趋于复杂化、高度集成化，使外部冲击、电磁干扰及背景噪声情况加重，而峭度容易受到这些干扰影响，严重影响FK的频带选择[2]。

针对于此，许多研究集中于改进评价指标使其更稳定地衡量信号周期性与冲击性。Wang提出了子带平均峭度(Subband Averaging Kurtogram, SAK)，通过矩形窗截取信号计算平均峭度指导子带选择[3]。相关峭度由McDonald提出，关注特定频率下重复出现的脉冲成分，几乎不受随机冲击影响[4]。Liu将相关峭度与经验小波变换(Empirical Wavelet Transform, EWT)结合，提出了自适应相关峭度图方法[5]，但EWT自适应分割边界的方式对信号移位敏感，小幅相位的变化可能使分解结果发生剧变，限制了构造相关峭度图的准确性；Zhang提出的小波包分解相关峭度图(Wavelet Packet Decomposition Correlated Kurtogram, WPDCK)方法通过小波包分解(Wavelet Packet Decomposition, WPD)结合相关峭度确定最优频带[6]，虽然WPD通过对高频与低频子带同时进行分解，实现了全频谱均匀分解，但仅能提供实部方向的选择，在非平稳信号方面分解效果不佳。双树复小波变换利用一对近似满足希尔伯特变换关系的实部与虚部滤波器组，提供了更好的方向选择性与近似平移不变性，有效避免频率混叠，解决了传统小波变换移位敏感、频谱泄露的问题。在此基础上，DTCWPT (Dual-Tree Complex Wavelet Packet Transform, DTCWPT)额外补充了实部与虚部高频子带的分解，减少了信息缺失，提供更为精细的频率分辨率。

为解决偶发性脉冲干扰FK、SAK频带选择，进而导致故障诊断失效的问题，本文提出了基于DTCWPT的相关峭度图方法。利用近似平移不变性和频谱泄露程度低的DTCWPT分解信号，以相关峭度为指标计算还原成与输入信号长度一致的重构信号进行最优频带的选择，最后通过包络解调完成故障诊断。通过模拟信号与实验结果与FK、SAK、WPDCK对比，验证了所提方法抗干扰性强、鲁棒性更好。

2. FK、SAK、CK简述

振动信号$X\left(n\right)$ 可以表示为一个非平稳的瞬态故障成分$Y\left(n\right)$ 和一个与其无关的平稳高斯噪声$B\left(n\right)$ 的混合形式：$X\left(n\right)=Y\left(n\right)+B\left(n\right)$ 。

峭度反映了信号数值概率分布的尖峭程度，是检测信号中脉冲特征的重要指标，$X\left(n\right)$ 的峭度表示为：

$\text{Kurtosis}=\frac{{\displaystyle {\sum }_{i=1}^N{\left(x_i-{\mu }_x\right)}^4/N}}{{\left({\displaystyle {\sum }_{i=1}^N{\left(x_i-{\mu }_x\right)}^2}/N\right)}^2}-3$ (1)

式中：$N$ 和${\mu }_x$ 为信号$x\left(n\right)$ 的长度和均值。

谱峭度(SK)用于衡量信号在特定频率下非高斯特性的统计指标[7]，定义如下：

$S{K}_x\left(f\right)=\frac{\langle {\left|H\left(n,f\right)\right|}^4\rangle }{{\langle {\left|H\left(n,f\right)\right|}^2\rangle }^2}-2$ (2)

式中：$\langle \cdot \rangle$ 表示平均运算符，${\left|H\left(n,f\right)\right|}^4$ 和${\left|H\left(n,f\right)\right|}^2$ 分别是$x\left(n\right)$ 的四阶和二阶瞬时矩。

FK通过1/3-二叉树结构的多分辨滤波器组分解信号，生成不同中心频率和带宽的窄带信号，输出复包络谱峭度最大的信号。但FK计算峭度采用的是仅通过降采样得到的缩短长度的子带信号，虽然显著减少了计算量、提高计算效率，但存在关键信息缺失的风险，如峭度与信号峰值和尾部特征密切相关，但降采样后可能会使某些瞬态特征的时间细节信息丢失，这增大了峭度指标计算误差进而降低FK准确性。

SAK通过DTCWPT分解不同层级的滤波信号，利用滑动矩形窗口按预设分段数求解信号的平均峭度。减小了随机性干扰误导子带选择的影响，但这在本质上只是峭度的平均化处理，并没有从根源解决峭度指标的抗噪性差、易受干扰的问题。

相关峭度(CK)最初的提出是为解决自适应FIR滤波器的设计易受高强度、多个引起不同共振频率的随机冲击而影响滤波性能的问题，而CK同时考虑到旋转机械故障特征的周期性和冲击性且几乎不受随机脉冲和谐波信号的干扰，出色地鼓励滤波器围绕周期提取故障信息，相关峭度针对预设周期$T$ ，同时考虑了故障信号周期性和冲击性的特点，定义如下：

$C{K}_M\left(T_s\right)=\frac{{\displaystyle {\sum }_{i=1}^N{\left({\displaystyle {\prod }_{m=0}^M{x}_{i-m{T}_s}}\right)}^2}}{{\left({\displaystyle {\sum }_{i=1}^N{x}_i^2}\right)}^{M+1}}\left(M\in {\mathbb{Z}}^{+}\right)$ (3)

式中：$M$ 为移位阶数，$N$ 为信号$x$ 的信号长度，$T_s$ 是信号周期长度，$T_s=T\cdot f_s$ 。

作为稀疏性指标，峭度对信号中随机冲击成分十分敏感。当信号包含周期性冲击时，峭度值会显著减小。与之相比，相关峭度抗噪性更强，抑制随机脉冲干扰，是关注信号重复性与脉冲性的指标。

3. 基于DTCWPT的相关峭度图

3.1. DTCWPT

DTCWT是由一对彼此近似形成Hibert变换对的小波滤波器组成，其最大特点是几乎不存在频谱混叠和近似平移不变性。WPD仅能支持实部的低频与高频分解，与之相比，DTCWT提供了更好的方向选择性与近似平移不变性，但由于仅进行低频分解，提供的频率分辨率不足。针对于此，Selesnick以不同低通与高通滤波器按实树与虚树的小波包平行排布[8]滤波结构提出了DTCWPT方法，在保留了DTCWT优势的基础上利用小波包变换理论，通过分解高频子带，获得实树与虚树的高频输出，提供了更为精细的频率分辨率，更适合复杂信号的分析。三种方法构成子带，分解信号的方式如图1所示。

Figure 1. Calculation flowchart of filtered signal and reconstructed signal. (a) WPD, (b) DTCWT, (c) DTCWPT

图1. 滤波信号与重构信号计算流程图。(a) WPD，(b) DTCWT，(c) DTCWPT

信号$x\left(t\right)$ 与小波滤波器组的实部和虚部卷积得到小波系数$d_k\left(m\right)$ 和尺度系数$c_k\left(m\right)$ ：

$\{\begin{array}{l}{d}_k^{*}\left(m\right)=2^{k/2}{\displaystyle {\int }_{-\infty }^{+\infty }X\left(t\right){\psi }^{*}}\left(2^{k}t-m\right)\text{d}t\\ c_k^{*}\left(m\right)=2^{k/2}{\displaystyle {\int }_{-\infty }^{+\infty }X\left(t\right){\varphi }^{*}}\left(2^{k}t-m\right)\text{d}t\end{array}$ (4)

式中：k是分解尺度，$\psi \left(t\right)$ 与$\varphi \left(t\right)$ 是双树复小波包的小波与尺度函数，“$*$ ”是实部或虚部。

重构信号$\tilde{x}\left(t\right)={\displaystyle {\sum }_{k=1}^n{d}_k\left(t\right)}+{\displaystyle {\sum }_{k=1}^n{c}_k\left(t\right)}$ 。

$\{\begin{array}{l}{d}_k\left(t\right)=2^{\left(k-1\right)/2}{\displaystyle {\sum }_{i=1}^n{d}_k^{*}\left(i\right){\psi }^{*}\left(2^{k}t-i\right)}\\ c_k\left(t\right)=2^{\left(k-1\right)/2}{\displaystyle {\sum }_{i=1}^n{c}_k^{*}\left(i\right){\varphi }^{*}\left(2^{k}t-i\right)}\end{array}$ (5)

式中：k是分解尺度，$\psi \left(t\right)$ 与$\varphi \left(t\right)$ 是双树复小波包的小波与尺度函数，“$*$ ”是实部或虚部。

逐层下采样滤波使信号长度迅速减少，存在丢失关键信息的隐患。对固有频率为50 Hz的周期性脉冲按周期长度截取信号，峭度(K)、平均峭度(SAK)、相关峭度(CK)三种指标随信号长度的变化情况如图2所示。

Figure 2. Relationship diagram of indicator value vs. signal length

图2. 指标值–信号长度关系图

四种指标均随信号长度增加逐渐收敛。峭度变化平稳、收敛速度快，而相关峭度与平均峭度收敛缓慢，变化波动大。因此，采用下采样长短不一的滤波信号无法公平对比，故计算解卷积后上采样的重构信号，保证计算精准。

3.2. 所提方法流程

所提方法规避了FK、SAK指标易受随机脉冲干扰的问题，采用相关峭度评价重构信号的周期性与冲击性。算法采用重构还原后长度一致的重构信号计算指标，流程如下：

1) 加载振动信号$x\left(t\right)$ ，根据经验预设置分解层数Level及$T$ 、$M$ 、$f_s$ 等参数。

2) 利用DTCWPT将信号逐层分解到不同频带，进行降采样，生成滤波信号。

3) 对滤波信号上采样后重构为与原信号长度一致的重构信号，根据频带依序排列。

4) 计算重构信号相关峭度构建相关峭度图，输出最大CK值的重构信号包络分析。

4. 仿真验证

模拟轴承故障信号$x\left(t\right)$ 包含信噪比为−10的高斯白噪声$g\left(t\right)$ ，谐波振动$h\left(t\right)$ ，模拟轴承外圈故障引起的周期性脉冲$p\left(t\right)$ ，模拟外部敲击或电磁干扰引起的偶发性脉冲$s\left(t\right)$ ，设置采样频率为20,000 Hz，采样长度为20,000。$x\left(t\right)$ 由式(6)表述，各参数设置参考表1，$x\left(t\right)$ 及组成部分的时域波形如图3所示。根据图3(e)，观察$x\left(t\right)$ 的时域图，轴承外圈故障缺陷引起的周期性故障信号$p\left(t\right)$ 被淹没在复杂的背景噪声中，难以直接观测出有效成分。

Table 1. Parameters of simulated signal $x\left(t\right)$

表1. 仿真信号$x\left(t\right)$ 的参数

$f_P$	$f_{NP}$	$f_{NS}$	${\eta }_P$	${\eta }_S$	$A_{H1}$	$A_{H2}$
22	2000	5100	1000	1500	0.05	0.025

Figure 3. Time-domain waveforms of each component of the simulated signal. (a) p(t), (b) s(t), (c) h(t), (d) g(t), (e) x(t)

图3. 仿真信号各分量时域波形。(a) p(t), (b) s(t), (c) h(t), (d) g(t), (e) x(t)

$\{\begin{array}{l}x\left(t\right)=p\left(t\right)+s\left(t\right)+h\left(t\right)+g\left(t\right)\\ p\left(t\right)={\displaystyle {\sum }_{i=1}^{f_P}{\text{e}}^{-{\eta }_P\left(t-i/f_P\right)}\cos \left(2\pi f_{NP}\left(t-i/f_P\right)\right)}\\ s\left(t\right)=2{\displaystyle {\sum }_{i=1}^2{\text{e}}^{-{\eta }_S\left(t-t_S\right)}\cos \left(2\pi f_{NS}\left(t-j{t}_S\right)\right)}\\ h\left(t\right)=A_{H1}\cos \left(14\pi t\right)+A_{H2}\cos \left(25\pi t\right)\end{array}$ (6)

式中：$f_P$ 是故障频率，$f_N$ 是共振频率，$\eta$ 是阻尼系数，$t_S$ 是脉冲间隔，$A_H$ 是谐波振幅。

$x\left(t\right)$ 各组成部分的K、SAK、CK值如表2所示。$s\left(t\right)$ 的K、SAK指标显著高于其他信号，说明SAK没有在根源消除$s\left(t\right)$ 引起峭度值激增的影响。$p\left(t\right)$ 的CK值在数量级上优于其他信号而$s\left(t\right)$ 的CK值锐减，说明相关峭度抑制无循环性的偶发性脉冲。

Table 2. Three indicator values of common signals

表2. 常见信号的三个指标值

	$p\left(t\right)$	$s\left(t\right)$	$h\left(t\right)$	$g\left(t\right)$
K	66.7	1444.3	−0.8	−0.1
SAK	68.8	672.3	−0.7	−0.1
CK	$3.3\times {10}^{-3}$	$6.2\times {10}^{-62}$	$5.7\times {10}^{-5}$	$4.7\times {10}^{-5}$

计算$p\left(t\right)$ 、$p\left(t\right)+s\left(t\right)$ 、$p\left(t\right)+h\left(t\right)+g\left(t\right)$ 、$p\left(t\right)+h\left(t\right)+g\left(t\right)+s\left(t\right)$ 的三种指标，对比验证$s\left(t\right)$ 对K、SAK、CK指标的影响与相关峭度对$s\left(t\right)$ 的抑制效果，如表3所示。$p\left(t\right)+h\left(t\right)+g\left(t\right)$ 与$p\left(t\right)$ 在添加$s\left(t\right)$ 后，K与SAK值均显著高于原信号，说明SAK仅减轻而非根除偶发性脉冲引起峭度指标增加的问题；而CK在添加噪声后呈减小趋势，意味着添加$s\left(t\right)$ 的信号稀疏性与周期性变差。故CK评价故障信号的周期性准确。

Table 3. Indicator values of the signal after adding $s\left(t\right)$

表3. 添加s(t)后的信号指标值

信号成分	K	SAK	CK
$p\left(t\right)$	66.7	68.8	$3.3\times {10}^{-3}$
$p\left(t\right)+s\left(t\right)$	108.4	89.4	$2.0\times {10}^{-3}$
$p\left(t\right)+h\left(t\right)+g\left(t\right)$	0.6	0.6	$7.8\times {10}^{-5}$
$p\left(t\right)+h\left(t\right)+g\left(t\right)+s\left(t\right)$	1.7	1.5	$7.5\times {10}^{-5}$

首先采用所提方法分解$x\left(t\right)$ ，提取故障成分，构造相关峭度图，根据经验预设Level = 5、T = 1/22 s、移位阶数M = 1。根据图4(a)，所选最优重构信号CK值约为$2.1\times {10}^{-4}$ ，显著高于其他区域。信号频段[1250, 2500] Hz、中心频率1875 Hz，与准确值基本相符。图4(b)是重构信号的时域波形，可以清晰地观察到脉冲成分呈周期性排列，结合观察图4(c)平方包络谱(SES)，周期性故障及其倍频清晰可见。

形成鲜明对比的是：FK选择[5000, 5312.5] Hz的窄带信号，与准确值相悖，而与偶发性噪声成分$s\left(t\right)$ 的共振频率$f_{NS}$ 基本一致，这说明FK被$s\left(t\right)$ 误导而错误地提取故障特征。观察图4(d)与图4(e)，滤波信号几乎被噪声淹没、没有周期性特征，这意味着FK在偶发性冲击的背景下滤波效果变差、不能选择正确的频段。将SAK的平均分段数设为8，根据图4(f)，滤波信号脉冲位置与X(t)中偶发性脉冲发生位置基本一致，证明SAK受到影响过滤出主要成分为偶发性冲击$s\left(t\right)$ 的滤波信号而非理想的周期性特征，分解没有达到预期效果。

对比图4(a)与图5(a)，在参数设置相同且均进行信号重构的情况下，所提方法重构信号CK值显著高于WPDCK组，各子带间指标差距鲜明，利用WPD构建的相关峭度图各子带信号相关峭度值差距较小，难以满足复杂信号的分解需求。

Figure 4. Simulation tests result. (a) Correlation kurtogram of the proposed method, (b) Time-domain waveform of the reconstructed signal by the proposed method, (c) SES of the reconstructed signal, (d) Time-domain waveform of the FK stimulus signal, (e) FK-SES, (f) Time-domain waveform of the SAK stimulus signal, (g) SAK-SES

图4. 仿真测试结果。(a) 所提方法的相关峭度图，(b) 所提方法重构信号时域波形，(c) 重构信号SES，(d) FK滤波信号时域波形，(e) FK-SES，(f) SAK滤波信号时域波形，(g) SAK-SES

Figure 5. WPDCK simulation tests result. (a) Correlation kurtogram of WPDCK, (b) Time-domain waveform of WPDCK reconstructed signal, (c) WPDCK-SES

图5. WPDCK仿真测试结果。(a) WPDCK的相关峭度图，(b) WPDCK重构信号时域波形，(c) WPDCK-SES

5. 实验结果与分析

为进一步验证所提方法的有效性与实用性，采用美国凯斯西储大学滚动轴承实验数据对所提方法与FK、SAK进行对比验证。在SKF6205型号滚动轴承外圈预先生成直径为0.1778 mm的单点缺陷，$f_s$ = 12 kHz、电机转速1750 r/min，理论故障频率为104.63 Hz。

根据经验预设Level = 5、T = 1/105 s、移位阶数M = 3。根据图6(a)，在相关峭度图上，所提方法判断最优子带频段[2625, 3000] Hz，定位准确，频率分布随层数的增加不断精确频率范围，分解达到第5层时，由于CK指标的优越能力，避免了过度分解与带宽过窄引起的信息缺失的问题。观察时域图及其包络：故障成分按周期分布，接近外圈故障频率104.63 Hz及其倍频(约至7倍)，噪声比例小。据此判断滚动轴承发生故障，该方法可以有效提取故障特征。

由于实验环境下相较于模拟仿真存在的背景噪声及干扰成分更为复杂，峭度与子带平均峭度指标易受高频噪声的误导，FK与SAK不约而同地将故障频段锁定在高频区域，与实际情况不符。根据图6(d)~(g)，FK提取了掺杂明显噪声的微弱周期性故障，故障特征不明显，至4倍频处消失。SAK被背景噪声误导，未滤波出有效成分，提取周期性故障特征失败。

观察图7，由于高频噪声的干扰，WPDCK的相关峭度图中各子带分布不清晰，故障频段的相关峭度难以区分于干扰区域；重构信号的时域波形周期循环特征不明显，结合包络图中周期故障频率约至2倍频处难以直接识别的表现，验证了WPD对背景复杂的非平稳信号分解能力存在局限性从而影响构建相关峭度图的缺陷。

Figure 6. Experimental result. (a) Correlation kurtogram of the proposed method, (b) Time-domain waveform of the reconstructed signal using the proposed method, (c) SES of the reconstructed signal, (d) Time-domain waveform of the FK-filtered signal, (e) FK-SES, (f) Time-domain waveform of the SAK-filtered signal, (g) SAK-SES

图6. 实验结果。(a) 所提方法的相关峭度图，(b) 所提方法重构信号时域波形，(c) 重构信号SES，(d) FK滤波信号时域波形，(e) FK-SES，(f) SAK滤波信号时域波形，(g) SAK-SES

Figure 7. Experimental result of WPDCK. (a) Correlation kurtogram of WPDCK, (b) Time-domain waveform of WPDCK reconstructed signal, (c) WPDCK-SES

图7. WPDCK实验结果。(a) WPDCK的相关峭度图 (b) WPDCK重构信号时域波形 (c) WPDCK-SES

对比四种方法，所提方法在提取故障周期性特征方面鲁棒性更好，具有良好的有效性和实用性。

6. 结论

本文提出了一种基于DTCWPT和相关峭度图的滚动轴承故障诊断方法。引入抑制偶发性冲击的相关峭度对重构成与输入信号长度一致的滤波信号选择最优子带，解决了FK、SAK易受尖锐偶发性脉冲干扰、错误选择频带的问题。通过仿真信号、实验信号进行故障诊断并与FK、SAK、WPDCK对比分析，所提方法更为准确地提取故障的周期性特征，识别效果更好，验证了所提方法的有效性和鲁棒性。

基金项目

航空科学基金(20183368004)项目资助。

NOTES

^*通讯作者。