2. 中航工业中国航空综合技术研究所, 北京 100028
2. AVIC China Aero-polytechnology Establishment, Beijing 100028, China
当滚动轴承发生早期局部故障时，由于滚道上的缺陷，会激发一系列周期性或准周期性脉冲，周期性的脉冲特征被认为是滚动轴承故障的重要标志[1]，从振动信号中提取出周期脉冲特征是诊断过程中的关键步骤。然而，由于机械部件在运行过程中受到其他噪声的干扰以及传递路径的影响，故障特征非常微弱。为解决这一问题，最小熵解卷积(Minimum Entropy Deconvolution, MED)被引入以提升微弱的冲击特性。
MED是一种盲解卷积方法，它最早由Wiggins[2]提出以增强地震反射数据。它旨在最大化信号中的峭度值以提取微弱脉冲，同时最小化其他噪声分量的峭度。Sawalhi等[3]首次证明了MED检测滚动轴承故障的有效性，并且在提升信号脉冲分量方面有着良好的效果。随后，MED得到了广泛的应用[4-9]。但是在MED中，存在着一些问题。滤波器的目标是使得信号的峭度最大化，但是当滤波长度不合适时，它倾向于最大化信号中的随机单一脉冲成分。而MED的滤波长度选取目前没有明确的方法。文献[10]表明，滤波长度越长，信号的峭度值越大。但是，在其研究结果中，只是根据试验数据得到了一个滤波长度的选取范围，没有准确的滤波长度选取方案。Cheng等[11]在对MED进行优化时，使用了一个经验公式来确定滤波长度的大小，但是，在其原始文献[12]中，并没有对滤波长度的影响进行深入的探讨研究。
目前许多使用MED方法进行故障诊断的研究中，缺乏对MED滤波长度的深入讨论，例如Abboud等[13]在使用MED方法对轴承大量样本进行诊断时，滤波长度是根据经验值选取的。Cheng等[14]在比较MED与其他盲解卷积算法时，滤波长度的值是给定的，没有针对滤波长度进行对比分析。Jiang等[15]提出了通过l0范数优化MED的方法，然而该方法也将滤波长度设为定值，没有研究不同滤波长度对算法的影响。一些研究者在使用MED算法时，忽略了滤波长度对信号的影响[16-17]。
为了使MED尽可能地提升周期脉冲成分，避免提升单一随机噪声脉冲，提出一种MED滤波长度自动选择的新方法，结合自相关函数提出了一个能量判定标准，来自适应地实现MED中最优滤波长度的选取问题。不同滤波长度下的能量判定标准是不同的，该标准可以使得周期性脉冲信号和噪声信号最大化地分离。在最优滤波长度下，能够提升连续周期性脉冲成分。通过仿真信号详细地比较了不同滤波长度下MED输出信号中所包含的周期脉冲信息。探讨了不同数据长度与滤波长度之间的关系，并给出了一个最优滤波长度自适应选取方案。最后，将该方法应用于滚动轴承故障诊断，通过滚动轴承全寿命疲劳加速试验和远离轴承故障源的滚动轴承故障试验，充分验证了该方法在诊断轴承早期微弱故障中的有效性。
1 MED方法简介 1.1 算法回顾当轴承表面发生局部缺陷时，滚动体通过缺陷部位会激发一系列周期或准周期脉冲，这些脉冲信号与系统噪声混杂在一起，由于结构传输路径和传感器位置的影响，脉冲信号会非常微弱，因此，在检测弱信号方面，使用MED来近似地恢复这些脉冲是有效的手段之一。MED系统模型如图 1所示。
|
| 图 1 MED的系统模型 Fig. 1 System mode of MED |
图 1中，s=[s1 s2…sN]代表故障周期脉冲, N代表数据长度；n=[n1 n2…nN]为噪声干扰；受到系统谐波的影响，h=[h1 h2…hN]；传感器采集到的观测信号x=[x1 x2…xN]可表示为
| $ \mathit{\boldsymbol{x}} = (\mathit{\boldsymbol{s}} + \mathit{\boldsymbol{n}}) * \mathit{\boldsymbol{h}} $ | （1） |
式中：“*”代表卷积。MED假定系统输入s的熵值较小，经过系统后得到的x序列熵值增加。因此，解卷积就是需要寻找一个大小为L的有限脉冲响应滤波器f=[f1 f2…fL]T，经过滤波器后的输出y=[y1 y2…yN]能够近似于原系统的输入s，即
| $ \begin{array}{*{20}{l}} {y(j) = \sum\limits_{l = 1}^L f (l)x(j - l) \approx s(j)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j = 1,2, \cdots ,N} \end{array} $ | （2） |
MED的实现，主要有特征向量法和目标函数法。应用较多的是目标函数法，可以通过优化峭度目标函数得到最优效果，其形式为
| $ {O_4}(\mathit{\boldsymbol{f}}) = \frac{{\sum\limits_{j = 1}^N {{y^4}} (j)}}{{{{[\sum\limits_{j = 1}^N {{y^2}} (j)]}^2}}} $ | （3） |
令其一阶导数为0可以得到最佳滤波器：
| $ \frac{{{\rm{d}}{O_4}(\mathit{\boldsymbol{f}})}}{{{\rm{d}}\mathit{\boldsymbol{f}}}} = 0 $ | （4） |
式(2)的矩阵形式为
| $ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{X}}_0^{\rm{T}}\mathit{\boldsymbol{f}} $ | （5） |
式中：
| $ {\mathit{\boldsymbol{X}}_0} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_1}}&{{\mathit{\boldsymbol{x}}_2}}&{{\mathit{\boldsymbol{x}}_3}}& \cdots &{{\mathit{\boldsymbol{x}}_N}}\\ 0&{{\mathit{\boldsymbol{x}}_1}}&{{\mathit{\boldsymbol{x}}_2}}& \cdots &{{\mathit{\boldsymbol{x}}_{N - 1}}}\\ 0&0&{{\mathit{\boldsymbol{x}}_1}}& \cdots &{{\mathit{\boldsymbol{x}}_{N - 2}}}\\ \vdots & \vdots & \vdots &{}& \vdots \\ 0&0&0& \cdots &{{\mathit{\boldsymbol{x}}_{N - L + 1}}} \end{array}} \right]_{L \times N}} $ |
将式(3)和式(4)代入式(5)可得：
| $ \begin{array}{l} \mathit{\boldsymbol{f}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\sum\limits_{j = 1}^N {y_j^2} }}{{\sum\limits_{j = 1}^N {y_j^4} }}{({\mathit{\boldsymbol{X}}_0}\mathit{\boldsymbol{X}}_0^{\rm{T}})^{ - 1}}{\mathit{\boldsymbol{X}}_0}{[\mathit{\boldsymbol{y}}_1^3\quad \mathit{\boldsymbol{y}}_2^3\quad \cdots \quad \mathit{\boldsymbol{y}}_N^3]^{\rm{T}}} \end{array} $ | （6） |
MED算法具体实现流程为
步骤1 初始化滤波器f0，代入信号x得到X0。
步骤2 设定迭代次数上限mmax、迭代终止阈值S和滤波长度L。
步骤3 根据式(5)，代入X0和滤波器参数fm计算出信号ym，m是当前迭代次数。
步骤4 根据式(6)迭代得到fm+1。
步骤5 根据式(3)计算O4(fm)和O4(fm+1)，计算迭代误差ΔE=O4(fm+1)-O4(fm)。
步骤6 若m < mmax且ΔE < S，则进入步骤3继续循环迭代，否则输出最终滤波器参数fend。
步骤7 代入滤波器参数和x，根据式(5)计算得到y作为s的近似。
1.2 滤波长度对MED的影响为了直观说明滤波长度L在MED算法中的重要性，建立仿真振动信号：
| $ x(t) = s(t) + n(t) + h(t) $ | （7） |
如图 2所示，图 2(a)为故障周期脉冲成分s(t)，故障采样间隔为50个点；n(t)为高斯白噪声，s(t)与n(t)信噪比为0.416，其混合之后的波形如图 2(b)所示；图 2(c)为观测信号x(t)，其中谐波成分为
|
| 图 2 模拟故障信号 Fig. 2 Simulated fault signal |
| $ \begin{array}{*{20}{l}} {h(t) = 0.1{\rm{sin}}(2\pi {f_1}t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.2{\rm{sin}}(2\pi {f_2}t) + {\rm{sin}}(2\pi {f_3}t)} \end{array} $ | （8） |
信号长度为2 000采样点，f1=2f2=4f3=0.067。
使用MED对仿真信号进行滤波，滤波长度分别取L=50和L=150，其输出结果如图 3所示。从结果可以看出，当L=50时，输出信号可以看到连续脉冲成分，但是信号中伴随着较大的噪声干扰。然而在L=150的情况下，输出信号是单一的随机脉冲，对于旋转机械的故障诊断，人们期望从微弱信号中恢复出一系列周期性或准周期性故障脉冲特征。因此，这种信号是无意义的。此类现象说明，滤波长度对于MED的输出结果有很大的影响。选择合适的滤波长度是MED方法中关键的一步，为避免单一脉冲最大化的发生，应当研究明确的选择方法，凭借经验值给定滤波长度值很可能会干扰故障诊断的结果。
|
| 图 3 滤波长度为50和150时MED的输出结果 Fig. 3 MED outputs when the filter lengths are 50 and 150 |
由前面分析结果可知，MED在不同滤波长度下，输出结果可能截然不同，因此，为了充分发挥MED提升微弱周期脉冲的潜力，需要建立衡量MED中输出信号周期性的方法。本文结合自相关函数方法来解决滤波长度选取问题。
2.1 自相关函数分析自相关分析反映了信号本身在不同时刻或阶段的相似性, 是判定信号是否具有周期性的有效方法。设a(t)为被测周期振动信号，b(t)为宽带高斯白噪声振动信号，可观测到的振动信号为
| $ c(t) = a(t) + b(t) $ | （9） |
若c(t)为正弦函数，a(t)叠加了不相关的噪声b(t)，即
| $ c(t) = a(t) + b(t) = A{\rm{sin}}(\omega t + \varphi ) + b(t) $ | （10） |
则其自相关可表示为
| $ \begin{array}{*{20}{l}} {{R_c}(\tau ) = {R_a}(\tau ) + {R_b}(\tau ) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathop {{\rm{lim}}}\limits_{T \to \infty } \frac{1}{{2T}}\int_{ - T}^T a (t)a(t - \tau ){\rm{d}}t + {R_b}(\tau ) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{A^2}}}{2}{\rm{cos}}(\omega t) + {R_b}(\tau )} \end{array} $ | （11） |
式中：T为a(t)信号的周期；A为信号幅度；ω为角频率大小；τ为时间t的延迟；φ为信号初始角度大小。当b(t)为宽带噪声时，Rb(τ)集中表现在τ=0附近，结果如图 4(a)和图 4(b)所示。如果a(t)为周期函数，则Ra(τ)仍为周期性函数。在图 5(a)和图 5(b)中，当τ较大时，Rc(τ)只反映Ra(τ)的情况。这样就可以由τ较大时的Rc(τ)测量出a(t)的幅度和频率。
|
| 图 4 宽带高斯白噪声信号及其自相关函数 Fig. 4 Gaussian white noise signal with wideband and its autocorrelation function |
|
| 图 5 含噪声的周期信号及其自相关函数 Fig. 5 Periodic signal with noise and its autocorrelation function |
为了使MED尽可能地提升微弱信号中的连续周期性脉冲成分，结合自相关函数，提出了一个能量判定指标，避免出现使用MED后提升单一随机脉冲的现象。假设在给定先验滤波长度l的情况下，MED的输出信号序列为yl，那么原始信号中剩余的其他噪声成分则可以表示为
| $ {\mathit{\boldsymbol{D}}_l} = \mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{y}}_l} $ | （12） |
式中：x代表原始振动信号。最优滤波判定标准可以通过式(13)选取：
| $ {L_\mu } = \frac{{\sum\limits_{j = 1}^N {R_{{\mathit{\boldsymbol{y}}_l}}^2} (j)}}{{\sum\limits_{j = 1}^N {R_{{\mathit{\boldsymbol{D}}_l}}^2} (j)}} $ | （13） |
式中：Ryl(j)代表输出信号的自相关函数值，RDl(j)代表剩余信号的自相关函数值。当MED出信号接近周期性信号时，由于滤波器使得信号的峭度值最大，微弱的周期脉冲成分得到提升，而系统谐波和噪声成分在此刻得到抑制，时域波形上的周期脉冲幅值增大，其输出信号的自相关函数也接近周期性波形。当MED输出信号为单一脉冲时，除了在时域波形某时刻为较大幅值的冲击特征外，剩余部分为噪声成分，所以在其自相关函数的表现中，幅值主要反映在0点处，其余时刻的信号幅值接近0。因此可以通过衡量输出信号自相关函数的能量来区分滤波效果。信号能量比值的大小可以作为输出结果是否为周期性脉冲的判断依据。为了排除干扰，将自相关函数在0时刻的幅值去除。
在计算Lμ之前，需要给出一个滤波长度L的范围。针对式(7)仿真信号，计算滤波长度L从2到500之间的Lμ值以及每个滤波长度下输出信号的整体峭度值，结果如图 6所示。
|
| 图 6 不同滤波长度下仿真信号的峭度和Lμ值 Fig. 6 Kurtosis andLμof simulated signal at different filter lengths |
从图 6(a)可以看出，随着滤波长度的增加，信号的峭度值也增加，这和文献[10]给出的结果一致，这个现象很容易理解，因为MED的目标函数就是将信号的峭度最大化。但是，仅凭借信号的峭度信息去选取滤波长度是不可行的。因为峭度的变化与输出信号是否为连续脉冲没有关联。图 6(b)反映了不同滤波长度下Lμ值的变化，其中A、B、C和D四点的滤波长度分别为50、104、166和167。其滤波结果如图 7所示。从图 6和图 7可以看出，当Lμ值较大时，例如B点和C点，输出信号的结果中周期故障脉冲更加清晰，而在Lμ较小的A点，滤波结果中噪声成分较多。但是在D点，输出结果为单一脉冲，此时相比于A点，Lμ值更小。值得注意的是，C点和D点的滤波长度值相差仅仅为1，而信号的输出结果却截然不同，这再次说明，如果L的选取不合适，将会产生错误的结果，因此使用MED时确定合适的滤波长度是非常必要的。在D点之后，Lμ趋于下降的稳定趋势。这意味着当滤波长度大于等于167时，MED的输出结果将会是单一脉冲，它将会失去提升周期性微弱信号的能力。
|
| 图 7 滤波长度在A、B、C和D点时的输出结果 Fig. 7 Output results when filter lengths are A, B, C and D |
在使用Lμ确定MED最优滤波长度之前，需要预先给定L一个范围，如果预先选择的L范围过大，会增加计算成本，如果L范围过小，则可能找不到最佳的滤波长度。为了确定Lμ在不同数据长度下的适用性以及研究L的先验范围，接下来通过增加仿真信号的长度来说明。
当仿真信号长度为2 000采样点时，其结果已经在图 6(b)中展示，在D点之后，MED将会失去其提升周期脉冲的能力，因此将D点定义为滤波长度失效点(注意失效点不代表在此之前所有滤波长度都是合适的，在失效点前有可能也会有些滤波长度失效，例如图 6(b)中在C点之前的一小段长度)。若失效点滤波长度值表示为Lf，则MED的滤波先验有效范围与信号长度的比值为
| $ K = \frac{{{L_{\rm{f}}}}}{N} \times 100\% $ | （14） |
图 8展示了仿真信号长度从2 000采样点等间隔变化到20 000采样点的情况下，滤波长度先验有效范围与数据长度比值的变化趋势。从结果来看，滤波长度的先验有效范围与信号长度的比值在6%~10%之间，其中最大的范围比为9.7%，最小的则为6.1%。并且随着信号长度的增加，比值有略微下降的趋势，有效先验滤波长度范围倾向于变小，一方面，这表明MED在提升微弱周期脉冲成分时，滤波长度的选取范围不能无限大，仅仅是在很小的一个范围内。另一方面也表明Lμ的结果是较为稳定的，波动范围是可接受的。因此，为了尽可能达到最好的效果，使用Lμ确定MED最优滤波长度时，建议使用先验有效滤波长度范围不超过原始数据长度的10%。
|
| 图 8 有效滤波长度与数据长度比值 Fig. 8 Ratio of effective filter length to signal length |
通过上述分析，本文提出的MED最优滤波长度自适应选择方法在原先MED方法流程上主要的变化为
1) 在滤波前，根据原始信号长度N，确定先验有效滤波长度范围不大于0.1N。
2) 在先验滤波长度范围中，计算Lμ的值，判断Lμ的是否达到最大值，选取Lμ最大值对应的滤波长度作为最优结果。
3) 在最优滤波长度下得到输出信号，从而使MED自适应地提升微弱信号中的周期故障脉冲成分。具体流程如图 9所示。
|
| 图 9 MED最优滤波长度选取流程图 Fig. 9 Flow chart of selecting MED optimal filter length |
实际工况下，滚动轴承信号较为复杂，故障信息受到噪声和传递路径影响较为严重，为了充分说明所提方法的有效性，通过两组情形下的滚动轴承故障试验进行验证：
1) 基于辛辛那提试验中心滚动轴承全寿命疲劳加速试验。目的是验证所提方法在早期微弱故障检测方面的适用性。
2) 基于航空发动机转子试验器机匣信号的滚动轴承故障试验。目的是验证所提方法在远离故障源的情况下，检测微弱故障的能力。
3.1 滚动轴承全寿命早期微弱故障检测全寿命疲劳加速试验数据由NSFI/UCR智能维护系统中心(IMS)[18]提供。试验所用轴承型号为Rexnord ZA-2115，轴承参数在表 1中给出。试验台由4个安装在轴上的滚动轴承组成，通过摩擦带连接到电机转动，径向受载为267 kN，转速恒定为2 000 r/min。试验台结构如图 10所示。PCB 353B33高灵敏度ICP加速度传感器放置在每个轴承座上，采样频率为20 480 Hz，每个样本采样点数为20 480，采样间隔时长为10 min。
| 型号 | 节径/mm | 滚子直径/mm | 滚子数 | 接触角/rad |
| ZA-2115 | 71.5 | 8.4 | 16 | 0.264 8 |
1号轴承在超过设计寿命后发生外圈故障，其中数据记录为984个样本。1号轴承的有效值RMS(Root Mean Square)变化趋势如图 11所示，计算前200个正常运行样本有效值的均值μ和方差σ，将失效阈值设为μ+3σ，从结果来看，在533样本时，有效值超出阈值范围。说明轴承状态可能发生了变化。因此，选取533时刻的振动信号作为检测的样本。该样本时域波形如图 12所示。在文献[13]的研究中，同样将533时刻样本作为早期故障的标志，然而遗憾的是在应用其所提的方法MED+SK (Spectral Kurtosis)+SES (Square Envelope Spectrum)对533样本进行诊断时，图 13结果中并没有检测出轴承的外圈故障特征频率(Ball Pass Frequency on Outer Race, BPFO)成分。文献[13]中作者将MED的滤波长度根据经验设置为1 000。
|
| 图 11 1号轴承有效值RMS的变化 Fig. 11 Change in RMS of No.1 bearing |
|
| 图 12 1号轴承533样本时域波形 Fig. 12 Time domain waveform of record 533, No. 1 bearing |
使用本文所提的方法对533样本进行诊断，首先根据数据长度20 480，预先设置先验有效滤波长度为2 000。图 14中分别展示了先验滤波范围内，信号峭度值与Lμ的变化情况。从图 14(a)可以看出，峭度值是与滤波长度正相关的。然而图 14(b)中，滤波长度失效值为210，因此，对于此样本，其真实的滤波有效范围占该样本长度的1%。在失效点之后，Lμ呈稳定的低值状态，可以看出，最优滤波长度取值应为61，而在滤波长度等于1 000时刻，显然此时已经处于MED失效阶段。在图 15中分别绘制了滤波长度为61和1 000情况下MED的输出结果。从时域波形上看，在1 000滤波长度下，时域结果变成单一随机脉冲，并且脉冲幅值非常大，此时MED将信号中的随机单一脉冲的峭度最大化，而有用的故障信息则被视为其他的噪声成分受到了抑制。相比于滤波长度为61情况下，尽管时域波形中还存在着噪声，但是有用信息还保留在信号中。从图 16两者的包络谱中可以直观地看出区别，在最优滤波长度61时，包络谱中BPFO (236 Hz)非常突出，并且其倍频成分也都能从噪声中甄别出来。而在滤波长度为1 000的情况下，BPFO淹没在噪声频率中。因为此时MED的输出信号已经失去了提升周期故障脉冲的能力。尽管在输出信号的峭度上表现为更大的值，但是在追求高峭度值的同时会导致输出结果变成单一随机大脉冲。而对于SK方法，它对于峭度值非常敏感，随机的单一大脉冲会干扰其诊断结果，因此这也是文献[13]中为什么使用MED+SK+SES方法没有检测到故障特征的原因。在最优滤波长度下，MED+SK+SES方法同样能诊断出微弱的故障特征，其结果在图 17中展示。
|
| 图 14 不同滤波长度下533样本的峭度值和Lμ Fig. 14 Kurtosis and Lμof record 533 at different filter lengths |
|
| 图 15 不同滤波长度下533样本的输出结果 Fig. 15 Output results of record 533 under different filter lengths |
|
| 图 16 不同滤波长度下533样本的包络谱 Fig. 16 Envelope spectra of record 533 under different filter lengths |
|
| 图 17 滤波长度L=61时MED+SK+SES的结果 Fig. 17 Results of MED+SK+SES when filter length L=61 |
通过上述分析可知，在最优滤波长度下，MED可以提升轴承微弱的周期性故障脉冲特征，为了进一步证明本文所提方法的有效性，选取1号轴承正常运行时的样本进行MED最优滤波，以验证该方法在提升微弱故障冲击的准确性。
图 18为正常运行情况下第200样本的Lμ值变化趋势，可以看到，在滤波长度为5时，Lμ值最大，此后的阶段，基本上维持在较小的平稳趋势。这是因为正常运行样本中没有周期冲击成分，因此在相关之后的滤波信号与噪声信号的能量比值上差异不大。该样本经过MED最优滤波之后的包络谱如图 19所示，频谱中没有故障特征频率的存在。
|
| 图 18 不同滤波长度下正常运行样本的Lμ Fig. 18 Lμ of normal record at different filter lengths |
|
| 图 19 正常运行样本MED最优滤波之后的包络谱 Fig. 19 Envelope spectrum of normal record after being filtered by MED at optimal filter length |
在一些情况下，振动传感器往往不能安装到靠近轴承的位置。例如在航空发动机中，振动传感器常放置于机匣外壳。这会给滚动轴承诊断带来新的挑战，由于传递路径及其他结构噪声的影响，轴承故障信息传递到机匣后会变得更加微弱。为了模拟远离轴承振动源的情况，通过带机匣的航空发动机转子试验器来进行滚动轴承故障试验，从而验证最优滤波长度下MED的效果。
航空发动机转子试验器如图 20所示，其结构接近某型真实发动机，比例为1:3。该试验器结构简化，转子支撑结构为0-2-0结构；两个带叶片的盘分别代表压气机盘和涡轮盘，并且分别由两个轴承进行支撑；压气机端为滚子轴承，涡轮端为球轴承；其支撑结构为弹性结构，并且支撑刚度可以调节；没有燃烧室结构，是单转子系统模型，由电机带动支撑。
|
| 图 20 航空发动机转子试验器及剖面图 Fig. 20 Aero-engine rotor tester and its sectional view |
故障轴承放置在涡轮端，为6206型球轴承，并且通过人工线切割在外圈和内圈上植入损伤, 图 21为植入故障的球轴承。其主要参数如表 2所示。分别在轴承座和机匣上放置B&K 4805 ICP振动传感器，如图 22所示，机匣上传感器位置有水平测点方向和垂直测点方向。数据采集器为NI-USB 9234，其中数据采样频率为10.24 kHz，数据长度为8 192，转速为1 500 r/min。
| 型号 | 节径/mm | 滚珠直径/mm | 滚珠数 | 接触角/rad |
| 6206 | 46 | 9.5 | 9 | 0 |
|
| 图 21 外圈故障和内圈故障轴承 Fig. 21 Rolling bearings with outer race and inner race faults |
|
| 图 22 机匣测点及轴承座测点 Fig. 22 Casing measuring points and bearing housing measuring point |
在1 500 r/min垂直测点方向上，图 23(a)为轴承座时域波形，图 23(b)为机匣信号时域波形。从时域上可以看到，轴承座故障信号冲击特征明显，振幅大。而在机匣信号中，冲击特征被大量噪声掩盖，幅值非常微弱。单从时域信号中，无法看出故障信息。在两者包络谱中，轴承座信号中的BPFO (92 Hz)以及其倍频非常明显，而在机匣信号中的包络谱中，转子的传输路径使得故障特征衰减，再加上旋转过程中结构振动的噪声干扰，虽然可以看到BPFO，但由于大量的无关频率成分的干扰，无法清楚地确定轴承是否故障。因此检测远离轴承振动源信号中的故障特征较为困难。
|
| 图 23 1 500 r/min时外圈故障轴承时域波形及包络谱 Fig. 23 Time domain waveforms and envelope spectra of outer race fault bearing at 1 500 r/min |
对图 23(b)中机匣信号进行最优MED滤波，确定其最优滤波长度。图 24反映了不同滤波长度下Lμ的变化趋势。其中在L=155时Lμ达到最大值。并且从200之后，Lμ呈稳定的较小值。这表明，对于1 500 r/min垂直测点的机匣信号，有效的滤波长度在155~200之间。MED最优滤波输出结果如图 25(a)所示，从时域信号中，可以看出，连续周期性冲击特性从噪声中突显出来，并且在其包络谱图 25(b)上，BPFO及其倍频成分非常清晰，与图 23(c)相比，在3BPFO和4BPFO故障特征倍频表现上更为突出。然而，在图 25(c)中，L=154的情况下，输出结果却截然相反，MED出现了随机单一脉冲，并且在其包络频谱图 25(d)中，无法确定其故障特征。
|
| 图 24 1 500 r/min时外圈故障机匣垂直测点信号不同滤波长度下Lμ的变化 Fig. 24 Variation of Lμ at different filter lengths of outer race fault casing vertical measuring point signal at 1 500 r/min |
|
| 图 25 1 500 r/min时外圈故障机匣垂直测点信号不同滤波长度下的结果 Fig. 25 Output results of outer race fault casing vertical measuring point signal at different filter lengths |
对于内圈故障，在1 500 r/min转速下机匣水平测点信号上，图 26展示了轴承座故障信号与机匣故障信号的对比情况。可以发现，图 26(a)轴承座时域信号中故障冲击明显，振幅较大。在其对应的包络谱图 26(c)中，可以看到清晰的内圈故障特征频率135 Hz(Ball Pass Frequency on Inner Race, BPFI)及其倍频成分，并且伴随着显著的转频调制现象(fr=25 Hz)。然而在图 26(b)机匣信号的时域波形上，噪声成分多，振幅小，故障冲击极其微弱，转速谐波为主要成分。在图 26(d)对应的包络谱中，无法识别BPFI，多为噪声频率。图 27为Lμ的变化趋势。最优滤波长度在L=100处，其输出结果如图 28(a)所示，可以看到机匣信号中的连续冲击特性较为明显地从噪声中提升出来，在其包络频谱图 28(b)中，BPFI、2BPFI和3BPFI都得以凸显，并且转速调制现象也非常明显。然而在L=101时，Lμ的值较低，MED的效果将大打折扣，无法提升出明显的周期冲击特征，在图 28(c)时域信号中表现为单一的大脉冲。在图 28(d)包络谱中BPFI及其倍频非常微弱。
|
| 图 26 1 500 r/min时内圈故障轴承时域波形及包络谱 Fig. 26 Time domain waveforms and envelope spectra of inner race fault bearing at 1 500 r/min |
|
| 图 27 1 500 r/min时内圈故障机匣水平测点信号不同滤波长度下Lμ的变化 Fig. 27 Variation of Lμ at different filter lengths of inner race fault casing horizontal measuring point signal at 1 500 r/min |
|
| 图 28 1 500 r/min时内圈故障机匣水平测点信号不同滤波长度下的结果 Fig. 28 Output results of inner race fault casing horizontal measuring point signal at different filter lengths at 1 500 r/min |
由上述结果可知，滤波长度的选取对于MED的输出结果是极其重要的。合适的滤波结果可以充分提升连续的周期性冲击特征，而错误的滤波结果会使得MED丧失其提升微弱信号中冲击的能力，这会对诊断轴承微弱故障带来不利的影响。本文所提方法可以避免此现象的发生。
为了进一步说明本文所提方法的优越性，选取最大相关峭度解卷积(Maximum Correlated Kurtosis Deconvolution, MCKD)方法对两组试验的故障信号进行诊断，从而与所提方法进行比较。MCKD方法由McDonald等[19]提出, 其主要思想是在MED的基础上，将算法中的峭度目标函数变换为相关峭度目标函数，避免MED容易提升单一脉冲的现象，增强指定周期的冲击信号。现已广泛应用于旋转机械的故障诊断中[20-21]。限于篇幅，MCKD具体算法请参照文献[19]，本文不再赘述。
MCKD中，需要提前设置故障提升周期T、时移周期参数M以及滤波长度L。其中T对算法的影响最大，T决定了信号的提升冲击周期。因此为了使得MCKD得到最佳的效果，故障提升周期T取相应的轴承故障特征周期。对于533样本，T=86.7。对于机匣信号，外圈故障T=111.3；内圈故障T=75.8。参数M取5，滤波长度与所提方法的选取相同。
图 29展示了MCKD对两组不同故障试验的诊断结果。其中图 29(a)和图 29(b)分别为对533样本滤波后的时域波形和包络谱。从结果可以看出，MCKD可以提取出相应的BPFO，但其在提取故障特征的倍频上，与图 16(a)相比表现较弱。图 29(c)为MCKD滤波后外圈故障机匣信号的时域波形，与本文所提方法的结果图 25(a)相比，其时域波形上噪声成分较多，故障冲击特征不明显。其对应的包络谱图 29(d)中，尽管可以看到BPFO及其倍频，但是噪声频率成分比较多，显然图 25(b)的诊断效果更好。同样的对于内圈故障机匣信号，图 29(e)为MCKD滤波后的时域波形，其包络谱如图 29(f)所示。从结果来看，图 29(f)中可以看到相应的BPFI，但是与图 28(b)的结果相比，显然本文所提方法中的BPFI及其倍频更加突出和明显，相应的调制频率也更加清晰。
|
| 图 29 使用MCKD诊断两组故障试验信号的结果 Fig. 29 Diagnosis results of two sets of failure testing signals using MCKD |
1) 滚动轴承全寿命试验结果表明，凭借经验值选取的滤波长度很容易在诊断过程中产生错误的结果，而本文所提出的自适应方法可以避免此类问题并且成功地在轴承发生故障的早期阶段识别并提升微弱故障特征。
2) 远离故障源的滚动轴承外圈和内圈故障试验结果表明，所提方法可以消除传递路径的影响，成功地提升航空发动机转子试验器机匣信号中微弱的周期性故障冲击特征。
3) 与MCKD方法比较的结果表明，本文所提方法在诊断两组试验数据的效果上更有优势。
本文所提的方法有利于提升MED在微弱冲击信号处理方面的正确性和高效性。
| [1] | RANDALL R B, ANTONI J. Rolling element bearing diagnostics-A tutorial[J]. Mechanical Systems and Signal Processing, 2011, 25(2): 485-520. |
| Click to display the text | |
| [2] | WIGGINS R A. Minimum entropy deconvolution[J]. Geoexploration, 1978, 16(1-2): 21-35. |
| Click to display the text | |
| [3] | SAWALHI N, RANDALL R B, ENDO H. The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[J]. Mechanical Systems and Signal Processing, 2007, 21(6): 2616-2633. |
| Click to display the text | |
| [4] |
王宏超, 陈进, 董广明. 基于最小熵解卷积与稀疏分解的滚动轴承微弱故障特征提取[J]. 机械工程学报, 2013, 49(1): 88-94. WANG H C, CHEN J, DONG G M. Fault diagnosis method for rolling bearing's weak fault based on minimum entropy deconvolution and sparse decomposition[J]. Journal of Mechanical Engineering, 2013, 49(1): 88-94. (in Chinese) |
| Cited By in Cnki (95) | Click to display the text | |
| [5] |
陈海周, 王家序, 汤宝平, 等. 基于最小熵解卷积和Teager能量算子直升机滚动轴承复合故障诊断研究[J]. 振动与冲击, 2017, 36(9): 45-50, 73. CHEN H Z, WANG J X, TANG B P, et al. Helicopter rolling bearing hybrid faults diagnosis using minimum entropy deconvolution and Teager energy operator[J]. Journal of Vibration and Shock, 2017, 36(9): 45-50, 73. (in Chinese) |
| Cited By in Cnki (26) | Click to display the text | |
| [6] |
林桐, 陈果, 滕春禹, 等. 基于机匣振动信号的滚动轴承故障协同诊断技术[J]. 航空动力学报, 2018, 33(10): 2376-2384. LIN T, CHEN G, TENG C Y, et al. Rolling bearing collaborative fault diagnosis technology for casing vibration signal[J]. Journal of Aerospace Power, 2018, 33(10): 2376-2384. (in Chinese) |
| Cited By in Cnki (3) | Click to display the text | |
| [7] |
姚成玉, 来博文, 陈东宁, 等. 基于最小熵解卷积-变分模态分解和优化支持向量机的滚动轴承故障诊断方法[J]. 中国机械工程, 2017, 28(24): 3001-3012. YAO C Y, LAI B W, CHEN D N, et al. Fault diagnosis method based on MED-VMD and optimized SVM for rolling bearings[J]. China Mechanical Engineering, 2017, 28(24): 3001-3012. (in Chinese) |
| Cited By in Cnki (16) | Click to display the text | |
| [8] |
崔伟成, 张征. 基于局部特征尺度分解与最小熵解卷积的轴承故障诊断[J]. 轴承, 2018(5): 51-55. CUI W C, ZHANG Z. Fault diagnosis for bearings based on LCD-MED[J]. Bearing, 2018(5): 51-55. (in Chinese) |
| Cited By in Cnki (4) | Click to display the text | |
| [9] |
郭家宇, 熊炘, 刘浩. 基于VMD和MED的滚动轴承故障诊断方法[J]. 轴承, 2018(6): 50-54. GUO J Y, XIONG X, LIU H. Fault diagnosis method for rolling bearings based on VMD and MED[J]. Bearing, 2018(6): 50-54. (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [10] | BSRSZCZ T, SAWALHI N. Fault detection enhancement in rolling element bearings using the minimum entropy deconvolution[J]. Archives of Acoustics, 2012, 37(2): 131-141. |
| Click to display the text | |
| [11] | CHENG Y, ZHOU N, ZHANG W H, et al. Application of an improved minimum entropy deconvolution method for railway rolling element bearing fault diagnosis[J]. Journal of Sound and Vibration, 2018, 425: 53-69. |
| Click to display the text | |
| [12] | MIAO Y H, ZHAO M, LIN J, et al. Application of an improved maximum correlated kurtosis deconvolution method for fault diagnosis of rolling element bearings[J]. Mechanical Systems and Signal Processing, 2017, 92: 173-195. |
| Click to display the text | |
| [13] | ABBOUD D, ELBADAOUI M, SMITH W A, et al. Advanced bearing diagnostics:A comparative study of two powerful approaches[J]. Mechanical Systems and Signal Processing, 2019, 114: 604-627. |
| Click to display the text | |
| [14] | CHENG Y, WANG Z W, ZHANG W H, et al. Particle swarm optimization algorithm to solve the deconvolution problem for rolling element bearing fault diagnosis[J]. ISA Transactions, 2019, 90: 244-267. |
| Click to display the text | |
| [15] | JIANG X X, XU C, SHI J J, et al. A new l0-norm embedded MED method for roller element bearing fault diagnosis at early stage of damage[J]. Measurement, 2018, 127: 414-424. |
| Click to display the text | |
| [16] | ZHANG L, HU N Q. Fault diagnosis of sun gear based on continuous vibration separation and minimum entropy deconvolution[J]. Measurement, 2019, 141: 332-344. |
| Click to display the text | |
| [17] | WANG S H, XIANG J W, TANG H S, et al. Minimum entropy deconvolution based on simulation-determined band pass filter to detect faults in axial piston pump bearings[J]. ISA Transactions, 2019, 88: 186-198. |
| Click to display the text | |
| [18] | QIU H, LEE J, LIN J, et al. Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics[J]. Journal of Sound and Vibration, 2006, 289(4-5): 1066-1090. |
| Click to display the text | |
| [19] | MCDONALD G L, ZHAO Q, ZUO M J. Maximum correlated Kurtosis deconvolution and application on gear tooth chip fault detection[J]. Mechanical Systems and Signal Processing, 2012, 33: 237-255. |
| Click to display the text | |
| [20] |
宿磊, 黄海润, 李可, 等. 基于LCD-MCKD的滚动轴承故障特征提取方法[J]. 华中科技大学学报(自然科学版), 2019, 47(9): 19-24. SU L, HUANG H R, LI K, et al. Feature extraction of fault rolling bearings based on LCD-MCKD[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2019, 47(9): 19-24. (in Chinese) |
| Cited By in Cnki | Click to display the text | |
| [21] |
李政, 张炜, 明安波, 等. 基于IEWT和MCKD的滚动轴承故障诊断方法[J]. 机械工程学报, 2019, 55(23): 136-146. LI Z, ZHANG W, MING A B, et al. A novel fault diagnosis method based on Improved Empirical Wavelet Trans-form and Maximum Correlated Kurtosis Deconvolution for rolling element bearing[J]. Journal of Mechanical Engineering, 2019, 55(23): 136-146. (in Chinese) |
| Cited By in Cnki (2) | Click to display the text |