CN102156873A - 一种基于混沌的机械零部件早期单点故障检测与分类方法 - Google Patents

Abstract

本发明为一种基于混沌的机械零部件早期单点故障检测与分类方法，首先对机械零部件现有不同状态的样本故障信号进行处理，建立不同故障类型的检验区间；其次获取机械零部件的所有单点故障状态所对应的故障特征频率，构造Duffing混沌振子的频率矩阵；然后求出不同故障特征频率下所对应的周期策动力幅值的临界阈值，构建频率-阈值矩阵；最后，将待检测信号加入计算最大Lyapunov指数矩阵M，根据M中数据进行检验，若存在故障信号，计算待检测信号的关联维数，对照已建立的不同故障类型的关联维数区间，进行故障分类，确定故障模式。本发明实现了对机械零部件早期单点故障的检测与分类，抗噪声能力强，并且故障检测成功率非常高。

Description

一种基于混沌的机械零部件早期单点故障检测与分类方法

技术领域

本发明涉及一种基于混沌的机械零部件早期单点故障检测与分类方法，属于机械零部件故障诊断技术领域。

背景技术

现代工业生产中，生产装置向大型化、复杂化、高速化、自动化和智能化的方向发展，不仅每一个设备的不同部分之间互相关联，紧密耦合，而且不同设备之间也存在着紧密的联系，在设备的运行过程中形成一个完整的系统。对于那些通常凭直观很难把握其运行状态的大型复杂机电设备而言，能否保证一些关键设备的正常运行，直接关系到一个企业发展的各个层面，轻者造成巨大的经济损失，重者还会产生严重的甚至灾难性的人员伤亡和社会影响。由于复杂先进的机电设备不应轻易解体检查，所以必须采用先进的测试设备和科学的方法。如何对机电设备进行不解体的监测与诊断，从获取的信息中分析设备的运行状态从而完成设备故障模式的识别，是当今机电设备状态监测与故障诊断的高技术高难点所在。

混沌理论(Chaos)是非线性科学的重要成就之一，与相对论、量子力学一起成为20世纪物理学的三次重大革命，它彻底消除了拉普拉斯关于决定论式可预测性的幻想。混沌现象是指在可确定的非线性系统中，系统随时间变化的运动状态对系统的初始条件非常敏感，并且形似紊乱，实则有序，无固定周期的循环性行为或形态。它是确定性非线性系统的内在随机性，这种随机性是由系统对初值的敏感依赖性而产生的；同时，它并非是杂乱无章，一片混乱，而是存在着复杂细致的几何结构，包含有更多的内在规律性。事实上，工程实际中的众多机械设备都具有混沌特性。对于这类混沌系统的故障检测，传统的解决方案通常都是把混沌背景当作噪声信号进行处理，很难对混沌系统本身的特性进行分析和利用，从而使检测效果无法令人满意。特别是在故障信号相对于系统混沌背景比较微弱的情况下，故障信号更是难于检测。

现有的机械零部件早期单点故障的故障诊断记录显示，机械零部件出现损伤的早期是诊断的最佳时期。机械零部件故障发生的早期，故障信号非常微弱，常被强烈的背景噪声所淹没，故要实现早期故障检测，实际上就是实现在强噪声背景下的微弱信号检测。长期以来，实现在强噪声背景下的微弱信号检测，应用最多的是频谱分析和小波分析方法。但是，这两种方法所能检测到的微弱信号的信噪比有限，当背景噪声比较强烈而所检测信号比较微弱时，它们不能很好地完成信号检测的任务。

发明内容

本发明的目的是为了解决在进行机械零部件早期单点故障的检测和分类时，现有方法检测成功率低、难以实现早期预报的问题，直接观察相轨迹的方法工作效率低、无法自动检测的问题，以及用Lyapunov指数方法进行故障分类时大量的对频率进行人工调整工作和专业性要求高的问题，结合Lyapunov指数和关联维数两种检测方法各自的优点，提出一种基于混沌的机械零部件早期单点故障检测与分类方法。

本发明一种基于混沌的机械零部件早期单点故障检测与分类方法，具体包括以下步骤：

步骤一、建立不同故障类型的关联维数区间。对机械零部件现有的、不同状态的样本故障信号，计算相应的关联维数，再应用基于小样本自采样方法对得到的关联维数进行自助训练，依据自助训练后得到的数值，进行正态分布样本均值与标准差的参数估计；根据所得到参数，建立不同故障类型的检验区间。

步骤二、构造Duffing混沌振子的频率矩阵。获取机械零部件的所有单点故障状态所对应的故障特征频率，建立包含所有状态的频率矩阵P，将混沌振子检测模型中策动力角频率ω设为频率矩阵P。所述的混沌振子检测模型为：

