作为旋转机械设备中关键零部件之一，滚动轴承的运行健康状态关系到整个设备的运转情况. 对轴承的运行状态做出准确监测和故障诊断具有重要的理论和实际意义^[1-2]. 信号处理与分析是机械设备健康状态评价的重要依据. 傅里叶变换和功率谱估计方法都是基于平稳信号分析理论，不适用于表征非平稳信号在某一时刻的频率成分分布情况^[3]. 对于非线性、非平稳信号，时频分析可以很好地保留时间和频率信息，描述出任意时刻的频率成分. 各种特征提取及诊断方法不断地发展，其中短时傅里叶变换、小波分析、经验模态分解(empirical mode decomposition, EMD)等时频分析方法被广泛应用于非线性、非平稳信号分析领域，并且取得了一定的成果. 王涛等^[4]比较分析了多种时频分析方法，从时间分辨率、频率分辨率等方面阐述各时频方法的优缺点. Obuchowski等^[5]改进了短时傅里叶方法，利用一个在时间上可滑移的时间窗进行傅里叶变换，在时域和频域上都具有较好的局部性. 但是短时傅里叶变换的窗函数固定，无法适用于分析时频分辨率有变化的非平稳信号. 王宏超等^[6]基于二阶循环统计量对谱相关或谱相关密度分析方法加以改进，提出一种提取调制频率的时频分析方法，具有强抗噪能力. 李宏坤等^[7]结合EMD和时频分析，提出基于时频分析的欠定盲源分离方法进行设备微弱特征提取，具有较好的分离效果.

熵被广泛引入到机械设备状态监测与故障诊断领域. 林京等^[8]将信息熵引入到信号时频分析，按时频块划分能量，并用于齿轮的故障诊断. 申弢等^[9]提出小波空间熵概念，在时频联合分布上按小波尺度分段能量，用作评价振动状态的特征指标. 于德介等^[10]提出基于EMD的时频熵，成功用作分析齿轮的工作状态和缺陷类型的指标. 王小玲等^[11]结合幅值熵和时频分析方法，表征全频域随时间变化情况，并用作故障类型的分类. 刘韬等^[12]在王小玲研究结果的基础上，基于频带熵设计自适应滤波器，实现从噪声中提取轴承故障特征频率. 刘学等^[13]提出自适应多尺度时频熵，并运用到遥测领域，实现自适应细致划分时频信息，解决了传统时频熵方法不能敏锐反应信息量变化、局部时频熵依赖人为经验的问题.

本研究提出一种基于时频分析的频带多尺度复合模糊熵本研究的故障特征提取方法. 该方法利用多尺度复合模糊熵对微弱信号潜在信息的挖掘能力，结合时频分析，辨识出系统固有频率并解调故障特征频率，实现对故障冲击成分的特征提取.

1 频带多尺度复合模糊熵方法 1.1 时频原理

短时傅里叶变换(short-time Fourier transform, STFT)作为一种常规、惯用的时频分析方法，适用于分析非平稳、非线性信号. 为了探究信号的局部频率特性，STFT通过对信号添加一个宽度很窄的时间窗，并且沿时间轴移动进行傅里叶相关变换，计算信号局部性的频率和相位，绘制出频带分量随时间变化的时频分布谱图.

${\rm{STFT}}\{ x(t)\} (\tau ,\omega ) = \int_{ - \infty }^\infty {x(t)\omega (t - \tau ){{\rm e}^{ -{\rm j}\omega t}}{\rm d}t}. $

(1)

式中：j为虚数单位， $\omega (t)$ 为窗函数，通常是以零为中心的汉宁或高斯窗； $x(t)$ 为分析信号. STFT反映 $x(t)$ 在 $\tau $ 时刻、频率 $\omega $ 处的信号成分相对含量.

STFT是一种单一分辨率的时频分析方法^[14]. STFT采用固定的时间窗把分析信号截断成许多小的时间间隔，因此时频分辨率受Heisenberg测不准原理约束，不能同时兼顾时间分辨率和频率分辨率. 时间窗越窄，对应的时域上的时间分辨率越好，但相应的频率分辨能力越差，影响频率分辨效果.

