Volterra理论在滚动轴承内圈故障程度特征定量提取的研究

王海涛， 王 琨， 史丽晨

(西安建筑科技大学 机电工程学院， 西安 710055)

建筑施工机械结构复杂、工况恶劣，其内部滚动轴承发生故障时，振动信号非平稳非线性程度较为明显。采用合适的分析方法提取振动信号中的有用信息是故障特征提取的关键[1-2]。滚动轴承故障诊断中常用的分析方法有：谱分析、小波分析、EMD(Empirical Mode Decomposition)分解法等[3-4]，基于上述方法的故障特征提取和诊断研究已取得了丰硕的成果。Kankar等[5]基于小波变换和人工智能技术对滚动轴承故障进行了特征提取和诊断，对比验证了支持向量机(Support Vector Machine, SVM)具有更为良好的诊断性能。Yu等[6]提出了基于EMD方法和希尔伯特谱的故障特征提取方法，并与传统包络谱方法进行了对比，得出前者更能对滚动轴承故障特征进行有效提取。然而大量文献主要是对故障部位和类型进行区分，实际上设备的故障是一个从轻微到严重的演变过程，准确识别设备运行中的故障程度演变情况对保证设备平稳运行有着重要意义。近期，对故障程度的研究已经引起了广泛的关注，张龙等[7]基于多尺度熵对滚动轴承故障程度进行了识别，能够较早的发现其早期故障并与故障发展趋势保持了一定的一致性。王玉静等[8]基于EMD和SVD(Singular Value Decompostion)进行特征提取，实现了对滚动轴承不同故障位置及其性能退化程度的特征提取和智能识别。上述研究内容在滚动轴承故障程度定量识别方面取得了一定成效，但依然不能完全满足机械设备实际工况下非线性非平稳等因素导致故障程度特征定量提取相对困难这一问题。针对此问题，本文提出了基于Volterra级数建模和双谱切片分析的滚动轴承故障程度特征提取方法，通过实验对比验证了所述方法的优越性。

1 方法提出

本文利用美国Case Western Reserve University轴承输入输出信号确定Volterra级数模型[9-12]，应用多脉冲激励的方法对Volterra核函数进行求解，并利用广义频率响应函数对核函数进行辨识；建立滚动轴承故障模型，利用双谱分析对Volterra二阶核函数进行表征，针对其双谱三维图和等高线图表征缺少定量结果的缺点，对其双谱三维图进行切片，得到其切片谱，对相位耦合的频率点进行了提取和对比，对故障模型进行了定量分析，使得Volterra核函数表征更为直观，在滚动轴承不同程度故障的对比和诊断中取得了良好效果，基于Volterra级数和双谱分析的特征提取方法流程图如图1所示。

图1 基于Volterra级数和双谱分析的特征提取方法流程图

2 理论简述

2.1 Volterra级数理论

Volterra级数理论最早从数学角度提出，用于研究某些微积分方程；可以利用一个输入值的多维卷积积分的无限长代数和来具体表示系统输出，确定系统模型。对于一个线性系统，系统输入和输出关系在时域中可以表示为以下卷积运算

(1)

式中:μ(t),y(t)分别为系统输入和输出；h(t)为系统的脉冲响应函数。

同理，对任意连续的非线性系统，在初始条件为零的情况下，若系统的输入μ(t)能量有限，相对应的，其系统的响应可以用Volterra级数表示为

(2)

式中：hk(τ1,τ2,…,τk)为k阶Volterra级数的核，是线性系统的脉冲响应函数在高维空间的推广。通过线性系统和非线性系统表示形式的对照，对于线性系统，h(t)代表了线性系统的系统特性；同理，对于非线性系统，h1(τ),h2(τ1,τ2),…hn(τ1,…,τn)所组成的集合便代表了非线性系统的系统特性，确定非线性系统的系统特性便是对非线性系统的Volterra核函数进行确定。

若将式(2)在时域内进行离散化，则可以得到表达式为Volterra级数的离散时域函数

μ(k-m1)…μ(k-mn)

(3)

式中：μ(k),y(k)∈R分别为系统的输入和输出;hn(m1,…，mn)为离散Volterra级数的第n阶核函数。

为了更深入的研究Volterra核函数，仅仅在时域上分析是不够的，故而引入广义频率响应函数的概念，实现Volterra核函数在频域内的研究，是线性频率函数的延伸和拓展。其定义为Volterra核函数hn(τ1,τ2,…,τn)的多维傅里叶变换，如下所示