$\left\{\begin{array}{c}\stackrel{\′}{x}=y\\ \stackrel{\·}{y}=-0.5y+x^3-x^5+f\cos \mathrm{\ωt}+A\cos \mathrm{\ωt}+\mathrm{\σn}\left(t\right)\end{array}\right.$

其中，x、y为以时间t为自变量的函数，f为周期策动力幅值，ω为策动力角频率，n(t)为加性随机噪声，Acosωt+σn(t)为待检混合信号，σ为噪声平均功率σ²的正平方根。

步骤三、求出不同故障特征频率下所对应的周期策动力幅值f的临界阈值，构建频率-阈值矩阵。将步骤二中得到的混沌振子检测模型在不加入外部信号Acosωt+σn(t)＝0情况下，调节周期策动力幅值f，使混沌振子检测模型系统处于临界的混沌态，将Lyapunov指数曲线的最后一个过零点所对应的周期策动力幅值f作为临界阈值f_d，并与步骤二获取的故障特征频率建立相对应的频率-阈值矩阵。

步骤四、进行故障检测。将频率-阈值矩阵中的每对对应的值都分别代入混沌振子检测模型中周期策动力角频率ω与策动力幅值f，得到一个混沌振子检测的方程组。然后向此方程组加入待检测信号，此时Acosωt+σn(t)的值为待检验信号值，计算最大Lyapunov指数，将得到的所有最大Lyapunov指数组成最大Lyapunov指数矩阵M，判断最大Lyapunov指数矩阵M中数据是否全部大于零，若是则无故障信号存在，结束本次故障检测与分类过程；若不是全部大于零则存在故障，执行步骤五。

步骤五、进行故障分类。计算待检测信号的关联维数，对照步骤一建立的不同故障类型的关联维数区间，进行故障分类，确定故障模式。

本发明的优点与积极效果在于：

(1)充分利用了Duffing混沌振子具有对某些参数变动非常敏感的特性，成功的检测到强噪声背景下的微弱故障信号，实现了对机械零部件早期单点故障的检测与分类，抗噪声能力强，故障检测成功率非常高，效果显著；

(2)引入Lyaponov指数，克服了相轨迹图法的工作效率低、主观因素大等缺点，实现了对故障的自动检测；

(3)界定了最大Lyapunov指数的最后一个过零点来确定策动力幅值f的临界阈值，克服了混沌区到大尺度周期区之间过渡区对临界阈值选取的负面影响；

(3)采用基于关联维数的故障分类方法，解决了用Lyapunov指数方法进行故障分类时大量的对频率进行人工调整工作和专业性要求高的问题；

(4)在计算关联维数的故障检测区间时，应用基于小样本自采样方法进行自助训练，成功的通过小样本来确定具备较高可信度的检验区间；

(5)本发明方法利用样本，无需建立模型，即可实现故障的检测和识别，降低了专业要求，增加了工程应用性；

(6)与现有的故障检测和识别方法相比，显著提高了通用性和精度。

附图说明

图1是本发明故障检测与分类方法的整体步骤流程图；

图2是本发明故障检测与分类方法步骤一的步骤流程图；

图3是本发明实施例中样本数据的关联维数分布图；

图4是经过自助训练后所得的样本均值的正态分布检验图；

图5是经过自助训练后的样本标准差的正态分布检验图；

图6是本发明实施例关联维数检验区间图示意图；

图7是周期策动力幅值f与最大Lyapunov指数的关系。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明结合Lyapunov指数和关联维数两种检测方法各自的优点，提出一种基于混沌的机械零部件早期单点故障检测与分类方法。

基于混沌抑制的弱信号检测方法是混沌理论在信号分析上的一个重要分支，混沌抑制的方法很多，但实际应用中多集中在基于Holmes型Duffing振子的检测上，Holmes型Duffing方程适合于检测任意频率的微弱周期信号，对噪声具有一定程度的免疫力，而对与内部周期摄动力同频的周期信号具有相对较高的敏感性，检测性能达到了很低的信噪比。

Duffing系统所描述的非线性动力学系统表现出丰富的非线性动力学特性，包括振荡、分岔、混沌的复杂动态，已成为研究混沌的常用模型之一。Duffing方程具体形式为：

式(1)中k为阻尼比；f为周期策动力幅值；ω为策动力角频率；ax(t)+bx³(t)项为非线性恢复力，a、b为实数因子，函数x(t)以时间t为自变量。

Duffing系统是一个非线性动力系统，某些系数的摄动会引起其解的性态发生本质的变化。检测前调节系统参数使系统处于某种状态，把外加待测信号作为系统某种作用力的补充，改变了系统的参数，使系统输出时域波形或相图发生某种非常明显的变化，例如系统状态由周期状态变成混沌状态，从而检测出微弱信号。

利用Duffing振子检测微弱信号的方程为：

$\stackrel{\·\·}{x}\left(t\right)+k\stackrel{\·}{x}\left(t\right)-(x^3-x^5)=f\cos \mathrm{\ωt}+\mathrm{As}\left(t\right)+\mathrm{\σn}\left(t\right)---\left(2\right)$

式(2)中As(t)为待检测信号，n(t)为加性随机噪声，σ为噪声平均功率σ²的正平方根，x³-x⁵为非线性恢复力；fcosωt为周期策动力，系数k＝0.5。其等价系统为：

$\left\{\begin{array}{c}\stackrel{\′}{x}=y\\ \stackrel{\·}{y}=-0.5y+x^3-x^5+f\cos \mathrm{\ωt}+A\cos \mathrm{\ωt}+\mathrm{\σn}\left(t\right)\end{array}\right.---\left(3\right)$

此式中，待检混合信号为Acosωt+σn(t)。

本发明的单点故障检测与分类方法主要依据混沌振子检测模型的式(3)来作为解决问题的理论依据。

本发明是一种针对机械零部件早期单点故障的采用Lyapunov指数和关联维数相结合的故障检测与分类的方法，如图1所示，具体步骤如下：

步骤一、建立不同故障类型的检验区间。

具体建立检验区间如图2所示，对机械零部件现有的、不同状态的样本故障数据，采用G-P算法计算相应的关联维数，再应用基于小样本自采样方法对关联维数进行自助训练，依据自助训练后得到的数值，进行正态分布的参数估计，所述的正态分布的参数为样本均值与标准差；根据所得到参数，计算出样本置信度为95％的分布区间，以此分布区间为不同故障类型的检验区间，为后面的故障分类提供依据。

所述的小样本自采样方法即自助法，并且在进行自助训练前，根据统计学原理检验样本分布是否是正态分布。一般情况下，机械零部件现有的、不同状态的样本故障数据都符合正态分布。

在不同运行状态下，因刚度非线性摩擦力等的影响，滚动轴承系统表现出不同的非线性特性，分形维数是用来定量刻画混沌吸引子“奇异”程度的一个重要参数，也可以用来刻画滚动轴承的故障状态。分形的维数有许多种类，最具代表性的是关联维数。关联维数能够体现未知系统的固有特性，并且机制不同的各类故障通常也具有不同的关联维数，可以作为系统故障特征量检测并区分滚动轴承的故障状态与故障模式。相同工作状态下的信号具有相近的关联维数，不同故障模式下的关联维数具有不同的数值，存在明显的可分性。而且，计算关联维数不需要建立系统方程，仅依靠一段样本数据就可以计算出该样本的关联维数，具备很强的通用性。实际应用中，样本数据往往比较匮乏，还需要采用自助法对样本进行重采样，将小样本问题转换为大样本问题来模仿未知分布。

自助法的基本原理是：在总体中抽取M个样本构成初始样本，之后随机、等概、独立、有放回地抽取M个样本单元，构成一个新的点集，即一个自助样本。数学描述为：设随机样本X＝(x₁，x₂…x_n)是来自于某未知的总体分布F(x)，θ＝θ(F(x))为总体分布F的某个未知参数，F_n(x)为抽样分布函数，

为θ的估计。记估计误差为：

记为从F_n(x)中抽样获得的再生样本，

是由X^*所获得的抽样分布函数。记：

$R_n^{*}=\widehat{\mathrm{\θ}}\left(F_n^{*}\left(x\right)\right)-\widehat{\mathrm{\θ}}\left(F_n\left(x\right)\right)---\left(5\right)$

称为T_n的Bootstrap统计量。在给定抽样分布函数F_n(x)的条件下，取所有统计量

的均值

去模仿估计误差T_n，则总体分布F(x)的参数

关联维数的分布通常符合正态分布，故而总体分布的参数有两个，分别为样本均值与标准差。根据自助训练得到的参数θ可以进一步得到样本均值与标准差。

目前计算关联维数的方法最主要的就是G-P算法，G-P算法是Grassberger和Procassi于1983年提出的一种比较容易实现的从实验数据中估算关联维数的算法。G-P算法如下：

对于时间序列{x(i)|i＝1，2，…n-l，n}，n是序列的长度，在本发明中，时间序列中的各值就是机械零部件现有的、不同状态的样本故障数据。

首先，将时间序列嵌入到m维欧氏空间R^m，得到n_m个样本点，这n_m个样本点用{y(i)|i＝1，2，…n_m-1，n_m}表示，其中n_m＝n-(m-1)τ，其中τ为延迟时间，m为嵌入维数，然后计算关联积分C(m，n，r，t)：

$C(m,n,r,\mathrm{\τ})=\frac{2}{n_m(n_m-1)}\underset{i=1}{\overset{n_m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n_m}{\mathrm{\Σ}}}H(r-D(i,j))---\left(6\right)$

r为相空间超球半径，H(r-D(i，j))根据Heaviside函数H(x)计算：

$H\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.---\left(7\right)$

m维欧氏空间R^m中的欧式距离为：

$\left|\right|y_i-y_j\left|\right|={\left[\underset{k=1}{\overset{m-1}{\mathrm{\&Sigma;}}}{(x_{i+\mathrm{k\&tau;}}-x_{j+\mathrm{k\&tau;}})}^2\right]}^{\frac{1}{2}}---\left(8\right)$

y_i、y_j表示n_m个样本点中的样本点y(i)与y(j)。对于给定的时间序列，r＞0足够小时，可定义关联维数：

$d(m,\mathrm{\&tau;})=\underset{r\&RightArrow;0}{\lim }\frac{\ln C_m\left(r\right)}{\ln r}---\left(9\right)$

画出标度曲线lnC-lnr图，取标度线中的线性区域的斜率来近似代替这个关联维数：

$d(m,\mathrm{\&tau;})=\frac{\ln C_m\left(r\right)}{\ln r}---\left(10\right)$

式(10)就是采用G-P算法最终得到的关联维数。

步骤二、构造Duffing混沌振子的频率矩阵。

获取机械零部件的所有单点故障状态所对应的故障特征频率，建立包含所有状态的频率矩阵P，将如式(3)中的混沌振子检测模型策动力角频率ω设为频率矩阵P。

步骤三、求出不同故障特征频率下所对应的f的临界阈值，构建频率-阈值矩阵。

将如式(3)中的混沌振子检测模型中策动力角频率设为频率矩阵P，不加入外部信号，即Acosωt+σn(t)＝0时，调节周期策动力幅值f，使如式(3)中的混沌振子检测模型系统处于临界的混沌态，如图7所示，为周期策动力幅值f与最大Lyapunov指数的关系，图7中，周期策动力幅值f的取值范围为[0.71，0.75]，步长step即两个连续点之间的距离为0.0001，将Lyapunov指数曲线的最后一个过零点所对应的策动力幅值f作为系统的临界阈值fd，并建立与频率相对应的频率-阈值矩阵。

在未加待检混合信号前，即Acosωt+σn(t)＝0时，调节周期策动力幅值f，使系统处于临界的混沌态。然后注入待检混合信号，

步骤四、进行故障检测。将频率-阈值矩阵中的每对对应的值都分别代入如式(3)中的混沌振子检测模型中周期策动力角频率ω与策动力幅值f，得到一个混沌振子检测的方程组。此时方程组的每个模型都处于临界的混沌态，然后向此方程组加入待检测信号，此时Acosωt+σn(t)的值为待检验信号值，如果待检信号中含有一定幅值的弱正弦信号，则会使得系统由临界的混沌态转为大尺度周期态。

采用直接观察相轨迹变化的方法，可以判断系统是否为混沌，但却需要人的肉眼观察，工作效率低。为了刻画系统是否处于混沌状态，引入Lyapunov指数。Lyapunov指数用来度量在相空间中初始条件不同的两条相邻轨迹随时间按指数率收敛或发散的程度，这种轨迹收敛或发散的比率，称为Lyapunov指数。在混沌状态判据中，Lyapunov指数起着非常重要的作用。一个系统是否是混沌的，可以由它的Lyapunov指数是否有正值来确定，这种方法比别的方法更准确，它给出了一个定量分析的标准。Lyapunov指数法确定阈值的基本思想是：通过计算系统的最大Lyapunov指数，根据它的符号变换(由正变负)，来确定系统由混沌状态变为周期状态的临界阈值，将系统置于临界状态，加入待检测信号，若最大Lyapunov指数由大于零变为小于零则证明该故障信号的存在。以此可以进行自动地故障检测。

计算最大Lyapunov指数，将得到的所有最大Lyapunov指数组成最大Lyapunov指数矩阵M，判断最大Lyapunov指数矩阵M中数据是否全部大于零，若M中数据全部大于零则判定无故障信号存在，若M中矩阵数据不是全部大于零则判定有故障信号，进入基于关联维数的故障分类计算。

Lyapunov指数的定义，是基于自治系统，而不是非自治系统的。因此，我们在求Lyapunov指数的时候，都需要将原来的非自治系统转化为自治系统，再进行求解。Holmes型简化的Duffing方程(3)为二维非自治系统，通过代换Z＝t，可以写成相应的三维自治系统：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=y\\ \stackrel{\&CenterDot;}{y}=-0.5y+x^3-x^5+f\cos z\\ \stackrel{\&CenterDot;}{z}=1\end{array}\right.---\left(11\right)$

针对混沌系统，存在定理1：非自治系统的N个Lyapunov指数等于其相应自治系统的前N个Lyapunov指数，且该自治系统的第N+1个Lyapunov指数等于零。

根据定理1，非自治系统(3)的两个Lyapunov指数等于自治系统(11)的前两个Lyapunov指数，因此现在只需考虑自治系统的Lyapunov指数求解问题。

对于自治系统(11)的线性变分方程为

$\stackrel{\&CenterDot;}{Y}\left(t\right)=J\left(t\right)Y\left(t\right),Y\left(0\right)=I_3---\left(12\right)$

Y∈R^3×3，I₃是3×3的单位矩阵；J(t)为三维自治系统(11)的雅可比(Jacobi)矩阵，可以看出，Jacobi矩阵是奇异的，因为它的第三行元素都为零。这意味着Y(t)的第三行元素都为常量。再利用初始条件Y(0)＝I₃，知Y(t)的第三行元素分别为Y₃₁＝0，Y₃₂＝0以及Y₃₃＝1。由于Y是非奇异的，将它进行QR分解，记为Y(t)＝Q(t)R(t)，R(t)为上三角矩阵，其对角元素均为正数。得到表达式(13)：

$Y\left(t\right)=\left[\begin{array}{ccc}{Y}_{11}& Y_{12}& Y_{13}\\ Y_{21}& Y_{22}& Y_{23}\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ 0& R_{22}& R_{23}\\ 0& 0& 1\end{array}\right]---\left(13\right)$

根据前面Lyapunov指数的定义，若δ_i为系统雅可比矩阵的第i个特征值，则Lyapunov指数

因为R₃₃≡1，可知三维自治系统(11)的一个Lyapunov指数始终为零，这也反映了定理1的内容。那么另外Lyapunov指数的计算只需求取R₁₁(t)和R₂₂(t)。这样就仅考虑三维自治系统(11)的二维自治子系统就可以了。二维自治子系统的变分方程可以表示为：

$\stackrel{\stackrel{\&CenterDot;}{~}}{Y}\left(t\right)=\tilde{J}\left(t\right)\tilde{Y}\left(t\right),\tilde{Y}\left(0\right)=I_2---\left(14\right)$

同理，

I₂为2×2的单位矩阵。为雅可比矩阵，其形式为：

$\tilde{J}\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\frac{\&PartialD;(\mathrm{ky}+x^3-x^5+f\cos z)}{\&PartialD;x}& -\frac{\&PartialD;(\mathrm{ky}+x^3-x^5+f\cos z)}{\&PartialD;y}\end{array}\right]---\left(15\right)$

将式(15)简记为

$\tilde{J}\left(t\right)=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]---\left(16\right)$

