滚动轴承振动的非线性超混沌特性研究

李兆飞，任小洪，黄臣程

(1.四川理工学院 自动化与电子信息学院，四川 自贡 643000；2.人工智能四川省重点实验室，四川 自贡 643000)

研究滚动轴承振动的非线性演化动态特性有助于认识其内在变化规律，并利用这些规律对其进行诊断和监测，具有一定的普遍意义。随着对轴承振动非线性机理的认识逐渐加深，采用关联维数[1]、最大Lyapunov指数[2-3]、Kolmogorov熵及相关复杂性测度[4-6]等方法对振动的非线性混沌特性进行判断，很好地解释了轴承振动现象的复杂运动现象，判定轴承振动是一种在有限相空间轨道永不重复、形态复杂的运动。然而，这些特征量只能表现轴承振动的某方面特性，而且受信号长度及噪声影响，这些特征量的计算结果会有所偏差，导致判断错误。因此，尝试采用递归图、CLY方法和功率谱3种方法对轴承不同振动状态实测信号的混沌特性进行试验分析。

如前所述，目前对轴承振动的混沌特性研究多是基于混沌的辨别，但对其振动超混沌特性的分析还未见报道，分析轴承振动的超混沌特性有助于研究振动在相空间多个不同方向的复杂演化规律。通常通过计算实测滚动轴承振动信号的Lyapunov指数谱判断其超混沌特性。Lyapunov指数谱的计算往往要通过重构相空间，之后再判定未知系统方程状态。判断未知系统方程通常采用人工神经网络[7-9]，存在训练时间较长、易陷入局部极小点及过拟合，且确定网络结构困难的问题[10]。因此，采用训练耗时较小且较少出现过拟合现象的最小二乘支持向量机[11](Least Squares Support Vector Maehine，LS-SVM)判别系统方程计算信号Lyapunov指数谱的方法，分析轴承振动的超混沌特性。

1 轴承振动信号采集及相空间重构

1.1 轴承振动信号的采集

采用Case West Reserve University滚动轴承实验室实测的内、外圈和滚动体的不同故障程度与转速振动数据库，该数据库采用深沟球轴承6205-2RS JEM SKF，电动机空载，采样频率fs=12 kHz，转速为1 797 r/min，轴承故障由电火花加工机在球轴承内、外圈及钢球上模拟损伤性故障，故障直径分别为0.18，0.36，0.54和0.72 mm(文中选取内圈模拟点蚀故障并在后续分析中分别对应简化表示为f1～f4)，故障深度为0.28 mm。正常状态和内圈故障振动信号的时域波形如图1所示，并分别取1 024个点研究不同故障类型及不同程度故障轴承振动的混沌特性。

(b)内圈故障

1.2 相空间重构技术

相空间重构技术是混沌特性研究的前提。混沌吸引子是混沌系统的特征之一，体现了混沌系统的运动规律，相空间重构目的就是为了在高维空间中恢复混沌吸引子。系统运动是有n个变量的动力系统，即一组n个变量的一阶微分方程

(1)

由连续变量坐标x(t)及其(n-1)阶导数构成的n维相空间就表示了系统随时间演变的状态空间，即

X(t)=[x(t),x1(t),…,xn-1(t)]T。

(2)

Ruelle用离散序列x(t)及其(n-1)时滞位移将其代替，构成一个新的n维嵌入相空间，即

X(t)= [x(t),x(t+τ),…,x(t+(m-1)τ)]T

,(3)

式中：τ为时延；m为重构相空间的维数，m≥2d+1(d为状态空间的关联维数，也称拓扑维)。对xn(n=1,2,…,N)，取定τ，其重构相空间为

Xl= [xl,xl+τ,…,xl+(m-1)τ]T=

(4)

式中：l=N-(m-1)τ。求取τ使原序列时延后能作为独立坐标，可根据平均位移[14]、自相关[13]、互信息[12]及(去偏)复自相关技术等计算。好的重构相空间是使重构后的吸引子与系统真正的吸引子拓扑等价。m是嵌入维，可根据预测误差最小[16]、几何不变量[14]、虚假邻近点[17]及其改进后的Cao方法[18]计算。确定时延与嵌入维有2种争议：一种认为可先求得时延再选择嵌入维；另一种则认为时延与嵌入维是相关的。研究者提出较常用的C-C[19]及其改进法[20-21]能同时计算时延与嵌入窗。因为不能得到混沌系统先验知识，这些方法均有一定的主观性[22]。

分析中采用自相关法[13]求τ，根据xn与xn+τ有些独立又不完全不相关的特性，为使其能在重构相空间中作为独立坐标处理，计算自相关函数为

(5)

则C(τ)首次下降到初值的1-1/e时的延时为最优τ，不同嵌入维下正常及内圈故障振动信号的自相关函数分别如图2所示，从图中可以看出，正常状态最佳延时为3，不同内圈故障状态的最佳延时均为2。

(b)内圈故障

再用虚假邻近点法[17]计算m，先选定m，求得Xl的最近邻点Xη(l)，即

‖Xl-Xη(l)‖m=min{‖Xl-Xj‖:j=1,2,…,l,l≠j}，