e-j(ω1τ1+ω2τ2+…+ωnτn)dτ1dτ2…dτn

(4)

利用广义频率响应函数(Generalized Frequency Response Function, GFRF)可以在对系统进行非线性分析辨识时，直观的描述出许多系统的频域特征。

如上文所述，Volterra核函数代表了非线性系统的系统特性，确定系统Volterra核函数便可对系统特性进行表述。对于Volterra核函数的求解，目前主要使用最小二乘法、神经网络等方法，但却暴露出了此类方法计算量大、复杂性高的缺点。针对此类问题，本文提出了改进的多脉冲激励法对Volterra级数核进行求解，更为方便的确定系统核函数。

2.2 模型辨识

2.2.1 建模与求解

在对系统进行建模时，可将m,n取有限值。当不计算直流分量y0的影响时，式(3)可以改写为具有截断形式的离散时域形式

(5)

式中：N为非线性系统的最高阶次；M为脉冲响应的记忆长度；e(k)为非线性系统被截断成有限项后产生的截断误差。

假设存在一非线性系统，而且可以用二阶Volterra级数进行表述。如式(6)所示

(6)

将一个单脉冲激励作为系统输入，即μ0(t)=aδδ(t)，则通过以下公式对系统一阶核进行求解

μ0(t)=aδδ(t)→y0(t)

y0(t)=aδh1(t)

(7)

式中：aδ为脉冲增益参数。

对式(7)移项后可得

(8)

现假设此非线性系统Volterra级数核h1,h2均存在，将两个不同幅值的脉冲激励作为输入进行求解

(9)

(10)

将式(9)，式(10)利用矩阵进行表示为

(11)

则h1(t),h2(t,t)可表示为

(12)

2.2.2 算例验证

通过上述理论研究，假设一非线性系统如下

(13)

vt=0.5μt+2.5μt-1+1.7μt-2

(14)

假设系统记忆长度为M=3，那么对yt在μt点进行泰勒级数展开，并将式(14)代入后得

yt=0.5+1.25(0.5μt+2.5μt-1+1.7μt-2)-

(15)

从式(15)中可以得到系统的一阶和二阶核函数原始数值如下

(16)

对系统进行两次输入，分别为σ(t)、σ(t)+σ(t-T)；系统输出分别记为y0和y1。利用式(12)可以计算得到系统的Volterra核函数估计值为

(17)

为了清晰地反映出辨识结果对上述系统的逼近程度，利用式(4)的离散形式求出核函数的广义频率响应，并做出系统一阶和二阶核函数理论值和辨识值的频率响应图。

通过对图2和图3中理论值和估计值的频率响应图的对比后发现，一阶核最大绝对误差幅值和相位的频率响应约为0.4dB和1°，二阶核最大绝对误差幅值和相位的频率响应约为0.04dB和0.3°，说明通过多脉冲激励的方法对非线性系统有良好的辨识效果，可以利用此方法对非线性系统进行建模。

当然，上述理论情况是建立在忽略外界噪音的基础上的。在实际应用中，噪声的影响难以忽略，直接利用系统核函数的频率响应图进行对比，难以消除噪声以及各频率点的频率耦合问题。虽然二阶核GFRF频率响应图可以表示系统特性，但由于二阶核频响图是三维谱图，对其中有用信息往往没有得到较好的利用[13]。本文利用双谱有效抑制高斯噪声的特性，能够对轴承故障特征进行有效提取，并对其进行了切片处理，将其有用信息投影到了二维图形上，对系统特征的定量提取和直观表达起到了重要作用。

图2 一阶核理论值和辨识值对比

Fig.2 Comparison of first-order theoretical and identification values

(a) 理论二阶核

(b) 辨识二阶核

(c) 误差对比

2.3 双谱及其切片谱

高阶累积量谱常称为高阶谱，假设k阶累积量ckx(τ1,τ2,…,τk-1)绝对可和，则k阶累积量谱定义为其k-1维傅里叶变换，即为

τk-1)e-j2π(f1τ1+f2τ2+…+fk-1τk-1)

(18)

由于在实际应用中累积量及其谱的值需要估计，但随着阶数的增加，计算量也随之大大增加，所以一般只考虑k≤4的情况。当k=3时，对应累积量的谱为三谱，即所谓的双谱。

三阶谱(双谱)表达式为

(19)

将式(19)转换到频域，利用频率特性函数进行表达

B(ω1,ω2)=γ3aH(ω1)H(ω2)H*(ω1+ω2)

(20)

