固有时间尺度分解(Intrinsic Time scale Decomposition, ITD)是由Frei和Osorio[3]提出的另一种自适应时频分析方法(ITD的方法原理本文不再赘述，请参考文献[3])，其能将复杂信号分解为若干相互独立且瞬时频率具有物理意义的固有旋转分量(Proper Rotation Component, PRC)之和。与EMD相比[4-5]，ITD采用线性插值取代了基于三次样条曲线的基线信号构建方法，得到的PRC保留了原始信号的极值时刻，不会出现过包络和欠包络现象，同时端点效应被严格限制在两端。程军圣等[6]在ITD的基本原理上定义了内禀尺度分量(Intrinsic Scale Component, ISC)，提出了基于ITD的改进算法，与能量矩相结合有效识别出齿轮的不同故障类型。考虑到行星齿轮箱信号的复合调制性，Feng等[7]提出了一种联合幅值和频率调制的分析方法，利用ITD将信号分解成若干分量，再对敏感分量的瞬时频率和瞬时幅值进行傅里叶变换得到幅值和频率的解调谱。结合EMD和ITD，Hu等[8]构造出了EITD方法，有效抑制了ITD的端点效应和信号失真。
原始信号经ITD分解之后，进一步突显了故障信号的局部特征信息(即冲击特征)，但噪声抑制作用不明显，高频PRC仍含有大量噪声，直接对高频分量进行故障分析，噪声干扰会严重影响诊断结果的准确性。为了能准确提取故障特征，本文提出一种基于ITD和改进形态滤波(Morphological Filter, MF)的故障诊断方法。首先，利用ITD分解细化原始信号，然后，对有效PRC进行改进形态滤波滤除噪声干扰，准确提取出故障冲击成分，最后，进行分量信号重构避免细节信息的损失。形态滤波是基于数学形态学变换的一种非线性滤波方法[9]，其利用结构元素与信号的形态特征进行匹配，可以有效地提取信号的边缘轮廓和形状特征，运算简单高效。
本文提出的改进形态滤波方法是基于滚动轴承故障信号受谐波和随机噪声干扰严重的特点和形态学基本理论，根据腐蚀、膨胀、开、闭四大基本形态学算子的性质，将对冲击成分具有相同抑制作用但作用方式不同的算子合并加强脉冲抑制效果，得到的2类算子组合作差反向提取出正负冲击，由此构造出一种新的组合差值形态滤波(Combination Difference morphological Filter, CDIF)算子，通过试验分析和比较，选择其中2种CDIF算子的平均作为滤波输出，即为平均组合差值形态滤波(Average Combination Difference morphological Filter, ACDIF)的最终结果。
1 数学形态学 1.1 基本形态学理论 数学形态学是建立在积分几何和随机集论基础上的一种数学分析方法。为了实现滤波和特征提取，对具有滤波窗作用的结构元素和待分析信号进行匹配，当信号局部形态特征和结构元素相匹配时信息得以保留[10-11]。膨胀、腐蚀、开和闭是数学形态学的4种基本算子。经典的形态滤波器有形态梯度(Morphological Gradient, MG)滤波和形态差值滤波(morphological Difference Filter, DIF)。
设f(n)和g(m)分别为定义在F=(0, 1, …, N－1)和G=(0, 1, …, M－1)上的离散函数, 且N≥M。f(n)是一维离散信号, g(m)是结构元素。则f(n)关于g(m)的腐蚀和膨胀运算定义为
 | (1) |
 | (2) |
f(n)关于g(m)的开运算和闭运算分别为
 | (3) |
 | (4) |
式中：n、m表示信号点数；Θ表示腐蚀运算；⊕表示膨胀运算；“°”表示开运算；“·”表示闭运算。
1.2 平均组合差值形态滤波 腐蚀算子具有抑制正脉冲、平滑负脉冲的能力，减小信号峰值，加宽谷域；膨胀算子具有平滑正脉冲、抑制负脉冲的能力，增大信号谷值，扩展峰顶；开算子具有抑制正脉冲、保留负脉冲的能力，滤除信号上尖峰；闭算子具有保留正脉冲、抑制负脉冲的能力，滤除信号下尖峰[12]。考虑到4种基本算子对脉冲的作用及其不同的作用方式，将基本算子按脉冲抑制作用分为2类。将对冲击成分具有相同抑制作用但作用方式不同的算子合并加强抑制效果，再将2类算子组合作差反向提取出正负冲击，由此得到的正负冲击更加准确突出。
首先定义膨胀算子和闭算子的结合：膨胀-闭(Dilation-Closing, DC)或闭-膨胀(Closing-Dilation, CD)滤波算子，定义如下：
 | (5) |
 | (6) |