(6)

将m增加到m+1，当满足

(7)

Xη(l)即为Xl的虚假邻近点。(7)式说明当相空间从m维演变到m+1维时，2个相点相差很大距离是因为高维混沌吸引子中2个不相邻的点投影到低维轨道上时变成相邻的2个临近虚假点造成的。求取最佳嵌入维数m时一般取RT=10[17]，统计每个m时其虚假邻近点占所有重构向量数的比例P(m)，直到P(m)不再随m的增加而减小或小于某个值(如l%或5%)时，吸引子被认为完全展开。

在不同嵌入维下正常及故障状态统计的P(m)如图3所示。从图中可以看出，当m=5和m=4时，正常状态及内圈故障状态的P(m)不再明显减少，因此正常状态最佳嵌入维取5，故障状态最佳嵌入维均取4。

(a)正常状态

根据嵌入维和延时，采用相空间重构技术绘制正常和内圈故障振动信号的混沌吸引子如图4所示。从图中可以看出，混沌吸引子在有限的空间不断的缠绕，呈现一定自相似的状态。

图4 滚动轴承振动信号混沌吸引子

2 滚动轴承振动的混沌特性分析

首先采用递归图对实测轴承不同振动状态的相关性和确定性进行检验，然后使用CLY方法判定轴承振动的混沌特性，最后结合功率谱对振动进行频谱分析。

2.1 递归图

为分析序列周期及非周期成分的相对大小，揭示动力系统内部相似性的先验知识，文献[23]提出了递归图(Recurrence Plot，RP)方法。周期系统信号被噪声污染时可能表现出非周期性，这表明系统不会准确地回到以前的状态变量，但可能返回很接近以前某时刻的状态，因此可用来表现周期及非周期的程度。由(4)式采用嵌入技术把信号重构为向量Xl(l=1,2,…,N-(m-1)τ)。然后求Xl的任意2个时刻i,j(i,j=1,2,…,N-(m-1)τ)的距离，即

dij=‖Xi-Xj‖ 。

(8)

如果dij≤r，以(i,j)为纵坐标做一点，即绘制出递归图，表现重构的轨线重复或递归其自身的信息，表示了系统的时间关联情况。试验中，计算正常状态及内圈故障状态信号时分别取r=2和r=1，得到的递归图如图5所示。

图5 滚动轴承振动信号的递归图

与文献[24]中周期序列、混沌序列及随机序列的递归图进行对比分析发现：轴承正常状态的递归图存在较明显的平行于对角线的线段，但带状线不是连续的，说明正常振动为非严格周期状态，存在随机性因素(部分为噪声干扰)；内圈故障状态递归图也存在较明显的平行于对角线的线段，与混沌系统类似，表明其也具有周期成分，但介于周期与随机之间，振动为混沌特征，且故障状态越大，递归图与白噪声信号越相似；但轴承正常和内圈故障状态周期性均不明显，由于绘制递归图受噪声和r的影响，其周期和非周期性大小只可作为定性的判断，还需要参考其他方法进行综合的分析判定。

2.2 CLY方法

为判别信号的混沌特性，以虚假邻域方法为基础，文献[25]提出了CLY算法：由(4)式将xn(n=1,2,…,N)嵌入到m维的相空间，对选定的m，由 (6) 式计算每个Xl(m)(l=1,2,…,N-(m-1)τ)的最近邻点Xη(l)(m)，得到平均预测误差

。(9)

定义函数E(m)=E*(m+1)/E*(m)。对选择的任意m，通常随机信号的E(m)=1或在其上下波动；而对混沌系统，随m的增大稳定在1附近。故根据m-E(m)曲线可简单的判别信号的混沌特性。正常及内圈故障状态振动信号的m-E(m)曲线如图6所示，从图中可以看出，当m≥5时E(m)趋于1，可简单判别出不同振动状态都存在混沌特性。

2.3 功率谱

对xn(n=1,2,…,N)加上周期条件xn=xn+j，求其自相关函数(离散卷积)

(10)

然后对Cn进行离散Fourier变换，计算Fourier变换系数[26]

(11)

图6 滚动轴承振动信号的m-E(m)曲线

图7 不同信号的功率谱

正常和内圈故障状态振动信号的功率谱如图8所示，从图中可以看出：正常状态振动信号有较明显的独立尖峰谱线，表明出较强的周期性，而存在连续的谱线，又表明其也具有混沌性，部分谱线分布类似于随机信号谱线，则对应于其中由噪声引起的随机成分；内圈故障状态振动信号的功率谱类似于Lorenz序列的功率谱图，其功率谱线不同于周期和随机信号的离散谱，表现为宽带连续谱，还表现出在高频段随频率指数衰减的特性，因此故障状态存在很强的混沌性，但周期性并不明显。综合分析可知，轴承不同状态振动信号的频率分布特性存在显著的差异，正常状态具有明显的周期性，而内圈故障状态具有明显的混沌特性，与递归图判定结论基本一致。

图8 滚动轴承振动信号的功率谱

3 滚动轴承振动的超混沌特性分析