将

进行QR分解，即为

代入式(14)有

$\stackrel{\stackrel{\&CenterDot;}{~}}{Q}\tilde{R}\left(t\right)+\tilde{Q}\stackrel{\stackrel{\&CenterDot;}{~}}{R}=\tilde{J}\tilde{Q}\tilde{R},\tilde{Q}\left(0\right)\tilde{R}\left(0\right)=I_2---\left(17\right)$

再将式(17)左乘

右乘

得到：

${\tilde{Q}}^T\stackrel{\stackrel{\&CenterDot;}{~}}{Q}+\stackrel{\stackrel{\&CenterDot;}{~}}{R}{\tilde{R}}^{-1}={\tilde{Q}}^T\tilde{J}\tilde{Q},\tilde{Q}\left(0\right)=I_2,\tilde{R}\left(0\right)=I_2---\left(18\right)$

由于

矩阵是斜对称矩阵，并且由于矩阵是上三角矩阵，所以

也是上三角矩阵，因此，矩阵K(t)中的元素K_ij(t)为：

$K_{\mathrm{ij}}\left(t\right)=\left\{\begin{array}{cc}{\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ij}},& i>j\\ 0& i=j\\ -{\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ji}}& i<j\end{array}\right.---\left(19\right)$

从中可以得到

