本发明涉及一种基于混沌的机械零部件早期单点故障检测与分类方法，属于机械零部件故障诊断技术领域。
背景技术：
现有的机械零部件早期单点故障的故障诊断记录显示，机械零部件出现损伤的早期是诊断的最佳时期。机械零部件故障发生的早期，故障信号非常微弱，常被强烈的背景噪声所淹没，故要实现早期故障检测，实际上就是实现在强噪声背景下的微弱信号检测。长期以来，实现在强噪声背景下的微弱信号检测，应用最多的是频谱分析和小波分析方法。但是，这两种方法所能检测到的微弱信号的信噪比有限，当背景噪声比较强烈而所检测信号比较微弱时，它们不能很好地完成信号检测的任务。
技术实现要素：
本发明的目的是为了解决在进行机械零部件早期单点故障的检测和分类时，现有方法检测成功率低、难以实现早期预报的问题，直接观察相轨迹的方法工作效率低、无法自动检测的问题，以及用lyapunov指数方法进行故障分类时大量的对频率进行人工调整工作和专业性要求高的问题，结合lyapunov指数和关联维数两种检测方法各自的优点，提出一种基于混沌的机械零部件早期单点故障检测与分类方法。
本发明一种基于混沌的机械零部件早期单点故障检测与分类方法，具体包括以下步骤：
步骤一、建立不同故障类型的关联维数区间。对机械零部件现有的、不同状态的样本故障信号，计算相应的关联维数，再应用基于小样本自采样方法对得到的关联维数进行自助训练，依据自助训练后得到的数值，进行正态分布样本均值与标准差的参数估计；根据所得到参数，建立不同故障类型的检验区间。
步骤二、构造duffing混沌振子的频率矩阵。获取机械零部件的所有单点故障状态所对应的故障特征频率，建立包含所有状态的频率矩阵p，将混沌振子检测模型中策动力角频率ω设为频率矩阵p。所述的混沌振子检测模型为：
x′=yy·=-0.5y+x3-x5+fcosωt+acosωt+σn(t)]]>
步骤三、求出不同故障特征频率下所对应的周期策动力幅值f的临界阈值，构建频率-阈值矩阵。将步骤二中得到的混沌振子检测模型在不加入外部信号acosωt+σn(t)＝0情况下，调节周期策动力幅值f，使混沌振子检测模型系统处于临界的混沌态，将lyapunov指数曲线的最后一个过零点所对应的周期策动力幅值f作为临界阈值fd，并与步骤二获取的故障特征频率建立相对应的频率-阈值矩阵。
步骤四、进行故障检测。将频率-阈值矩阵中的每对对应的值都分别代入混沌振子检测模型中周期策动力角频率ω与策动力幅值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方程具体形式为：
duffing系统是一个非线性动力系统，某些系数的摄动会引起其解的性态发生本质的变化。检测前调节系统参数使系统处于某种状态，把外加待测信号作为系统某种作用力的补充，改变了系统的参数，使系统输出时域波形或相图发生某种非常明显的变化，例如系统状态由周期状态变成混沌状态，从而检测出微弱信号。
利用duffing振子检测微弱信号的方程为：
x··(t)+kx·(t)-(x3-x5)=fcosωt+as(t)+σn(t)---(2)]]>
式(2)中as(t)为待检测信号，n(t)为加性随机噪声，σ为噪声平均功率σ2的正平方根，x3-x5为非线性恢复力；fcosωt为周期策动力，系数k＝0.5。其等价系统为：
x′=yy·=-0.5y+x3-x5+fcosωt+acosωt+σn(t)---(3)]]>
此式中，待检混合信号为acosωt+σn(t)。
本发明的单点故障检测与分类方法主要依据混沌振子检测模型的式(3)来作为解决问题的理论依据。
本发明是一种针对机械零部件早期单点故障的采用lyapunov指数和关联维数相结合的故障检测与分类的方法，如图1所示，具体步骤如下：
步骤一、建立不同故障类型的检验区间。
具体建立检验区间如图2所示，对机械零部件现有的、不同状态的样本故障数据，采用g-p算法计算相应的关联维数，再应用基于小样本自采样方法对关联维数进行自助训练，依据自助训练后得到的数值，进行正态分布的参数估计，所述的正态分布的参数为样本均值与标准差；根据所得到参数，计算出样本置信度为95％的分布区间，以此分布区间为不同故障类型的检验区间，为后面的故障分类提供依据。
所述的小样本自采样方法即自助法，并且在进行自助训练前，根据统计学原理检验样本分布是否是正态分布。一般情况下，机械零部件现有的、不同状态的样本故障数据都符合正态分布。
在不同运行状态下，因刚度非线性摩擦力等的影响，滚动轴承系统表现出不同的非线性特性，分形维数是用来定量刻画混沌吸引子“奇异”程度的一个重要参数，也可以用来刻画滚动轴承的故障状态。分形的维数有许多种类，最具代表性的是关联维数。关联维数能够体现未知系统的固有特性，并且机制不同的各类故障通常也具有不同的关联维数，可以作为系统故障特征量检测并区分滚动轴承的故障状态与故障模式。相同工作状态下的信号具有相近的关联维数，不同故障模式下的关联维数具有不同的数值，存在明显的可分性。而且，计算关联维数不需要建立系统方程，仅依靠一段样本数据就可以计算出该样本的关联维数，具备很强的通用性。实际应用中，样本数据往往比较匮乏，还需要采用自助法对样本进行重采样，将小样本问题转换为大样本问题来模仿未知分布。
自助法的基本原理是：在总体中抽取m个样本构成初始样本，之后随机、等概、独立、有放回地抽取m个样本单元，构成一个新的点集，即一个自助样本。数学描述为：设随机样本x＝(x1，x2…xn)是来自于某未知的总体分布f(x)，θ＝θ(f(x))为总体分布f的某个未知参数，fn(x)为抽样分布函数，为θ的估计。记估计误差为：
记为从fn(x)中抽样获得的再生样本，是由x*所获得的抽样分布函数。记：
rn*=θ^(fn*(x))-θ^(fn(x))---(5)]]>
称为tn的bootstrap统计量。在给定抽样分布函数fn(x)的条件下，取所有统计量的均值去模仿估计误差tn，则总体分布f(x)的参数
关联维数的分布通常符合正态分布，故而总体分布的参数有两个，分别为样本均值与标准差。根据自助训练得到的参数θ可以进一步得到样本均值与标准差。
目前计算关联维数的方法最主要的就是g-p算法，g-p算法是grassberger和procassi于1983年提出的一种比较容易实现的从实验数据中估算关联维数的算法。g-p算法如下：
c(m,n,r,τ)=2nm(nm-1)σi=1nmσj=1nmh(r-d(i,j))---(6)]]>
r为相空间超球半径，h(r-d(i，j))根据heaviside函数h(x)计算：
h(x)=1x&greaterequal;00x<0---(7)]]>
m维欧氏空间rm中的欧式距离为：
||yi-yj||=[σk=1m-1(xi+kτ-xj+kτ)2]12---(8)]]>
d(m,τ)=limr→0lncm(r)lnr---(9)]]>
画出标度曲线lnc-lnr图，取标度线中的线性区域的斜率来近似代替这个关联维数：
d(m,τ)=lncm(r)lnr---(10)]]>
式(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指数矩阵m，判断最大lyapunov指数矩阵m中数据是否全部大于零，若m中数据全部大于零则判定无故障信号存在，若m中矩阵数据不是全部大于零则判定有故障信号，进入基于关联维数的故障分类计算。
lyapunov指数的定义，是基于自治系统，而不是非自治系统的。因此，我们在求lyapunov指数的时候，都需要将原来的非自治系统转化为自治系统，再进行求解。holmes型简化的duffing方程(3)为二维非自治系统，通过代换z＝t，可以写成相应的三维自治系统：
x·=yy·=-0.5y+x3-x5+fcoszz·=1---(11)]]>
针对混沌系统，存在定理1：非自治系统的n个lyapunov指数等于其相应自治系统的前n个lyapunov指数，且该自治系统的第n+1个lyapunov指数等于零。
根据定理1，非自治系统(3)的两个lyapunov指数等于自治系统(11)的前两个lyapunov指数，因此现在只需考虑自治系统的lyapunov指数求解问题。
对于自治系统(11)的线性变分方程为
y·(t)=j(t)y(t),y(0)=i3---(12)]]>
y∈r3×3，i3是3×3的单位矩阵；j(t)为三维自治系统(11)的雅可比(jacobi)矩阵，可以看出，jacobi矩阵是奇异的，因为它的第三行元素都为零。这意味着y(t)的第三行元素都为常量。再利用初始条件y(0)＝i3，知y(t)的第三行元素分别为y31＝0，y32＝0以及y33＝1。由于y是非奇异的，将它进行qr分解，记为y(t)＝q(t)r(t)，r(t)为上三角矩阵，其对角元素均为正数。得到表达式(13)：
y(t)=y11y12y13y21y22y23001=q11q120q21q220001r11r12r130r22r23001---(13)]]>
根据前面lyapunov指数的定义，若δi为系统雅可比矩阵的第i个特征值，则lyapunov指数因为r33≡1，可知三维自治系统(11)的一个lyapunov指数始终为零，这也反映了定理1的内容。那么另外lyapunov指数的计算只需求取r11(t)和r22(t)。这样就仅考虑三维自治系统(11)的二维自治子系统就可以了。二维自治子系统的变分方程可以表示为：
y~·(t)=j~(t)y~(t),y~(0)=i2---(14)]]>
同理，i2为2×2的单位矩阵。为雅可比矩阵，其形式为：
j~(t)=01-&partiald;(ky+x3-x5+fcosz)&partiald;x-&partiald;(ky+x3-x5+fcosz)&partiald;y---(15)]]>
将式(15)简记为
j~(t)=j11j12j21j22---(16)]]>
将进行qr分解，即为代入式(14)有
q~·r~(t)+q~r~·=j~q~r~,q~(0)r~(0)=i2---(17)]]>
再将式(17)左乘右乘得到：
q~tq~·+r~·r~-1=q~tj~q~,q~(0)=i2,r~(0)=i2---(18)]]>
由于矩阵是斜对称矩阵，并且由于矩阵是上三角矩阵，所以也是上三角矩阵，因此，矩阵k(t)中的元素kij(t)为：
kij(t)=(q~t(t)j~(t)q~(t))ij,i>j0i=j-(q~t(t)j~(t)q~(t))jii<j---(19)]]>
从中可以得到
q~·(t)=q~(t)k(t),q~(0)=i2---(20)]]>
由式(19)和式(20)就可以确定了关于矩阵的微分方程。由于是斜对称的，由式(18)得到
r~·i′(t)=r·ii(t)rii(t)=(q~t(t)j~(t)q~(t))ij,ri′=0,i=1,2---(21)]]>
σ1=limt→∞r1′(t)t=limt→∞σ1(t)σ2=limt→∞r2′(t)t=limt→∞σ2(t)---(22)]]>
σ1和σ2就是自治系统(11)的另两个lyapunov指数，也即duffing方程式(3)的两个lyapunov指数。
将如式(3)中的混沌振子检测模型中周期策动力幅值f与策动力角频率ω每代入频率-阈值矩阵中的两个对应值，根据上面计算lyapunov指数的方法得到两个lyapunov指数，选择最大的lyapunov指数输出为最大lyapunov指数矩阵m中一个值，最终得到完整的最大lyapunov指数矩阵m。
步骤五、进行故障分类。若待检测信号含有故障信号，计算其关联维数。对照步骤一建立的不同故障类型的关联维数区间，进行故障分类，确定故障模式。关联维数的计算同步骤一中计算关联维数的方法相同。
实施例：
本实例采取华盛顿天主教大学轴承数据中心提供的滚动轴承故障信号进行验证。分别使用正常、内环故障、外环故障和滚动体故障四种状态下的样本信号对本发明针对机械零部件早期单点故障的lyapunov指数和关联维数相结合的故障检测和分类的方法进行检测验证，具体步骤如下：
步骤一、建立不同故障类型的关联维数区间。
根据已有的滚动轴承实际故障数据进行关联维数的计算，如表1所示。
表1滚动轴承不同状态振动信号的关联维数
对机械零部件现有的、不同状态的样本故障信号，计算相应的关联维数，如图3所示；根据统计学原理，对样本进行正态性检验；再应用基于小样本自采样方法进行自助训练，例如内环故障信号样本的均值和标准差分别如图4与图5所示；然后进行正态分布的参数估计，得到；根据所得到样本均值与标准差的估计值，计算出样本置信度为95％的分布区间。如图6所示为数据经过自助训练后分别计算出的检验区间图，分别为滚动体故障[2.3968，2.7567]，外环故障[3.0674，3.3534]，内环故障[3.4643，3.7571]。
步骤二、构造duffing混沌振子的频率矩阵
根据滚动轴承的尺寸及转速可以直接求出滚动轴承的各种故障特征频率。表2为华盛顿天主教大学轴承数据中心提供的用于计算滚动轴承的单点故障特征频率的实验值。
表2用于计算滚动轴承单点故障特征频率的实验值表
用采集信号时的轴的旋转角速度，乘以表2中的数值，即可计算出滚动轴承各种故障状态下对应的故障特征频率，进一步建立包含所有状态的频率矩阵p。
x′=ωyy·=ω(-0.5y+x3-x5+fcosωt+acosωt+σn(t))---(23)]]>
将策动力角频率ω分别设为频率矩阵p中的数值，不加入外部信号，调节f，使如式(3)中的混沌振子检测模型系统处于临界的混沌态，将lyapunov指数曲线的最后一个过零点所对应的策动力幅值f作为系统的临界阈值fd，并建立与频率相对应的频率-阈值矩阵。
步骤四、进行故障检测。
将如式(3)中的混沌振子检测模型中策动力幅值与角频率分别设为频率-阈值矩阵中的各个值；然后加入待检测信号计算最大lyapunov指数，输出最大lyapunov指数矩阵m；之后进行检验：若m中数据全部大于零则判定无故障信号存在，若m中矩阵数据不是全部大于零则判定有故障信号，进入基于关联维数的故障分类计算。
步骤五、进行故障分类。
若待检测信号含有故障信号，计算其关联维数，对照步骤一建立的不同故障类型的关联维数区间，如图6所示，关联维数落在哪一个区间，即为哪一种故障类型。
采取华盛顿天主教大学轴承数据中心提供的滚动轴承故障信号进行验证，分别采取正常、内环故障、外环故障和滚动体故障四种状态下20组样本信号进行检测，检测结果如表3：
表3检测结果
通过表3可以看出正常轴承信号和滚动体故障的判断全部正确，合格率为100％，内环故障有一个误判为外环故障，外环故障同样有一个误判为内环故障，这可能是两种故障的检验区间比较接近的缘故，其检测合格率为95％，可见本发明的方法能够实现机械零部件早期单点故障检测与分类，故障检测成功率高，具有明显的实际应用价值。