1.2 模糊熵

熵定义为新信息产生率，可以用于描述时间序列的无规则程度和不确定性，并被广泛用于机械设备故障诊断领域^[15]. 模糊熵的主要计算过程如下。

1) 对时间序列 ${ X}= \{ {x_1},{x_2}, \cdots ,{x_N}\} $ ，设定模式维数 $m$ ，定义 ${ X}(i) = [{x_i},{ x_{i + 1}}, \cdots ,{x_{i + m - 1}}],1 \leqslant i \leqslant N - m$ .

2) 定义 ${ X}(i)$ 与 ${ X}(j)$ 之间的距离 ${d_{ij}}$ 为两者对应元素之差中绝对值的最大值，即

$d_{ij}^m = \mathop {\max }\limits_{k = 1, \cdots ,m - 1} (\left| {x(i + k) - x(j + k)} \right|). $

(2)

3) 引入模糊隶属度函数

$A(x)= \left\{ \begin{aligned}& 1, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, x = 0 \text{；} \\ & \exp \, \left[ { - \ln \, \left(2{{\left( {\frac{x}{r}} \right)}^2}\right) }\right], \,\,\,\,\,\,\,\,\,\,\, x >0 .\end{aligned} \right.$

(3)

式中： $r$ 为相似容限参数， $r = R \cdot\rm SD$ (R为对序列噪声水平的估计参数，SD为原始数据的标准差).

向量 ${{{X}}_m}(i)\text{、}{{{X}}_m}(j)$ 间相似度定义为

$A_{ij}^m = \left\{ \begin{aligned}& 1, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{ij}^m = 0 \text{；} \\ & \exp \, \left[ { - \ln \, \left(2{{\left( {\frac{{d_{ij}^m}}{r}} \right)}^2}\right)} \right], \,\,\,\,\,\,\,\,\,\,\, d_{ij}^m >0 .\end{aligned} \right.$

(4)

4) 对于维数 $m$ ，对每一个 $j$ 值统计与其余矢量的互相相似概率：

$C_i^m(r) = \frac{1}{{N - m}}\sum\limits_{j = 1,j \ne i}^{N - m + 1} {A_{ij}^m} . $

(5)

对所有的 $C_i^m(r)$ 取平均，平均值记为

${B^m}(r) = \frac{1}{{N - m + 1}}\sum\limits_{i = 1}^{N - m + 1} {C_i^m}(r) . $

(6)

5) 对于维数 $m + 1$ ，重复1)~4)，得到 ${B^{m + 1}}(r)$ .

6) 最后得到模糊熵的表达式为

${\rm{FEn}}(m,r,N) = \ln {B^m}(r) - \ln {B^{m + 1}}(r).$

(7)

式中： $N$ 为时间序列长度.

1.3 基于方差的多尺度复合模糊熵

郑近德等^[16]基于模糊熵，利用粗粒化获得多尺度时间序列，提出一种多尺度模糊熵，用于描述时间序列在不同时间尺度上的复杂性. 虽然多尺度模糊熵有效克服了模糊熵只能从单一尺度反映时间序列信息的不足，但是对故障特征相似的不同类型故障仍无法准确辨识^[17]. Azami^[18]等基于多尺度模糊熵，引入滑移间隔因子，提出了多尺度复合模糊熵. 在粗粒化序列时，采用方差替代原来的均值计算，提取序列中隐藏的微弱信号. 具体算法如下。

1) 设定相似容限r、尺度因子 $\tau \; (\tau \geqslant 1 )$ 和模式维数 $m$ ，对长度为 $N$ 的时间序列 $X = \{ {x_1},{x_2}, \cdots ,{x_N}\} $ 进行粗粒化，基于粗粒序列的方差建立新的时间序列，并根据间隔因子将粗粒时间序列不断向后滑移，构造出滑移粗粒时间矩阵.

$\left.\begin{aligned}{y_{k,j}}(\tau ,p) = & {\rm{Var}}({x_i})\left| {_{(j - 1)\tau + k}^{j\tau + k - 1}} = \frac{1}{\tau }\sum\limits_{i = (j - 1)\tau + k}^{j\tau + k - 1} {x_i^2 - {{(\overline x )}^2}} \right.;\\ & {1 \leqslant j \leqslant \left\lfloor {\frac{{N - k}}{\tau }} \right\rfloor = c} , {1 \leqslant k \leqslant \left\lfloor {{\frac{\tau }{p}}} \right\rfloor = r}. \end{aligned} \right\}$