$\stackrel{\stackrel{\&CenterDot;}{~}}{Q}\left(t\right)=\tilde{Q}\left(t\right)K\left(t\right),\tilde{Q}\left(0\right)=I_2---\left(20\right)$

由式(19)和式(20)就可以确定了关于矩阵

的微分方程。由于

是斜对称的，由式(18)得到

${\stackrel{\stackrel{\&CenterDot;}{~}}{R}}_i^{\&prime;}\left(t\right)=\frac{{\stackrel{\&CenterDot;}{R}}_{\mathrm{ii}}\left(t\right)}{R_{\mathrm{ii}}\left(t\right)}={\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ij}},R_i^{\&prime;}=0,i=\mathrm{1,2}---\left(21\right)$

这里R′_i(t)＝ln(R_ii(t))，Lyapunov指数的时间演化式为σ_i(t)＝R′_i(t)/t，所以可以得到Lyapunov指数分别为

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1=\underset{t\&RightArrow;\&infin;}{\lim }\frac{R_1^{\&prime;}\left(t\right)}{t}=\underset{t\&RightArrow;\&infin;}{\lim }{\mathrm{\&sigma;}}_1\left(t\right)\\ {\mathrm{\&sigma;}}_2=\underset{t\&RightArrow;\&infin;}{\lim }\frac{R_2^{\&prime;}\left(t\right)}{t}=\underset{t\&RightArrow;\&infin;}{\lim }{\mathrm{\&sigma;}}_2\left(t\right)\end{array}\right.---\left(22\right)$