再定义腐蚀算子和开算子的结合：腐蚀-开(Erosion-Opening, EO)或开-腐蚀(Opening-Erosion, OE)滤波算子，定义如下：
 | (7) |
 | (8) |
为了同时提取正负冲击，取FDC或FCD与FEO或FOE的差值作为新的CDIF算子。CDIF的滤波输出形式主要有以下4种：
 | (9) |
 | (10) |
 | (11) |
 | (12) |
由文献[9, 13]可知，当f(n)满足一定条件时，(f·g⊕g)(n)=(f⊕g)(n)，(f°gΘg)(n)=(f Θg)(n)，即FCD-OE(f(n))=FMG(f(n))。对于同一信号f(n)，每种运算滤波效果各异。通过实验分析比较，发现FCD-EO和FDC-EO的输出信号的冲击幅值较高，相对于其余2种滤波，整体最优。所以，本文决定采用两者的平均滤波输出——ACDIF作为最终结果，即
 | (13) |
FACDIF整合了FCD-EO和FDC-EO的优缺点，既可以有效抑制随机噪声，又能够突出信号冲击特征，兼顾细节信息。
结构元素和形态运算方式是影响形态滤波的2个最主要因素，结构元素的特性与其形状、高度和长度有关。结构元素的形状对滤波效果影响较小[14]，如何选择结构元素的最佳长度成为关键。本文采用结构简单、计算效率高的直线型结构元素，其高度为零。参考文献[14]定义冲击特性幅值：
 | (14) |
式中：Ai为幅值谱中故障频率的i倍频对应的幅值；N为描述A的最高倍频，本文实验中N=5。不同结构元素长度下，滤波后振动信号的A值越大，说明故障特征越明显，结构元素长度的取值越佳。
2 基于ITD-ACDIF的故障诊断方法 ITD能够自适应地将复杂多分量信号分解为若干个PRC，PRC按分解出来的顺序，其包含的频率成分呈现明显的下降趋势。高频部分PRC信号中包含故障信息，但也包含大量噪声。直接将所有高频分量重构进行频率分析，由于噪声污染严重很难提取到故障特征信息，所以需要对PRC进行筛选和滤波去噪。
针对轴承故障振动信号中的特征成分极易被噪声信号淹没而不能及时检测的问题，结合ITD的自适应性，保留信息完整性，突出信号局部故障特征等方面的优势，以及ACDIF在消除随机噪声成分和同时提取正负冲击上的优良特性，提出基于ITD-ACDIF的故障诊断方案。方案流程如图 1所示，具体步骤如下：
 |
| 图 1 基于ITD-ACDIF的故障诊断流程 Fig. 1 Fault diagnosis process based on ITD-ACDIF |
| 图选项 |
步骤1 ?对采样信号x进行ITD分解，得到多个PRC。
步骤2 ?利用峭度指标筛选出包含故障信息最多的分量为有效PRC。
步骤3 ?以冲击特性幅值为依据计算出每个分量信号的最佳结构元素长度Li。
步骤4 ?对有效PRC分别进行ACDIF抑制随机噪声，突出其冲击成分。
步骤5 ?重构形态滤波后的各个分量得到合成信号x′。
步骤6 ?对合成信号x′进行快速傅里叶变换(Fast Fourier Transform, FFT)得到幅值谱，从中提取出故障特征频率，实现故障诊断。
3 仿真分析 为了验证本文提出的ACDIF滤波的可行性和有效性，采用仿真信号x(t)=x1(t)+x2(t)+x3(t)进行分析，其中，x1(t)=cos(40πt)+2cos(110πt)是频率为20 Hz和55 Hz的谐波信号，x2(t)是频率为12 Hz的周期性指数衰减冲击信号，再加上高斯白噪声干扰x3(t)，使得原始信号信噪比为-15 dB，采用高噪声方式以验证本文方法在早期故障下对轻微冲击的提取效果。信号采样频率为2 048 Hz, 采样点数为2 048。x(t)的时域波形和幅值谱如图 2所示。可以看到，冲击信号x2(t)的频率12及24、36、48、60、72 Hz等倍频，谐波信号x1(t)的20 Hz和55 Hz的频率较为突出，尤其55 Hz的频率明显高于冲击信号的频率，高斯噪声的干扰同样严重。
 |