式中:H(ω)为频率特性函数；γ3a为歪斜度；a(n)为一独立同分布的非高斯随机过程。

双谱Bx(f1,f2)具有如下重要性质

Bx(f1,-f1-f2)=Bx(-f1-f2,f1)=

Bx(f2,-f1-f2)

(21)

式中：“*”为复共轭。由式(21)可知，双谱Bx(f1,f2)的对称线包括f1=f2、f1=0、f2=0、2f1=-f2、2f2=-f1、f1=-f2六条，从而将双谱的定义域分为12个区域，其中区域Q={(f1,f2)|0≤f1,f2≤f1}称为双谱Bx(f1,f2)在平面(f1,f2)的主域，从而对称性是双谱的重要性质，只要知道主域内的双谱，便可通过对称关系对其他双谱值进行求解。

显然，双谱反映了信号的三阶矩(信号偏度)在双频率平面内的分布情况，且双谱是一个复数，所以它既包含了幅值信息，又包含了相位信息，可以用于描述系统的二次相位耦合程度，这是功率谱所不具备的。以此来对不同的Volterra故障模型进行比对，实现故障诊断的目的。

3 实验论证

本文将Volterra级数与双谱切片谱引入滚动轴承故障诊断，对正常与不同程度故障轴承的双谱及其切片谱图进行对比，验证此方法的可行性和优越性。

本实验的数据信号来自于Case Western Reserve University轴承故障中心，试验中驱动端输入轴轴承型号为SKF6205,轴承规格见表1。对轴承内圈采用电火花加工单点点蚀故障。实验中使用加速度传感器采集振动信号，通过使用磁性底座将传感器安放在电机壳体上。加速度传感器分别安装在电机壳体的驱动端和风扇端12点钟的位置。振动信号是通过16位通道的DAT记录器采集的，并且后期在MATLAB环境中处理。数字信号的采样频率为12 kHz。本文选取正常轴承与不同损伤程度的内圈故障轴承进行实验，故障直径分别为0.18 mm和0.36 mm，故障深度为0.28 mm，驱动端负载取0，驱动端转速为1 797 r/min。对系统给予输入信号分别为σ(t)、σ(t)+σ(t-1)，运行系统并利用加速度传感器采集振动输出信号，对实验观测数据进行抽样，截取2 000组数据进行实验。利用多脉冲激励方法在MATLAB中对输入输出数据进行分析，分别得到正常轴承和内圈故障轴承系统的Volterra核函数h1,h2，利用二阶核h2便可以反映出正常轴承和故障轴承系统的特性。

表1 SKF6205深沟球轴承规格

做出正常轴承和内圈故障轴承二阶核的双谱三维图和等高线图，如图4、图5所示，可以观察到正常轴承和损伤程度不同的故障轴承的双谱三维图明显不同，谱峰的形状及能量随着损伤程度的增加而逐渐增多。损伤程度越大的轴承二次耦合频率点更多、更分散，这一点在其等高线上可以更明显的看出。利用双谱三维图和等高线图揭示出了其频率二次耦合的位置，表征出了正常和不同损伤程度的故障轴承之间系统特性的区别。

(a) 正常轴承

(b) 损伤直径0.18 mm

(c) 损伤直径0.36 mm

(a) 正常轴承

(b) 损伤直径0.18 mm

(c) 损伤直径0.36 mm

由于双谱三维图和等高线图只可以表征出正常轴承和故障轴承不同的系统特性，即仅可以定性的对正常轴承和故障轴承进行比较，而缺乏更为定量的系统特征描述；并且从双谱三维图中可以看出，其显示的特征较为复杂。为了更加明了地观察出正常轴承系统和不同程度故障轴承系统的区别，对双谱进行对角切片分析。选取f=f1=f2的位置进行切片，如图6所示，既提取出了双谱三维图的有用信息，又有效地减少了运算量，并且包含了双谱的全部信息，可以对正常和不同损伤程度的故障轴承有一个更为清晰的对比，从而达到故障诊断的目的。

(a) 正常轴承

(b) 损伤直径0.18 mm

(c) 损伤直径0.36 mm

分析双谱切片图可知，正常轴承其双谱切片图只有一个峰值，出现在f=107 Hz这一频率点处；而损伤直径为0.18 mm的轴承双谱切片图峰值分布分散，在f=52 Hz、f=77 Hz、f=101 Hz等处都出现了明显峰值；对于损伤直径为0.36 mm的轴承双谱切片图分布更为分散，峰值位置在f=41 Hz、f=54 Hz、f=73 Hz、f=99 Hz等处明显。通过以上分析说明二次相位耦合位置和频率接近的点随着故障损伤程度变化，耦合现象越来越严重，与正常轴承有着明显的差异。