σ₁和σ₂就是自治系统(11)的另两个Lyapunov指数，也即Duffing方程式(3)的两个Lyapunov指数。

将如式(3)中的混沌振子检测模型中周期策动力幅值f与策动力角频率ω每代入频率-阈值矩阵中的两个对应值，根据上面计算Lyapunov指数的方法得到两个Lyapunov指数，选择最大的Lyapunov指数输出为最大Lyapunov指数矩阵M中一个值，最终得到完整的最大Lyapunov指数矩阵M。

步骤五、进行故障分类。若待检测信号含有故障信号，计算其关联维数。对照步骤一建立的不同故障类型的关联维数区间，进行故障分类，确定故障模式。关联维数的计算同步骤一中计算关联维数的方法相同。

实施例：

本实例采取华盛顿天主教大学轴承数据中心提供的滚动轴承故障信号进行验证。分别使用正常、内环故障、外环故障和滚动体故障四种状态下的样本信号对本发明针对机械零部件早期单点故障的Lyapunov指数和关联维数相结合的故障检测和分类的方法进行检测验证，具体步骤如下：

步骤一、建立不同故障类型的关联维数区间。

根据已有的滚动轴承实际故障数据进行关联维数的计算，如表1所示。

表1滚动轴承不同状态振动信号的关联维数