| 图 2 仿真信号时域波形及其幅值谱 Fig. 2 Time-domain waveform and amplitude spectrum of simulated signal |
| 图选项 |
首先对仿真信号x(t)进行经典的形态差值滤波(DIF)[15]，得到滤波后的时域波形及其幅值谱如图 3所示，时域波形中仍含有较多噪声，从幅值谱中也仅看出12 Hz的冲击频率。而原信号经本文所提出的ACDIF(最佳结构元素为22)滤波后，滤波效果显著改善。从图 4中可以看出，ACDIF能够有效抑制原始信号中的谐波成分和随机噪声，同时将正负冲击提取出来，幅值谱中谐波的2个频率基本消除，冲击信号频率及其倍频突出。
 |
| 图 3 仿真信号DIF滤波结果 Fig. 3 DIF filtering results of simulated signal |
| 图选项 |
 |
| 图 4 仿真信号ACDIF滤波结果 Fig. 4 ACDIF filtering results of simulated signal |
| 图选项 |
为进一步验证ITD-ACDIF结合滤波的有效性，对仿真信号x(t)进行ITD分解，得到4个PRC分量和1个残余量R如图 5所示，ITD将噪声和有用信号分离，噪声主要集中在PRC1上。依据峭度准则选取峭度值最大的3个分量(PRC2、PRC3、PRC4)为有效分量，并进行ACDIF形态滤波，以主频率(12、24、36、48、60 Hz)幅值之和筛选出各分量对应的最佳结构元素长度(20、20、30)。在图 6中，有效PRC中的谐波成分和随机噪声都得到了抑制。将滤波后的3个分量相加得到合成信号，从图 7合成信号的时域波形及其幅值谱中可以发现，原始信号中谐波的2个频率得到了有效抑制，随机噪声基本消除。相比于单独使用ACDIF滤波，经过ITD分解细化信号之后再进行有效分量的ACDIF滤波可以保留更多的细节信息，波形较为光滑，冲击信号的频率12 Hz及其二倍频24 Hz更为突出。
 |
| 图 5 仿真信号的ITD分解结果 Fig. 5 ITD decomposition results of simulated signal |
| 图选项 |
 |
| 图 6 有效PRC ACDIF滤波后结果 Fig. 6 ACDIF filtering results of effective PRC |
| 图选项 |
 |
| 图 7 合成信号时域波形及其幅值谱 Fig. 7 Time-domain waveform and amplitude spectrum of synthetic signal |
| 图选项 |
4 应用实例 为了进一步验证本文所提方法在实际旋转机械故障诊断中抑制噪声和提取冲击特征的有效性，采用实际运行状态下轴承从正常到最终失效的全寿命振动数据作实例分析。加速轴承寿命试验机(ABLT-1A)由杭州轴承试验研究中心提供，其由交流电机驱动，在同一根轴上同时进行4个轴承的寿命试验，试验台总体如图 8所示，轴承各参数如表 1和表 2所示, 轴承型号为6308。
 |
| 图 8 试验台 Fig. 8 Test-bed |
| 图选项 |
表 1 轴承结构参数 Table 1 Structure parameters of bearing
| 参数 | 数值 |
| 负载/kN | 12.74 |
| 节径/mm | 65.5 |
| 钢球直径/mm | 15.08 |
| 钢球数 | 8 |
| 接触角/(°) | 0 |
表选项
表 2 轴承外圈故障数据 Table 2 Fault data of bearing outer race
| 参数 | 数值 |
| 采样频率/kHz | 20.0 |
| 采样点数 | 2 048 |
| 发动机转速/(r·min-1) | 4 000 |
| 转轴基频/Hz | 66.67 |
| 故障频率/Hz | 205.29 |
表选项
对外圈故障失效模式下的样本进行分析研究，取故障出现早期的2 048个数据作为原始振动信号进行分析，图 9(a)为外圈信号时域波形，波形较复杂且伴有大量噪声，仅通过时域波形无法了解故障信息。图 9(b)中，由于故障冲击的作用，系统的固有频率已被激起，并且边频成分丰富，调制现象明显，但由于噪声的影响，频率基本集中在中、高频处，低频处的转频和故障频率难以识别，故障特征频率谱线并不明显。
 |