在实际应用中，考虑到设备自身的影响，在上述实验的基础上，对采集到的信号加入不同信噪比SNR=10 dB和SNR=5 dB的外部噪声，以此来模拟平稳运行无明显冲击振动和非平稳运行设备自身对实验结果的影响。利用本文方法对数据进行建模分析，做出Volterra核函数的双谱切片图，如图7所示，并与上述实验结果进行对比。

(a) 正常轴承

(b) 损伤直径0.18 mm

(c) 损伤直径0.36 mm

对比分析图6和图7后发现，加入信噪比SNR=10 dB的外部噪声之后，应用本文所述方法对Volterra核函数进行双谱切片分析，相比较图6，图7中虽然出现了一定的波动，但整体上未对其峰值所对应的频率点产生影响，正常轴承所对应的双谱切片图中依旧在f=107 Hz出现了最高峰值，同样的，损伤直径为0.18 mm的轴承双谱切片图在f=53 Hz、f=77 Hz、f=101 Hz等附近也都出现了明显峰值；而对于损伤直径为0.36 mm的轴承双谱切片图，其峰值位置分布依旧出现在了f=41 Hz、f=54 Hz、f=70 Hz、f=102 Hz等附近，可以描述不同的二次相位耦合现象，对比不同损伤程度的故障轴承。

当加入信噪比为SNR=5 dB的外部噪声后，做出不同损伤程度的双谱切片图，其峰值排列杂乱无序，有用信息被噪声淹没，难以对故障特征进行提取，但是其更多的细节信息应隐含在三阶及以上的高阶Volterra核函数中，本文对此不作详细描述。

总体来讲，未加外部噪声和加入一定信噪比的外部噪声的分析结果并无明显出入，虽然在外部条件的影响下出现了一部分的杂乱信号，但对其二次频率耦合的频率点位置并无太大影响，说明在设备平稳运行无明显冲击产生的条件下，本文方法可以忽略一部分设备自身所产生的影响，对提取不同损伤程度的故障轴承特征有一定的适用性。

现采用小波和包络谱的方法对相同实验条件下的轴承进行分析，即对加入信噪比SNR=10 dB的外部噪声后的输出信号进行小波降噪和包络谱分析，并与本文所述方法进行对比。首先利用式(22)以及前文所述轴承数据求得内圈故障特征频率为finner=159.7 Hz。选取小波基为“db6”并并进行Hilbert包络谱分析，如图8所示，并与前文所述方法进行对比。

(22)

(a) 正常轴承

(b) 损伤直径0.18 mm

(c) 损伤直径0.36 mm

对比图8中正常轴承和不同损伤程度的轴承包络谱图，发现图8(b)和图8(c)中出现最高峰值的频率点为161.7 Hz，与理论上计算的故障频率点基本吻合，可以依据此频率点来区分正常轴承和内圈故障轴承。随着故障程度的加深，虽然图8(b)和图8(c)中除最高峰(161.7 Hz)外其余波动也随故障程度的加深而增大，但其分布却无规律可寻，难以提取隐含在其中的更多故障信息。与此方法相比，本文所述方法最突出的优点在于利用Volterra级数建模和以二阶核函数双谱切片图中所表示的二次相位耦合频率作为故障特征，从而对不同损伤程度的故障轴承进行对比和区分，不需要进行滤波处理。另外，得益于双谱分析对噪声的不敏感，能够在平稳运行无明显冲击影响的条件下对轴承进行特征提取，在一定程度上克服了设备自身的影响。

4 结 论

本文提出了一种基于Volterra级数理论和双谱切片分析提取故障特征的方法，利用Case Western University轴承输入输出数据进行了验证，通过改进的多脉冲激励法对Volterra核函数进行求解，而后采用双谱切片对二阶核函数进行表征，并与传统包络谱方法进行了对比。研究结果表明：由于Volterra核函数的确定是利用了系统输入和输出信号，可以更为完善的表征系统状态，不完全依赖滤波处理；并且利用双谱切片分析，可以直观的表示出核函数的特征，提取出故障轴承有用信息，实现正常轴承和不同损伤程度故障轴承的对比和辨识；相比较传统方法，本文所述方法能够在设备平稳运行无明显冲击振动的情况下消除一部分设备自身的影响，具有更好的抗噪性，能够实现滚动轴承内圈故障程度特征定量提取。