(8)

式中：p为滑移间隔因子. 当 $\tau = 1$ ， $p = 0$ 时， ${y_{k,j}}(1,0)$ 即为原始时间序列. 对于非零 $\tau $ 和p，原始时间序列{X_i}被构造成 $\left\lfloor {\tau /p} \right\rfloor $ 维的滑移粗粒时间矩阵(式中 $\left\lfloor {\,\, } \right\rfloor $ 是向下取整) ，其中第 $k$ 行滑移粗粒时间向量被分割成 $(N - k)/\tau $ 段、长度为 $N/\tau $ 的粗粒时间序列 $\{ {y_{k,j}}(\tau ,p)\} $ .

2)对于滑移粗粒时间矩阵中第 $k$ 行滑移粗粒时间向量 $\{ {y_{k,j}}(\tau ,p)\} $ ，计算 $B_k^{^m}(r)$ .

$B_k^{^m}(r) = \frac{1}{{N - m + 1}}\sum\limits_{i = 1}^{N - m + 1} {C_i^m}(r) . $

(9)

则对于 $\left\lfloor {\tau /p} \right\rfloor $ 维的滑移粗粒时间矩阵为

${B^m}(r) = \sum\limits_{k = 1}^{\tau /p} {B_k^{^m}(r)} . $

(10)

3) 对于维数 $m + 1$ ，重复1)~2)，得到 ${B^{m + 1}}(r)$ .

4) 对尺度因子为 $\tau $ 时的滑移粗粒时间矩阵计算模糊熵：

${\rm{FEn}}(\tau ,p) = \ln {B^m}(r) - \ln {B^{m + 1}}(r). $

(11)

5) 最后得到多尺度复合模糊熵的表达式：

$\begin{aligned}& {\rm {RCMFEn}}({\rm{scale}},p,m,r,N)=\\ & \left\{{\rm FEn}(1,p),{\rm FEn}(2,p), \cdots ,{\rm FEn}({\rm{scale}},p)\right\}.\end{aligned} $

(12)

式中： ${\rm{scale}}$ 为尺度因子. 尺度因子 ${\rm{scale}}$ 和滑移间隔因子 $p$ 的选取对多尺度复合模糊熵的影响显著. 模式维数和相似容限的选取也非常重要，根据文献[19, 20]的研究结果，当取 $m = 1$ 或 $ 2$ ， $0.1\rm SD \leqslant r \leqslant 0.25\rm SD$ 时所得模糊熵的统计特性及故障状态区分度更好.

1.4 基于时频分析的频带多尺度复合模糊熵

时频分析用来表征非平稳信号中所含的频率分量随时间的变化特性，但是无法体现信号的频率结构复杂度. 本文提出基于时频分析的频带多尺度复合模糊熵方法，用来度量信号中各频率分量分布均匀情况和随时间变化特性，反映信号能量分布的不确定性和频率复杂度. 图1为该方法的流程示意图.

对原始信号 $\{ {x_i}\} (i = 1,2, \cdots ,N)$ 做短时傅里叶时频变换，得到信号的时频分布矩阵TFA，分别由不同的频带向量 ${{{R}}_{i}}(i = 1,2, \cdots ,M)$ 组成，表征从小到大第 $i$ 个频率分量 ${f_i}$ 随时间变化的时域分布.

${{{\bf {TFA}}}} = \left[ {\begin{array}{*{20}{c}} {r_{1,1}}& \cdots &{r_{1,C}} \\ \vdots & {} & \vdots \\ {r_{M,1}}& \cdots &{r_{M,C}} \end{array}} \right] = {[\begin{array}{*{20}{c}} {{{R}}_1^{\rm{T}}},&{{{R}}_2^{\rm{T}}},& \cdots ,&{{{R}}_M^{\rm{T}}} \end{array}]^{\rm{T}}}.$

(13)