信号样本编号	正常	滚动体故障	内环故障	外环故障
					1	1.620	2.509	3.568	3.370
2	1.607	2.593	3.722	3.163
					3	1.666	2.693	3.641	3.139
4	1.667	2.637	3.692	3.274
					5	1.607	2.494	3.635	3.197
6	1.607	2.555	3.503	3.153
					7	1.683	2.661	3.541	3.215
8	1.673	2.709	3.493	3.253
					9	1.567	2.431	3.666	3.098
10	1.579	2.483	3.638	3.242
					均值	1.6166	2.5122	3.6096	3.2104

对机械零部件现有的、不同状态的样本故障信号，计算相应的关联维数，如图3所示；根据统计学原理，对样本进行正态性检验；再应用基于小样本自采样方法进行自助训练，例如内环故障信号样本的均值和标准差分别如图4与图5所示；然后进行正态分布的参数估计，得到；根据所得到样本均值与标准差的估计值，计算出样本置信度为95％的分布区间。如图6所示为数据经过自助训练后分别计算出的检验区间图，分别为滚动体故障[2.3968，2.7567]，外环故障[3.0674，3.3534]，内环故障[3.4643，3.7571]。

步骤二、构造Duffing混沌振子的频率矩阵

根据滚动轴承的尺寸及转速可以直接求出滚动轴承的各种故障特征频率。表2为华盛顿天主教大学轴承数据中心提供的用于计算滚动轴承的单点故障特征频率的实验值。

表2用于计算滚动轴承单点故障特征频率的实验值表

内环故障	外环故障	滚动体故障
			5.4152	3.5848	4.7135

用采集信号时的轴的旋转角速度，乘以表2中的数值，即可计算出滚动轴承各种故障状态下对应的故障特征频率，进一步建立包含所有状态的频率矩阵P。

步骤三、求出不同故障特征频率下所对应的周期策动力幅值的临界阈值f_d，构建频率-阈值矩阵。在实际应用中，一般将如式(3)中的混沌振子检测模型做时间尺度上的变换，令t＝ωτ，可得到针对任意频率正弦信号的检测数学模型为：