| 图 9 外圈故障信号时域波形及其幅值谱 Fig. 9 Time-domain waveform and amplitude spectrum of outer race fault signal |
| 图选项 |
利用ITD算法对信号进行分解，得到4个频率从高频到低频分布的PRC和1个残余量R，计算各分量峭度值列入表 3中，分量时域波形如图 10所示。利用主频率(205.1、410.2、615.2、820.3、1 016 Hz)幅值之和选出最佳结构元素长度组合(12、8、16)对峭度值最大的前3个PRC(PRC1、PRC2与PRC3)分别进行ACDIF滤波，滤波之后的分量时域波形如图 11所示，各个分量的正负冲击被清晰地提取出来，随机噪声的干扰消除。重构滤波后的PRC得到合成信号，在图 12的幅值谱中，中高频成分得到了抑制，频率峰值集中在低频处，205.1、410.2、615.2、1 016 Hz正好分别对应于外圈故障频率、二倍频、三倍频及五倍频，故障频率及其倍频的谱线清晰且突出，由此可以判断出轴承的外圈故障。
表 3 各PRC及R的峭度值 Table 3 Kurtosis values of PRC and R
| 分量 | 峭度 |
| PRC1 | 3.52 |
| PRC2 | 3.46 |
| PRC3 | 4.83 |
| PRC4 | 2.83 |
| R | 2.44 |
表选项
 |
| 图 10 外圈故障信号的ITD分解结果 Fig. 10 ITD decomposition results of outer race fault signal |
| 图选项 |
 |
| 图 11 外圈故障信号有效PRC ACDIF滤波后结果 Fig. 11 Filtering results of effective PRC of outer race fault signal |
| 图选项 |
 |
| 图 12 外圈合成信号时域波形及其幅值谱 Fig. 12 Time-domain waveform and amplitude spectrum of outer race synthetic signal |
| 图选项 |
利用传统的DIF形态滤波器对各个有效PRC进行形态滤波后合成，得到的最终信号及其幅值谱如图 13所示。与图 12比较可见，利用ACDIF滤波之后的波形更加光滑干净，经DIF处理后的波形仍然残留少量噪声，信号幅值较小，可以得到故障频率但频率幅值远小于ACDIF所得频率幅值，并且倍频难以识别。
 |
| 图 13 ITD-DIF结果信号的时域波形及其幅值谱 Fig. 13 Time-domain waveform and amplitude spectrum of ITD-DIF result signal |
| 图选项 |
以EMD替代ITD作为前置处理方法，先对外圈故障信号进行EMD分解，得到若干内禀模态函数(Intrinsic Mode Function, IMF)分量和残余量R，如图 14所示，噪声和能量都集中在高频IMF分量中，同样以峭度为标准选取前3个IMF分量为主要分量(见表 4)。分别对前3个分量进行ACDIF滤波后重构得到图 15所示的EMD-ACDIF滤波结果，与ITD-ACDIF滤波结果较为相似，冲击成分得以提取，但仅能分辨出故障频率及其二倍频。由此进一步证实了ITD与ACDIF结合的有效性。
 |
| 图 14 外圈故障信号的EMD分解结果 Fig. 14 EMD decomposition results of outer race fault signal |
| 图选项 |
表 4 各IMF分量及R的峭度值 Table 4 Kurtosis values of IMF components and R
| 分量 | 峭度 |
| IMF1 | 4.36 |
| IMF2 | 3.43 |
| IMF3 | 2.98 |
| IMF4 | 2.85 |
| R | 2.57 |
表选项
 |
| 图 15 EMD-ACDIF结果信号的时域波形及其幅值谱 Fig. 15 Time-domain waveform and amplitude spectrum of EMD-ACDIF result signal |
| 图选项 |
最后采用典型的Hilbert包络解调方法[16]对故障信号进行分析。Hilbert解调结果如图 16所示，冲击成分得以提取，但仍存在严重的噪声干扰，提取效果不如本文方法。综上，ITD-ACDIF算法能够有效滤除随机噪声和脉冲干扰，提取出淹没在强背景噪声中的故障特征，提取效果整体上优于其他方法。
 |