参 考 文 献

[1] 李娟，周东华，司小胜，等. 微小故障诊断方法综述[J]. 控制理论与应用, 2012, 29(12): 1517-1529.

LI Juan, ZHOU Donghua, SI Xiaosheng, et al.Review of incipient fault diagnosis methods[J].Control Theory & Applications, 2012, 29(12): 1517-1529.

[2] 王国彪，何正嘉，陈雪峰，等. 机械故障诊断基础研究“何去何从”[J]. 机械工程学报, 2013, 37(1): 63-72.

WANG Guobiao, HE Zhengjia, CHEN Xuefeng, et al. Basic research on machinery fault diagnosis-what is the prescription[J].Chinese Journal of Mechanical Engineering, 2013, 37(1): 63-72.

[3] 蔡艳平，李艾华，石林锁，等. 基于EMD与谱峭度的滚动轴承故障检测改进包络谱分析[J]. 振动与冲击, 2011, 30(2): 167-172.

CAI Yanping,LI Aihua,SHI Linsuo, et al. Roller bearing fault detection using poved anaiysis based on }n and spectxum envelope spectxum kurtosis[J]. Journal of Vibration& Shock, 2011, 30(2): 167-172.

[4] MENG Lingjie, XIANG Jiawei, ZHONG Yongteng, et al.Fault diagnosis of rolling bearing based on second generation wavelet denoising and morphological filter[J]. Journal of Mechanical Science & Technology, 2015, 29(8): 3121-3129.

[5] KANKAR P K， SHARMA S C, HARSHA S P. Rolling element bearing fault diagnosis using wavelet transform[J]. Neurocomputing, 2011, 74(10): 1638-1645.

[6] YU Dejie， CHENG Junsheng， YANG Yu. Application of EMD method and Hilbert spectrum to the fault diagnosis of roller bearings[J]. Mechanical Systems & Signal Processing, 2005, 19(2): 259-270.

[7] 张龙，黄文艺，熊国良. 基于多尺度熵的滚动轴承故障程度评估[J]. 振动与冲击，2014, 33(9): 185-189.

ZHANG Long,HUANG Wenyi,XIONG Guoliang.Assessment of rolling element bearing fault severity using multi-scale entropy[J]. Journal of Vibration& Shock， 2014, 33(9): 185-189.

[8] 王玉静，姜义成，康守强，等. 基于优化集合EMD的滚动轴承故障位置及性能退化程度诊断方法[J]. 仪器仪表学报，2013, 34(8): 1834-1840.

WANG Yujing,JIANG Yicheng,KANG Shouqiang et al. Diagnosis method of fault location and performance degradation degree of rolling bearing based on optimal ensemble EMD[J]. Chinese Journal of Scientific Instrument， 2013, 34(8): 1834-1840.

[9] 彭志科，程长明. Volterra级数理论研究进展与展望[J]. 科学通报, 2015(20): 1874-1888.

PENG Zhike, CHENG Changming.Volterra series theory: a state-of-the-art review[J]. Chinese Science Bulletin, 2015(20): 1874-1888.

[10] WANG Qisheng， WANG Keyan， CHEN Shaojun. Least squares approximation method for the solution of Volterra-Fredholm integral equations[J]. Journal of Computational & Applied Mathematics, 2014, 272(3): 141-147.

[11] 魏瑞轩，韩崇昭，张优云，等. 非线性系统故障诊断的Volterra模型方法[J]. 系统工程与电子技术， 2004, 26(11): 1736-1738.

WEI Ruixuan, HAN Chongzhao, ZHANG Youyun,et al. Volterra model method for fault diagnosis of nonlinear system[J]. Systems Engineering and Electronics, 2004, 26(11): 1736-1738.

[12] 韩海涛，马红光，曹建福，等. 多输入多输出非线性系统Volterra频域核的非参数辨识方法[J]. 西安交通大学学报, 2012, 46(10): 66-71.

HAN Haitao, MA Hongguang, CAO Jianfu, et al.A non-parametric identification method of volterra frequency domain kernels for MIMO nonlinear system[J]. Journal of Xi’an Jiaotong University, 2012, 46(10): 66-71.

[13] 窦唯，刘占生. 旋转机械故障诊断的图形识别方法研究[J]. 振动与冲击, 2012, 31(17): 171-175.

DOU Wei, LIU Zhansheng. A fault diagnosis method based on graphic recognition for rotating machinery[J]. Journal of Vibration& Shock,2012, 31(17): 171-175.