式中： $M$ 为频带点数， $C$ 为窗函数沿时间轴滑动后对应的时间序列长度.

${{\rm{En}}_f} = \left\{{{\rm{En}}_{f1}},{{\rm{En}}_{f2}}, \cdots ,{{\rm{En}}_{fM}}\right\},$

(14)

${{\rm{En}}_f}(k) = \min \,\, ({{\rm{En}}_f}).$

(15)

对于不同频带向量求解多尺度复合模糊熵，得到各个频率分量的时间复杂度 ${{\rm{En}}_f}$ ，揭示了信号在不同频带处沿时间轴的变化情况.

模糊熵是用来定量描述信号复杂度的方法，当某个频率分量变化平稳或规律性地蕴含在信号中，则该频率线对应较小的模糊熵值，反映出信号复杂度较低，信息平缓稳定、分布均匀；当频率分量在时间响应上波动频繁，则会对应较大的模糊熵. 该特性可以运用到设备的故障特征提取领域. 根据共振解调法理论^[21]，振动信号在一段时间内，全频带中最小模糊熵所对应的频率分量，可以认为是系统的固有频率 $f_{\rm{n}}$ . 信号所在频率分量 $f(k)$ 处的时间分布，即频带向量 ${{{R}}_k}$ ，认为是信号在固有频率 $f_{\rm{n}}$ 及谐波调制下随时间的变化响应，包含有丰富的故障信息和特征频率.

频带多尺度复合模糊熵基于时频联合分析，实现了定量表征信号的动态变化，反映出各个频率随时间分布的复杂度. 同时找出固有频率对应的频率分量 $f(k)$ ，通过对共振频带向量 ${{{R}}_k}$ 的重构分析，提取故障特征.

1.5 自适应带通滤波器

利用基于时频分析的频带多尺度复合模糊熵，设计自适应带基于通滤波器，可以实现轴承故障冲击的特征识别，有效提高分辨能力和信号信噪比.

自适应带通滤波器的设计过程如下：首先引入有限冲击响应(FIR) 数字滤波器，然后通过找到匹配的 $(f_{\rm{c}},\Delta f)$ 系数和滤波阶数来设计FIR滤波器 $G(f)$ ，其中 $f_{\rm{c}}$ 为中心通过频率， $\Delta f$ 为 $G(f)$ 的滤波带宽. 设计参数公式为

$\left. \begin{aligned}& {f_0} = f_{\rm{s}}/N_{\rm{w}},\\ & {{\rm{En}}_f}(k) = \min \, ({{\rm{En}}_f}(i)),i = 1,2 ,\cdots ,M,\\ & f_{\rm{c}} = k {f_0},\\ & \Delta f = 1.5{f_0}.\end{aligned} \right\}$

(16)

式中： $f_{\rm{s}}$ 为采样频率， $N_{\rm{w}}$ 为STFT窗函数长度.

FIR滤波器是通过对已知脉冲形状进行采样并利用这些样本以相反的阶数作为滤波器系数，来设计自适应滤波器的脉冲响应. 结合时频分析和频带多尺度复合模糊熵的特性，可实现既保证滤波器中心通过频率的识别准确率，又加强对故障冲击特征的提取、分辨能力.

2 实验验证分析 2.1 仿真分析

根据滚动轴承故障模型对滚动轴承外圈故障进行模拟仿真，来验证基于时频的频带多尺度复合模糊熵方法的有效性.

$\left. \begin{aligned}& x(t) = \sum\limits_{i = 1}^M {{A_i} s(t - iT - \tau _{i}) + n(t),} \\ & {A_i} = {A_0}\cos \,\, (2\text{π} {f_{\rm{m}}}t + {\varphi_{\rm{A}}}) + {C_{\rm{A}}},\\ & s(t) = {{\rm e}^{ - Bt}} \cos \,\, (2\text{π} {f_{\rm{n}}}t + {\varphi_{\rm{w}}}).\end{aligned} \right\}$

(17)

式中： ${\tau _i}$ 为第 $i$ 次冲击相对于平均轴T的微小波动. n(t)为代表背景噪声的随机信号，f_m为轴承信号的故障特征频率， ${\rm}{\varphi_{\rm{A}}}$ 、 ${\rm}{\varphi_{\rm{w}}}$ 和C_A为任意常数，f_n代表系统固有频率，B为系统阻尼系数。