$\left\{\begin{array}{c}\stackrel{\&prime;}{x}=\mathrm{\&omega;y}\\ \stackrel{\&CenterDot;}{y}=\mathrm{\&omega;}(-0.5y+x^3-x^5+f\cos \mathrm{\&omega;t}+A\cos \mathrm{\&omega;t}+\mathrm{\&sigma;n}\left(t\right))\end{array}\right.---\left(23\right)$

将策动力角频率ω分别设为频率矩阵P中的数值，不加入外部信号，调节f，使如式(3)中的混沌振子检测模型系统处于临界的混沌态，将Lyapunov指数曲线的最后一个过零点所对应的策动力幅值f作为系统的临界阈值f_d，并建立与频率相对应的频率-阈值矩阵。

步骤四、进行故障检测。

将如式(3)中的混沌振子检测模型中策动力幅值与角频率分别设为频率-阈值矩阵中的各个值；然后加入待检测信号计算最大Lyapunov指数，输出最大Lyapunov指数矩阵M；之后进行检验：若M中数据全部大于零则判定无故障信号存在，若M中矩阵数据不是全部大于零则判定有故障信号，进入基于关联维数的故障分类计算。

步骤五、进行故障分类。

若待检测信号含有故障信号，计算其关联维数，对照步骤一建立的不同故障类型的关联维数区间，如图6所示，关联维数落在哪一个区间，即为哪一种故障类型。

采取华盛顿天主教大学轴承数据中心提供的滚动轴承故障信号进行验证，分别采取正常、内环故障、外环故障和滚动体故障四种状态下20组样本信号进行检测，检测结果如表3：

表3检测结果

样本结果	正确	误判	漏判	检测合格率
					正常轴承信号	20/20	20/0	0	100％
内环故障	20/19	20/1	0	95％

外环故障	20/19	20/1	0	95％
					滚动体故障	20/20	20/0	0	100％

通过表3可以看出正常轴承信号和滚动体故障的判断全部正确，合格率为100％，内环故障有一个误判为外环故障，外环故障同样有一个误判为内环故障，这可能是两种故障的检验区间比较接近的缘故，其检测合格率为95％，可见本发明的方法能够实现机械零部件早期单点故障检测与分类，故障检测成功率高，具有明显的实际应用价值。

Claims (5)

1.一种基于混沌的机械零部件早期单点故障检测与分类方法，其特征在于，具体包括步骤：

步骤一、建立不同故障类型的检验区间：对机械零部件现有的、不同状态的样本故障信号，计算相应的关联维数，对结果进行正态性检验，再应用基于小样本自采样方法对得到的关联维数进行自助训练；依据自助训练后得到的数值，进行正态分布的样本均值与标准差的参数估计；根据所得到参数，建立不同故障类型的检验区间；

步骤二、构造Duffing混沌振子的频率矩阵：获取机械零部件的所有单点故障状态所对应的故障特征频率，建立包含所有单点故障状态的频率矩阵P，将混沌振子检测模型中策动力角频率ω设为频率矩阵P；所述的混沌振子检测模型为：

$\left\{\begin{array}{c}\stackrel{\&prime;}{x}=y\\ \stackrel{\&CenterDot;}{y}=-0.5y+x^3-x^5+f\cos \mathrm{\&omega;t}+A\cos \mathrm{\&omega;t}+\mathrm{\&sigma;n}\left(t\right)\end{array}\right.$

其中，x、y为以时间t为自变量的函数，f为周期策动力幅值，ω为策动力角频率，n(t)为加性随机噪声，Acosωt+σn(t)为待检混合信号；

步骤三、求出不同故障特征频率下所对应的周期策动力幅值f的临界阈值f_d，构建频率-阈值矩阵：将步骤二中得到的混沌振子检测模型在不加入待检测信号情况下，此时Acosωtσn(t)＝0，调节周期策动力幅值f，使混沌振子检测模型系统处于临界的混沌态，将Lyapunov指数曲线的最后一个过零点所对应的周期策动力幅值f作为系统的临界阈值f_d，并建立与步骤二获取的故障特征频率相对应的频率-阈值矩阵；

步骤四、进行故障检测：将频率-阈值矩阵中的每对对应的值都分别代入混沌振子检测模型中周期策动力幅值f与策动力角频率ω，然后在每一次代入值后加入待检测信号，此时Acosωt+σn(t)的值为待检验信号值，计算最大Lyapunov指数，将得到的所有最大Lyapunov指数组成最大Lyapunov指数矩阵M，判断最大Lyapunov指数矩阵M中数据是否全部大于零，若是则无故障信号存在，结束本次故障检测与分类过程；若不是全部大于零则存在故障，执行步骤五；

步骤五、进行故障分类：计算待检测信号的关联维数，对照步骤一建立的不同故障类型的检验区间，进行故障分类，确定故障模式。

2.根据权利要求1所述的一种基于混沌的机械零部件早期单点故障检测与分类方法，其特征在于，所述步骤一与步骤五中的关联维数采用G-P算法计算得到：

首先，将时间序列{x(i)|i＝1，2，…n-l，n}嵌入到m维欧氏空间R^m，得到n_m个样本点，这n_m个样本点用{y(i)|i＝1，2，…n_m-1，n_m}表示，其中n_m＝n-(m-1)τ，τ为延迟时间，m为嵌入维数，然后计算关联积分C(m，n，r，t)：

$C(m,n,r,\mathrm{\&tau;})=\frac{2}{n_m(n_m-1)}\underset{i=1}{\overset{n_m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n_m}{\mathrm{\&Sigma;}}}H(r-D(i,j))$

r为相空间超球半径，H(r-D(i，j))根据Heaviside函数

计算得到，时间序列中x(i)的值，在步骤一中，对应机械零部件现有的、不同状态的样本故障信号数据，在步骤五中，对应待检测信号数据，n为时间序列的长度；

m维欧氏空间R^m中，D(i，j)为欧式距离：

其中，y_i、y_j表示n_m个样本点中的样本点y(i)与y(j)，对于给定的时间序列，r＞0足够小时，有关联维数：

$d(m,\mathrm{\&tau;})=\underset{r\&RightArrow;0}{\lim }\frac{\ln C_m\left(r\right)}{\ln r}$

C_m(r)为关联积分C(m，n，r，t)的简写，画出标度曲线lnC-lnr图，取标度线中的线性区域的斜率来近似代替这个关联维数：

$d(m,\mathrm{\&tau;})=\frac{\ln C_m\left(r\right)}{\ln r}$

所得到的这个关联维数就是要求的关联维数。

3.根据权利要求1所述的一种基于混沌的机械零部件早期单点故障检测与分类方法，其特征在于，步骤一所述的小样本自采样方法具体为：

将机械零部件现有的、同种故障类型的信号所计算出来的关联维数作为初始样本X＝(x₁，x₂…x_n)，对初始样本，随机、等概、独立、有放回地抽取M个样本单元，构成一个新的点集，形成一个自助样本；该初始样本为该故障类型下关联维数的总体分布F(x)的一个随机样本，总体分布F(x)的未知参数θ＝θ(F(x))，参数θ的估计

F_n(x)为抽样分布函数，则估计误差

从抽样分布函数F_n(x)中抽样获得再生样本X^*＝(x₁ ^*，x₂ ^*…x_n)^*，则估计误差T_n的Bootstrap统计量为：

是由X^*所获得的抽样分布函数；

在给定抽样分布函数F_n(x)的条件下，取所有统计量

的均值

去模仿估计误差T_n，则总体分布F(x)的参数

4.根据权利要求1所述的一种基于混沌的机械零部件早期单点故障检测与分类方法，其特征在于，步骤一所述的检验区间，是采用样本置信度为95％的分布区间得到的。

5.根据权利要求1所述的一种基于混沌的机械零部件早期单点故障检测与分类方法，其特征在于，步骤四中的最大Lyapunov指数具体是根据以下过程得到的：

首先通过代换Z＝t，将混沌振子检测模型等价为相应的三维自治系统：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=y\\ \stackrel{\&CenterDot;}{y}=-0.5y+x^3-x^5+f\cos z\\ \stackrel{\&CenterDot;}{z}=1\end{array}\right.$

然后求解三维自治系统的Lyapunov指数

对于三维自治系统的线性变分方程：

Y(0)＝I₃，Y(t)的第三行元素分别为Y₃₁＝0，Y₃₂＝0以及Y₃₃＝1，其中，Y∈R^3×3，I₃是3×3的单位矩阵，J(t)为三维自治系统的雅可比矩阵；将Y(t)进行QR分解，有Y(t)＝Q(t)R(t)，R(t)为上三角矩阵，其对角元素均为正数，得到：

$Y\left(t\right)=\left[\begin{array}{ccc}{Y}_{11}& Y_{12}& Y_{13}\\ Y_{21}& Y_{22}& Y_{23}\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ 0& R_{22}& R_{23}\\ 0& 0& 1\end{array}\right]$

Lyapunov指数

δ_i为系统雅可比矩阵的第i个特征值，因为R₃₃≡1，则三维自治系统的一个Lyapunov指数始终为零；根据三维自治系统的二维自治子系统求取R₁₁(t)和R₂₂(t)，进一步求取另外两个Lyapunov指数：

三维自治系统的二维自治子系统：

I₂为2×2的单位矩阵，

为雅可比矩阵：

$\tilde{J}\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\frac{\&PartialD;(\mathrm{ky}+x^3-x^5+f\cos z)}{\&PartialD;x}& -\frac{\&PartialD;(\mathrm{ky}+x^3-x^5+f\cos z)}{\&PartialD;y}\end{array}\right]$

将

进行QR分解：

代入二维自治子系统，得到：

$\stackrel{\stackrel{\&CenterDot;}{~}}{Q}\tilde{R}\left(t\right)+\tilde{Q}\stackrel{\stackrel{\&CenterDot;}{~}}{R}=\tilde{J}\tilde{Q}\tilde{R},\tilde{Q}\left(0\right)\tilde{R}\left(0\right)=I_2$

再将上式左乘

右乘

得到：

${\tilde{Q}}^T\stackrel{\stackrel{\&CenterDot;}{~}}{Q}+\stackrel{\stackrel{\&CenterDot;}{~}}{R}{\tilde{R}}^{-1}={\tilde{Q}}^T\tilde{J}\tilde{Q},$ $\tilde{Q}\left(0\right)=I_2,$ $\tilde{R}\left(0\right)=I_2$

由于

矩阵

是斜对称矩阵，矩阵

是上三角矩阵，所以

也是上三角矩阵，因此，矩阵K(t)中的元素K_ij(t)为：

$K_{\mathrm{ij}}\left(t\right)=\left\{\begin{array}{cc}{\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ij}},& i>j\\ 0& i=j\\ -{\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ji}}& i<j\end{array}\right.$