Lyapunov指数谱λi(i=1,2,…,m-1)是表示相空间各方向轨道旁的膨胀及收缩的平均度量，由系统长时间演变决定，反映了初始时刻2个无限靠近的点随时间分离的情况，每个λi都是相空间不同方向上相对运动局部变形的平均，不同特征的运动与λi的关系见表1。

混沌吸引子的局部不稳定和某方向上的指数发散或膨胀相对应，有2个或2个以上正的λi系统就为超混沌，否则系统为非混沌。现采用LS-SVM判别系统方程方法得到轴承振动的Lyapunov指数谱[27-28]，研究滚动轴承不同状态振动的超混沌特性。

表1 不同运动与Lyapunov指数的关系

3.1 LS-SVM算法

由(4)式重构xn(n=1,2,…,N)在m维欧式空间的一条轨道Xl(l=1,2,…,N-(m-1)τ)，则Xl能再现原系统的动力学特性，设Xl到Xl+τ的映射关系为G[6]，即

(12)

显然Xi+τ的前m-1个数值即为Xi的后m-1个数值，只要确定Xl→Xl+mτ的关系F即可确定G，设

xl+mτ=F(Xl) 。

(13)

1)LS-SVM辨识未知系统

采用LS-SVM辨识F。先由原始信号xn(n=1,2,…,N)构造训练数据集(Xl,xj+τ)，其中l=1,2,…,N-mτ。利用高维特征空间的线性函数拟合样本集[29]

xl+mτ=ωTφ(Xl)+b，

(14)

则原空间中的约束条件和优化问题可表达为

s.t.xj+mτ=ωTφ(Xl)+b+el，

(15)

式中：el为误差，其组成的误差向量为e，偏差为b，权重为C。用Lagrange乘子al有

b+ej-xj+mτ}。

(16)

由KKT条件并消除el和ω，得到线性方程组

(17)

式中：Y=[x1+mτ,x2+mτ,…,xN]T；α=[α1,α2,…,αN-mτ]T；I=[1,1,…,1]T；E为N-mτ阶单位阵；l=(ki,j)(N-mτ)×(N-mτ)，采用径向基函数得到K的每个元素

(18)

解方程(17)得到α和b，则LS-SVM回归函数逼近 (12) 式中的F为

(19)

2)计算Lyapunov指数谱

由 (12) 式得到

ΔXl+τ=DG(Xl)·ΔXl，

(20)

(21)

式中：DG(Xn)为映射G的Jacobian矩阵；ΔXl表示重构的向量Xl在l时刻的微小变化。通过DG(Xl)矩阵，在l+τ时刻重构相空间的向量Xl+τ值的微小变化中将被表现出来。继续下去，该作用会累积到l+rτ(r∈N)时刻重构相空间的向量Xl+rτ取值的变化，其关系为

ΔXl+rτ=DG(Xl+(r-1)τ)…DG(Xl)·ΔXl。(22)

(23)

按大小排列λi(i=1,2,…,m)。因为 (20) 式中矩阵会产生分数幂和指数，所以，ΔXl+τ往往是病态阵，不能求得其全部精确特征值。可利用长乘积矩阵的分解方法对 (23) 式求解。首先定义

Ar=DG(Xl+(r-1)τ)…DG(Xl)=A(r)A(r-1)…A(1)

，(24)

对此长乘积矩阵，递归计算

A(i)Q(i-1)=Q(i)R(i),i=1,2…r，

(25)

式中：Q(i)为正交矩阵；R(i)为上三角矩阵；Q(0)为m阶单位阵。按 (25) 式方法将QR分解r次，则Ar的QR分解为Ar=Q(r)R(r)R(r-1) …R(1)，故 (22) 式的本征值可以通过下式求出

(26)

进而求得系统的Lyapunov指数谱为

(27)

3.2 计算结果分析

采用前文给出的实测滚动轴承正常和内圈故障状态的振动信号，分析其振动的超混沌特性。计算同样采用相关法和虚假邻近点法分别确定延时τ和嵌入维m。LS-SVM计算参数取为：C=1 000，a=1，计算出的Lyapunov指数谱见表2。

表2 不同状态振动的Lyapunov指数谱

由表可知，轴承正常及内圈故障状态振动均为超混沌，但正常状态振动的超混沌特征不明显。当故障增大时，轴承超混沌特性程度也增强，并且表现的更加明显。

此外，对轴承外圈、滚动体及保持架不同故障状态进行了相同的混沌及超混沌试验分析，均得到同样的结果。因此轴承振动的混沌程度与故障大小存在一定的对应关系，正常状态振动比较平稳有序，较大故障的振动状态会增加混乱程度。

4 结束语

应用4种方法，从多个方面对实测轴承正常及故障状态振动的非线性混沌特性进行了判定和分析，研究表明：轴承正常及故障状态振动的混沌特征具有明显的差异性，研究轴承振动的非线性混沌特性要结合不同的振动状态。上述研究结果有助于理解轴承振动的复杂演化规律，为今后使用其振动信号混沌特性提取故障特征及进行故障诊断奠定一定的基础。