| 图 16 Hilbert解调结果信号的时域波形及其幅值谱 Fig. 16 Time-domain waveform and amplitude spectrum of Hilbert demodulation result signal |
| 图选项 |
5 结论 1) 本文将ITD和ACDIF相结合，利用了ITD方法对信号的自适应分解能力和ACDIF的有效滤噪并保留冲击成分的优点，提高了滤波输出的准确度。
2) ACDIF滤波器不仅结构设计简单，而且可应用于振动信号的故障诊断中，通过该方法滤波之后的信号冲击特征被完整保留。
3) 数值仿真和轴承外圈故障振动信号的试验结果表明，本文所提出的ITD-ACDIF故障诊断方法能够较好地将信号的冲击特征从成分复杂的背景噪声中提取出来，实现有效的故障诊断。
参考文献
| [1] | YU J. Local and nonlocal preserving projection for bearing defect classification and performance assessment[J].IEEE Transactions on Industrial Electronics, 2012, 59(5): 2363–2376.DOI:10.1109/TIE.2011.2167893 |
| [2] | 李舜酩, 郭海东, 李殿荣. 振动信号处理方法综述[J].仪器仪表学报, 2013, 34(8): 1907–1915. LI S M, GUO H D, LI D R. Review of vibration signal processing methods[J].Chinese Journal of Scientific Instrument, 2013, 34(8): 1907–1915.(in Chinese) |
| [3] | FREI M G, OSORIO I. Intrinsic time-scale decomposition:Time-frequency-energy analysis and real-time filtering of non-stationary signals[J].Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 2007, 463(2078): 321–342.DOI:10.1098/rspa.2006.1761 |
| [4] | HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903–995.DOI:10.1098/rspa.1998.0193 |
| [5] | 张琳, 黄敏. 基于EMD和切片双谱的轴承故障诊断方法[J].北京航空航天大学学报, 2010, 36(3): 287–290. ZHANG L, HUANG M. Fault diagnosis approach for bearing based on EMD and slice bi-spectrum[J].Journal of Beijing University of Aeronautics and Astronautics, 2010, 36(3): 287–290.(in Chinese) |
| [6] | 程军圣, 李海龙, 杨宇. 改进ITD和能量矩在齿轮故障诊断中的应用[J].振动、测试与诊断, 2013, 33(6): 954–959. CHENG J S, LI H L, YANG Y. Based on the improved ITD and energy moment to diagnose the gear[J].Journal of Vibration, Measurement & Diagnosis, 2013, 33(6): 954–959.(in Chinese) |
| [7] | FENG Z, LIN X, ZUO M. Joint amplitude and frequency demodulation analysis based on intrinsic time-scale decomposition for planetary gearbox fault diagnosis[J].Mechanical Systems and Signal Processing, 2016, 72-73: 223–240.DOI:10.1016/j.ymssp.2015.11.024 |
| [8] | HU A, YAN X, XIANG L. A new wind turbine fault diagnosis method based on ensemble intrinsic time-scale decomposition and WPT-fractal dimension[J].Renewable Energy, 2015, 83: 767–778.DOI:10.1016/j.renene.2015.04.063 |
| [9] | MARAGOS P, SCHAFER R W. Morphological filters-Part Ⅰ:Their set-theoretic analysis and relations to linear shift-invariant filters[J].IEEE Transactions on Acoustics, Speech, and Signal Proceeding, 1987, 35(8): 1153–1169.DOI:10.1109/TASSP.1987.1165259 |
| [10] | 李兵, 张培林, 刘东升, 等. 基于形态梯度解调算子的齿轮故障特征提取[J].振动、测试与诊断, 2010, 30(1): 39–42. LI B, ZHANG P L, LIU D S, et al. Feature extraction for gear fault diagnosis using morphological gradient demodulation[J].Journal of Vibration, Measurement & Diagnosis, 2010, 30(1): 39–42.(in Chinese) |
| [11] | HU Z, WANG C, ZHU J, et al. Bearing fault diagnosis based on an improved morphological filter[J].Measurement, 2016, 80: 163–178.DOI:10.1016/j.measurement.2015.11.028 |
| [12] | 张平, 张小栋, 董晓妮, 等. 自调整复合级联形态滤波算法及应用[J].振动、测试与诊断, 2015, 35(3): 459–463. ZHANG P, ZHANG X D, DONG X N, et al. Self-adjusting composite cascade morphology filter algorithm and its application[J].Journal of Vibration, Measurement & Diagnosis, 2015, 35(3): 459–463.(in Chinese) |
| [13] | HARALICK R M, STERNBERG S R, ZHUANG X. Image analysis using mathematical morphology[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 1987, 9(4): 532–550. |
| [14] | 章立军, 杨德斌, 徐金梧, 等. 基于数学形态滤波的齿轮故障特征提取方法[J].机械工程学报, 2007, 43(2): 71–75. ZHANG L J, YANG D B, XU J W, et al. Approach to extracting gear fault feature based on mathematical morphological filtering[J].Journal of Mechanical Engineering, 2007, 43(2): 71–75.(in Chinese) |
| [15] | RAJ A S, MURALI N. Early classification of bearing faults using morphological operators and fuzzy inference[J].IEEE Transactions on Industrial Electronics, 2013, 60(2): 567–574.DOI:10.1109/TIE.2012.2188259 |
| [16] | FELDMAN M. Hilbert transform in vibration analysis[J].Mechanical Systems and Signal Processing, 2011, 25(3): 735–802.DOI:10.1016/j.ymssp.2010.07.018 |