设置对应仿真信号系统的采样频率 $f_{\rm{s}} $ =25 600 Hz，转频 $f_{\rm{r}}$ =12 Hz，外圈通过频率 $f_{\rm{o}}$ =40 Hz，系统固有频率 $f_{\rm{n}}$ =3 700 Hz. 信号包含信噪比为–15 dB的高斯白噪声.图2分别是仿真信号的波形图和包络谱图，时域图中噪声几乎将故障特征淹没，包络谱中充斥干扰噪声，无法准确识别特征谱线. 图3为对应的时频谱图.

针对图2所示的轴承外圈故障仿真信号，利用频带多尺度复合模糊熵进行故障特征识别与提取，如图4所示. 图4(a)是基于频带的单尺度模糊熵，图4(b)是对应的多尺度复合模糊熵，其中不同尺度对应的熵值线在f=3 700 Hz处最低，可以认为此处是系统的固有主频率. 从时频分析的时频矩阵TFA中，提取固有频率对应的时间分布系列，从图5(a)中可以看出时频分析后的重构信号对噪声取得较好的抑制，从包络谱图中可以轻易辨识出系统的故障特征频率为40 Hz.

2.2 实例分析

为了验证频带多尺度复合模糊熵方法的有效性。搭建人工故障滚动轴承试验台，如图 6所示. 试验台主要包括减速箱、伺服电机和磁粉制动器. 实验研究对象为30304型的圆锥滚子轴承，通过电火花加工、人工制造轴承微小故障来模拟轴承的真实失效. 图7是外圈轴承故障. 设置系统采样频率25.6 kHz，数据采样长度为6 400个数据点.

图8为采集的30304型轴承的振动信号和包络谱，图9为对应的时频谱图，该轴承故障的类型为外圈1 mm宽的横槽. 如图8所示，由于强背景噪声的干扰，无法从时域波形中识别出轴承故障的周期性冲击特征，包络谱图也被噪声频率干扰. 针对图8所示的振动信号，利用本文提出的方法进行故障特征提取，得到基于时频的频带模糊熵. 如图10所示，在f=7 500 Hz处对应的模糊熵值最小，认为是系统的固有主频率. 提取固有频率在时频分析对应的时间响应序列，得到系统的时间波形图和包络谱图。如图11所示，时域波形具有显著的周期性冲击特点，整体背景噪声得到有效抑制. 从包络谱可以看出，故障特征频率76 Hz及倍频非常明显，由此判断出轴承存在故障缺陷.

利用上文提出的自适应带通滤波方法，根据频带模糊熵确定滤波中心频率 ${f_{\rm{c}}} = 7 \,\, 500 \; {\rm{Hz}}$ ，计算滤波器带宽 $\Delta f = 150 \; {\rm{Hz}}$ . 图12为滤波后信号的时域波形图和包络谱图，可以看出，自适应滤波后信号降噪效果显著，能够轻易识别出故障特征频率. 为了对比分析，引入AR预测滤波器，对图8所示的振动信号进行滤波处理. 如图13所示，经过AR滤波器滤波后的时域波形，故障冲击特征表现微弱，不能准确地观察出周期性的冲击故障特征，无法有效抑制强背景噪声.

3 结　语

针对单尺度模糊熵不能敏锐反应信息量变化的不足，基于时频分析的频带多尺度复合模糊熵实现定量地表征非平稳信号的数据信息，反映出不同频带分量在时间轴上的变化特性. 结合自适应带通滤波器，成功实现对滚动轴承微弱信号的特征提取和故障识别，较传统滤波方法在降噪抑制方面具有优势. 仿真模拟和实验信号分析表明，方法能够有效描述振动信号的复杂度和不确定性，做到快速识别滚动承振动信号的冲击特征，具有良好的故障特征识别效果. 由于方法以时频分析为基础，识别精度和诊断效果与所用的时频分析方法紧密相连. 因此在提高方法时频分频率以及模糊熵参数选取方面有待进一步研究.