从

中得到：

$\stackrel{\stackrel{\&CenterDot;}{~}}{Q}\left(t\right)=\tilde{Q}\left(t\right)K\left(t\right),\tilde{Q}\left(0\right)=I_2$

由于

是斜对称的，于是得到：

${\stackrel{\stackrel{\&CenterDot;}{~}}{R}}_i^{\&prime;}\left(t\right)=\frac{{\stackrel{\&CenterDot;}{R}}_{\mathrm{ii}}\left(t\right)}{R_{\mathrm{ii}}\left(t\right)}={\left({\tilde{Q}}^T\left(t\right)\tilde{J}\left(t\right)\tilde{Q}\left(t\right)\right)}_{\mathrm{ij}},R_i^{\&prime;}=0,i=\mathrm{1,2}$

此处Lyapunov指数的时间演化式为σ_i(t)＝R′_i(t)/t，所以得到Lyapunov指数分别为：

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_1=\underset{t\&RightArrow;\&infin;}{\lim }\frac{R_1^{\&prime;}\left(t\right)}{t}=\underset{t\&RightArrow;\&infin;}{\lim }{\mathrm{\&sigma;}}_1\left(t\right)\\ {\mathrm{\&sigma;}}_2=\underset{t\&RightArrow;\&infin;}{\lim }\frac{R_2^{\&prime;}\left(t\right)}{t}=\underset{t\&RightArrow;\&infin;}{\lim }{\mathrm{\&sigma;}}_2\left(t\right)\end{array}\right.$

σ₁和σ₂就是步骤二中混沌振子检测模型所要求取的两个Lyapunov指数；选取这两个Lyapunov指数中的大者，作为此次周期策动力幅值f与策动力角频率ω代入值后得到的混沌振子检测模型的最大Lyapunov指数。

Images