CN114486252B - 一种矢量模极大值包络的滚动轴承故障诊断方法 - Google Patents

一种矢量模极大值包络的滚动轴承故障诊断方法 Download PDF

Info

Publication number
CN114486252B
CN114486252B CN202210108671.4A CN202210108671A CN114486252B CN 114486252 B CN114486252 B CN 114486252B CN 202210108671 A CN202210108671 A CN 202210108671A CN 114486252 B CN114486252 B CN 114486252B
Authority
CN
China
Prior art keywords
rolling bearing
maximum
coordinate system
vibration acceleration
vector
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210108671.4A
Other languages
English (en)
Other versions
CN114486252A (zh
Inventor
王俨剀
史浩宁
陈康
王彤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202210108671.4A priority Critical patent/CN114486252B/zh
Publication of CN114486252A publication Critical patent/CN114486252A/zh
Application granted granted Critical
Publication of CN114486252B publication Critical patent/CN114486252B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

一种矢量模极大值包络的滚动轴承故障信息融合方法，在滚动轴承圆周表面布置两个加速度传感器，用于采集滚动轴承的振动信号。对所采集信号进行正交化处理，去除与故障无关的常矢量，留下与故障有关的旋转矢量；取模，得到模长序列。通过得到的连续函数曲线中的按时间顺序从小到大排列的第一个极大值点、第二个极大值点，以及倒数第一个极大值点、倒数第二个极大值点的坐标的得出模长序列开端的极大值点和末端的极大值点。通过三次样条插值方法得到三次样条插值函数。再在此基础上，通过傅里叶变换，得到矢量模极大值包络谱。本发明能够得到更清晰的调制边频带，可以更明显地、更不失真地提取出故障特征频率，有利于滚动轴承局部故障诊断。

Description

一种矢量模极大值包络的滚动轴承故障诊断方法
技术领域
本发明涉及航空发动机状态检测及故障诊断领域，是一种基于多通道信息融合的矢量模极大值包络滚动轴承故障诊断分析方法。
背景技术
对于航空发动机滚动轴承乃至民用旋转机械的滚动轴承的状态检测及故障诊断，普遍的方法是，通过传感器测得时域信号，在对其进行分析，从而诊断出故障与异常。常用的方法可以从“时域分析”、“频域分析”、“时频域分析”三个角度进行归类。
“时域分析”是在时域内通过各种指标来判断是否出现故障。
“频域分析”是将时域信号通过傅里叶变换转换成频域信号，再在频域内分析频率组成成分和信号幅值变化，来判断出故障的类型。
“时频域分析”是通过时间和频率的二维联合分布来描述非平稳信号的统计特征随时间的变化。传统的时域和频域分析方法都是基于平稳理论而提出的，故适用于分析解决外在表现为平稳振动信号的故障。但滚动轴承更多的故障所产生的动态响应是非平稳过程，其振动信号是非平稳信号，此时，单纯地通过时域分析法或频域分析法进行分析并不适用。
一般地，对于一个滚动轴承，通常使用一个传感器对其振动信号进行收集。
但经仿真与实验发现，在同一时刻下，滚动轴承外环、内环和滚动体表面损伤产生的冲击响应在两个正交方向上的表现是不同的。这说明，仅仅依赖单一方向上的振动数据进行分析，提取到的滚动轴承故障振动特征可能非常微弱，甚至提取不到故障特征。
而且，径向静载荷分布、局部损伤点与传感器测点相对位置变化、以及转子不平衡引起的动载荷均会对滚动轴承局部损伤引起的冲击振动幅值产生影响，反映在包络谱中，就是在故障特征频率及其谐频周围出现边频。由于不同因素的影响强弱程度不同、引起的调制频率大小、相位不同，相互叠加的结果就是在故障特征频率及其谐频周围出现极为复杂的边频带，甚至使主频发生偏移。
上述原因导致利用故障特征频率以及边频特征无法准确判断滚动轴承是否发生故障，容易造成漏诊、误诊，引起故障诊断系统的虚警率和误警率上升。
造成以上后果的主要原因是能够突出滚动轴承局部故障特征的振动是二维的、在径向平面内旋转的矢量。利用单个传感器去采集这一平面振动矢量势必会引入虚假信息并且难以全面描述其振动特征。
在申请号为CN201611110937.X的发明创造中，提出了一种基于互补随机共振滤波器的微弱信号增强检测方法。通过采用更先进的双通道互补随机共振滤波器，经过加权谱峭度指标自适应地调整合适的系统参数，利用互补通道相位差异噪声去增强主通道的微弱周期信号。得到了实现滚动轴承微弱故障特征频率增强，提高滚动轴承故障诊断的准确度的效果。但该方法所使用的采集数据是单通道传感器得到的，并不是多通道的。只是在基于单通道的测试诊断基础上进行改良，而非跨越到双通道的层面。
在公开号为CN112766331A的发明创造中，提出了一种基于多通道CNN多元信息融合的滚动轴承故障诊断方法。通过将采集到的滚动轴承双测数据打标签处理后，输入进初始化的多信息融合的卷积神经网络模型中，得出满足数据集迭代训练次数后的网格结构，并将其应用于模型评估。得到了使得网络模型不仅利用到了多通道数据而且还没有增加网络模型的计算过程并且达到了比单通道网络模型更高且更稳定的准确率的结果。但该方法中并没有涉及对双侧数据的处理，而且其采用地打标签及需要一定次数训练的方式不适于对发动机状态进行实时监测。
发明内容
为了克服目前航空发动机常用的滚动轴承故障诊断方法中，存在的单通道测量数据的片面性以及多通道测量数据处理与使用过程复杂，运算量大的不足，本发明提出一种矢量模极大值包络的滚动轴承故障信息融合方法。
本发明的具体过程是：
步骤1：构造10000个复数信号矢量
Figure BDA0003494664630000021
在滚动轴承周向所在的同一平面内，存在一新坐标系，所述新坐标系的坐标轴为
Figure BDA0003494664630000022
Figure BDA0003494664630000023
所述新坐标系的坐标轴/>
Figure BDA0003494664630000024
的方向与原坐标系中y轴的方向一致，记为新坐标系中的竖直方向。所述新坐标系的坐标轴/>
Figure BDA0003494664630000025
的方向与新坐标系中/>
Figure BDA0003494664630000026
的方向垂直，记为新坐标系中的水平方向。
通过测量得到滚动轴承的水平方向的振动加速度信号x(t)与竖直方向的振动加速度信号y(t)，并通过所述x(t)与y(t)构造复信号
Figure BDA0003494664630000027
在滚动轴承圆周表面布置两个相同型号的加速度传感器，并使这两个传感器之间的夹角为25°～155°。以采样频率为10kHz的等时间采样方法，在滚动轴承外环静止且内环转速为2400rpm时，得到在3s时段内，所述两个加速度传感器分别测得的振动加速度信号
Figure BDA0003494664630000031
与/>
Figure BDA0003494664630000032
单位为g，且g＝9.8m/s2。其中/>
Figure BDA0003494664630000033
Figure BDA0003494664630000034
中各时刻分别与0s,0.0001s,0.0002s…2.9998s,2.9999s,3.0000s时刻一一对应。
选取3s时间段中的1.0001s时刻作为用于处理数据时的0.0001s时刻；选取3s时间段中的2s时刻作为用于处理数据时的1s时刻，以此排除接通传感器与关闭传感器时对数据采集的影响。故用于分析处理的采样数据包含的时长为0.9999s。
在所述0.9999s的时间段内，各加速度传感器能够测得在0.0001s,0.0002s,0.0003s…0.9998s,0.9999s,1.0000s共1000个时刻的振动加速度数据共计10000对。分别是：x(0.0001)，x(0.0002)，x(0.0003)…x(0.9998)，x(0.9999)，x(1.0000)，记为x(ti)；以及y(0.0001)，y(0.0002)，y(0.0003)…y(0.9998)，y(0.9999)，y(1.0000)，记为y(ti)。ti中i的取值为1，2，3…9998，9999，10000，所述ti对应的时刻分别是0.0001s,0.0002s,0.0003s…0.9998s,0.9999s,1.0000s。
该新坐标系对应的复平面内的复数坐标系为将新坐标系的纵轴
Figure BDA0003494664630000035
乘以虚数单位j后，其余不变，所得的复数坐标系。
若布置的两个相同型号的加速度传感器的夹角为90°时，则将x(ti)直接记为
Figure BDA0003494664630000036
表示10000个新坐标系内的水平方向的振动加速度信号，并将y(ti)直接记为/>
Figure BDA0003494664630000037
表示10000个新坐标系内的竖直方向的振动加速度信号。
在实际的航空发动机的结构中，由于加工、装配的因素，两个相同型号的加速度传感器如果法正交布置，即二者之间的夹角θ≠90°。则需通过公式(1)与(2)进行正交变换：
Figure BDA0003494664630000038
Figure BDA0003494664630000039
其中，式(1)是将原坐标系中的沿y轴方向的10000个振动加速度信号y(ti)作为新坐标系中的沿
Figure BDA00034946646300000310
轴方向的10000个振动加速度信号/>
Figure BDA00034946646300000311
式(2)是将原坐标系中的沿x轴方向的10000个振动加速度信号x(ti)变换为新坐标系中的沿/>
Figure BDA00034946646300000312
轴方向的10000个振动加速度信号/>
Figure BDA00034946646300000313
将新坐标系中的两方向相互垂直的10000对振动加速度信号
Figure BDA00034946646300000314
与/>
Figure BDA00034946646300000315
通过
Figure BDA0003494664630000041
与虚数单位j相乘，再与/>
Figure BDA0003494664630000042
加和，构造成10000个复数坐标系内的振动加速度复数信号矢量/>
Figure BDA0003494664630000043
用以表征10000个变换到复数坐标系里的振动加速度复数信号矢量形式。
步骤2，确定与故障有关的旋转矢量
Figure BDA0003494664630000044
通过公式：
Figure BDA0003494664630000045
得到10000个/>
Figure BDA0003494664630000046
振动加速度信号的平均值/>
Figure BDA0003494664630000047
所述平均值为新坐标系中的沿/>
Figure BDA0003494664630000048
轴方向的振动加速度信号平均值。
通过公式
Figure BDA0003494664630000049
得到10000个/>
Figure BDA00034946646300000410
振动加速度信号的平均值/>
Figure BDA00034946646300000411
该平均值为新坐标系中的沿/>
Figure BDA00034946646300000412
轴方向的振动加速度信号平均值。
将
Figure BDA00034946646300000413
与虚数单位j相乘后与/>
Figure BDA00034946646300000414
相加，构造出复数坐标系内的振动加速度复数信号矢量/>
Figure BDA00034946646300000415
用以表示复数坐标系内的振动加速度信号平均值的矢量形式，代表故障无关的常矢量。在得到的10000个矢量/>
Figure BDA00034946646300000416
中分别有一个常矢量/>
Figure BDA00034946646300000417
从各矢量/>
Figure BDA00034946646300000418
中分别减去该常矢量，分别得到10000个复数坐标系内与故障有关的旋转矢量/>
Figure BDA00034946646300000419
Figure BDA00034946646300000420
步骤3，确定振动旋转矢量的模长序列：
对步骤2中得到的10000个与故障有关的旋转矢量
Figure BDA00034946646300000421
通过实部平方与虚部平方之和开二次方的运算，进行取模，得到的模为：
Figure BDA00034946646300000422
所述取的模为10000个，该10000个模分别对应ti各时刻。
各模长均为L(ti)，分别是L(t1),L(t2)……L(t9998),L(t9999),L(t10000)；各模长值分别代表相应时刻的振动加速度的幅值；该模长值单位为g，为标量。
将得到系列模长值按时间从小到大的顺序t1＝0.0001s,t2＝0.0002s……t9998＝0.9998s，t9999＝0.9999s，t10000＝1.0000s排列在加速度振幅--时间图上；所述加速度幅值--时间图的横坐标是时间t，单位为s；纵坐标是加速度幅值A，亦即模长值，该加速度幅值的单位是g。得到一系列横坐标间隔为0.0001s的离散的点。将相邻的点用直线连接，得到连续函数曲线。
步骤4，确定连续函数曲线中的两端的用于线性插值的极大值点坐标：
所述用于线性插值的极大值点是分别位于所述得到连续函数曲线两端的各两个极大值，具体分别是该连续函数曲线中的按时间顺序从小到大排列的第一个极大值点、第二个极大值点，以及该按时间顺序排列的倒数第一个极大值点、倒数第二个极大值点。
得到所述连续函数曲线中的按时间顺序从小到大排列的第一个极大值点的坐标P1、第二个极大值点的坐标P2，以及倒数第一个极大值点的坐标Ps-1、倒数第二个极大值点的坐标Ps。分别为：
P1：(0.0005，4870.84)，
P2：(0.0009，3638.49)。
Ps-1：(0.9991，127.579)，
Ps：(0.9996，103.047)。
根据数学上对极大值点的定义得到s个极大值点：P1，P2，P3…Ps-2，Ps-1，Ps。所述s个极大值点的坐标分别为：P1(t1,A(t1))，P2(t2,A(t2)),P3(t3,A(t3))…Ps-2(ts-2,A(ts-2))，Ps-1(ts-1,A(ts-1))，Ps(ts,A(ts))。其中，0.0001s＜t1＜t2＜t3。。。＜ts-2＜ts-1＜ts＜1s。
步骤5，通过线性插值得出区间端点处的极大值点坐标：
在所述P1(0.0005，4870.84)与P2(0.0009，3638.49)之间构造直线，在加速度振幅--时间图上，A＝－3080875t+6411.2775。当t＝0.0001s时，代入A中得到A(0.0001)＝6103.19g。得到采样区间t＝0.0001s时刻的原图像中的极大值点，该极大值点的坐标为(0.0001，6103.19)。
在所述Ps-1(0.9991，127.579)与Ps(0.9996，103.047)之间构造直线，在加速度振幅--时间图上，为：A＝－49064t+49147.4214。当t＝1s时，代入A中得到A(1)＝83.4214g。得到采样区间t＝1.0000s时刻的原图像中的极大值点，该极大值点的坐标为(1，83.4214)。
步骤6，确定采样区间内的三次样条插值函数：
利用三次样条插值法，对采样区间内的各极大值点的坐标以及步骤5中插值得到的区间两端端点处的极大值点的坐标进行插值，得到三次样条插值函数。具有过程是：
Ⅰ采样区间内的极大值点与区间端点处的极大值点均为三次样条插值所用的数据节点。采样区间内的极大值点个数为s，通过插值得到的区间端点处的极大值点个数为2，故所述数据节点的个数为s+2，按横坐标时间顺序从小到大记为(t0，A0)，(t1，A1)，(t2，A2)......(ts+1，As+1)。各所述数据节点在横坐标方向上，将采样区间[0.0001s，1.0000s]分成了s+1个子区间：(t0，t1)，(t1，t2)，(t2，t3)…(ts，ts+1)。且三次样条插值函数的边界条件为自然边界条件，即区间两端端点的二阶导数均为0。
II计算步长hi＝ti+1-ti。其中，i＝0，1，2…s。
III将s+2个数据节点的坐标值，与自然边界条件：“m0＝0及ms+1＝0”入如下矩阵中：
Figure BDA0003494664630000061
并解出得到s+2个数据节点处的二阶导数的值：m0，m1，m2…ms+1。
IV计算样条曲线的系数：
ai＝Ai是各子区间中三次样条插值函数的常数项。
Figure BDA0003494664630000062
是各子区间[ti，ti+1]中三次样条插值函数的一次方项系数。i＝0，1，2...s。
Figure BDA0003494664630000063
是各子区间[ti，ti+1]s中，三次样条插值函数的平方项系数；i＝0，1，2…s。
Figure BDA0003494664630000064
是各子区间[ti，ti+1]中三次样条插值函数的立方项系数。i＝0，1，2…s。
V在各子区间[ti，ti+1]内构造曲线Ai(t)＝ai+bi(t-ti)+ci(t-ti)2+di(t-ti)3；则各子区间[ti，ti+1]，i＝0，1，2…s内的共计s+1条曲线：A0(t)，A1(t)…As-1(t)，As(t)在采样区间内s个极大值点的坐标处是二阶连续的。将所述s+1条在各子区间[ti，ti+1]内构造曲线记为函数A(t)，其中，t的取值范围：[t0，ts+1]。所述函数A(t)是通过三次样条插值方法得到的三次样条插值函数A(t)。
步骤7，确定矢量模极大值包络谱：
对得到的三次样条插值函数A(t)进行傅里叶变换，得到矢量模极大值包络谱。具体是：
对区间[t0，ts+1]内，进行傅里叶变换：
Figure BDA0003494664630000071
其中，A表示按所述数据处理方法得到的滚动轴承的加速度幅值，单位是g。f表示按所述数据处理方法得到的滚动轴承的频率，单位是Hz。对步骤7得到的三次样条插值函数A(t)进行傅里叶变换后，将振动加速度振幅--时间曲线变换为振动加速度振幅--频率曲线。F[A(t)]表示对A(t)进行傅里叶变换。
得到矢量模极大值包络谱。
步骤8，结合矢量模极大值包络谱对滚动轴承故障进行诊断与分析。
根据得到的所述矢量模极大值包络谱，在滚动轴承外环故障特征频率131Hz及滚动轴承外环故障特征频率两侧的边频116Hz与146Hz处具有较大加速度值，在滚动轴承外环故障特征频率的二倍频261Hz及滚动轴承外环故障特征频率两侧的二倍频的边频246Hz与276Hz处，以及滚动轴承外环故障特征频率的三倍频392Hz及滚动轴承外环故障特征频率的三倍频两侧的边频377Hz与407Hz处，均有较大的加速度值。判定该滚动轴承外环处存在损伤故障。
本发明提出了一种基于多通道信息融合的滚动轴承故障诊断方法。具体流程是：在滚动轴承圆周表面，布置两个相同型号的加速度传感器，并使这两个传感器之间的夹角在25°～155°之间，采集滚动轴承的振动信号，再对所采集信号进行正交化处理，去除与故障无关的常矢量，留下与故障有关的旋转矢量，再取模，得到模长序列，通过得到的连续函数曲线中的按时间顺序从小到大排列的第一个极大值点、第二个极大值点，以及，倒数第一个极大值点、倒数第二个极大值点的坐标的得出模长序列开端的极大值点和末端的极大值点，再通过三次样条插值方法得到三次样条插值函数。再在此基础上，通过傅里叶变换，得到矢量模极大值包络谱。
与现有技术相比较，本发明的有益效果在于：
通过融合两个传感器采集到的振动数据，得到能够全面表达滚动轴承故障特征信息的二元复信号，避免了单一通道测得的数据信号对与反应故障信息的片面性。通过对其进行的预处理，消除了振动矢量中不随时间变化的常矢，使得包含故障特征的信息更加明显，如图3a，3b所示。并且，提出的“矢量模极大值包络”算法较比传统的包络算法，能够得到更清晰的调制边频带，可以更明显地、更不失真地提取出故障特征频率，有利于滚动轴承局部故障诊断。更具体地，对于滚动轴承每一部件的故障的诊断增益性体现在如下3个部件的实例中。
对于外环故障特征频率为255.6Hz的滚动轴承的外环局部故障，当转速为2400rpm时，单独利用水平方向信号分析得到的包络谱如图5所示。单独利用竖直方向信号分析得到的包络谱如图6所示。利用本申请所提到的一种矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱如图7所示。可见，虽然利用单通道信号分析得到的包络谱与利用矢量模极大值包络算法得到的包络谱在频率成分分辨上没有明显差异，均可以看出外环故障特征频率255.6Hz及其左右侧的边频214.9Hz与294.6Hz，但是利用矢量模极大值包络算法得到的包络谱上的特征频率255.6Hz及其左右侧的边频214.9Hz与294.6Hz处具有更大的幅值，使得外环故障特征频率的提取更加明显。这说明，矢量模极大值包络方法能够有效防止因传感器主响应方向与振动响应方向不一致造成幅值削减的现象，从而增强外环故障特征信息，有效诊断出滚动轴承外环局部故障。对应的增益效果如下表格所示：
Figure BDA0003494664630000081
对于内环故障特征频率为171.4Hz的滚动轴承的内环局部故障，当转速为1200rpm时，单独利用水平方向信号分析得到的包络谱如图8所示。单独利用竖直方向信号分析得到的包络谱如图9所示。利用本申请所出的矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱如图10所示。可见：利用单通道信号分析得到的包络谱的内环故障特征频率171.4Hz处存在大量的调制频带，而利用矢量模极大值包络算法得到的包络谱调制的内环故障特征频率171.4Hz处周围的调制频带数量大大减少。并且，矢量模极大值包络谱中内环故障特征频率171.4Hz处对应的幅值更为显著。这说明，矢量模极大值包络方法能够有效防止因单一传感器与内环损伤点相对位置变化引入的调制情况，使得主频能量更为集中。该方法可将径向静载荷引起的调制与传感器与内环损伤点相对位置变化引入的调制分离开，增强内环故障特征信息，能够有效诊断出滚动轴承内环局部故障。对应的增益效果如下表格所示：
Figure BDA0003494664630000091
对于滚动体故障特征频率为156.9Hz的滚动轴承的滚动体局部故障，当转速为3000rpm时，单独利用水平方向信号分析得到的包络谱如图11所示。单独利用竖直方向信号分析得到的包络谱如图12所示。利用本申请所提到的一种矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱如图13所示。可见：利用单通道信号分析得到的包络谱在滚动体故障特征频率156.9Hz周围存在大量的调制频带，而利用矢量模极大值包络算法得到的包络谱在滚动体故障特征频率156.9Hz周围调制频带数量减少，并且，矢量模极大值包络谱中滚动体故障特征频率156.9Hz对应的的幅值更为显著。这说明，矢量模极大值包络方法能够有效防止因单一传感器与滚动体损伤点相对位置变化引入的调制情况，使得主频能量更为集中。该方法可将径向静载荷引起的调制与传感器与滚动体损伤点相对位置变化引入的调制分离开，增强滚动体故障特征信息，能够有效诊断出滚动轴承滚动体局部故障。对应的增益效果如下表格所示：
Figure BDA0003494664630000092
Figure BDA0003494664630000101
附图说明
图1为本发明技术方案示意图。
图2a为两传感器正交放置示意图；图2b为两传感器非正交放置示意图。
图3a，3b为假设构造的包含常矢与旋转矢的振动矢量图。
图4为线性插值得模长序列开端的极大值点和末端的极小值点示意图及插值后通过三次样条插值的方法获得的振动矢量模长序列的上包络信号图。
图5为对于存在外环故障的滚动轴承，单独利用水平方向信号分析得到的包络谱
图6为对于存在外环故障的滚动轴承，单独利用竖直方向信号分析得到的包络谱
图7为对于存在外环故障的滚动轴承，利用本申请所提到的一种矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱
图8为对于存在内环故障的滚动轴承，单独利用水平方向信号分析得到的包络谱
图9为对于存在内环故障的滚动轴承，单独利用竖直方向信号分析得到的包络谱
图10为7为对于存在内环故障的滚动轴承，利用本申请所提到的一种矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱
图11为对于存在滚动体故障的滚动轴承，单独利用水平方向信号分析得到的包络谱
图12为对于存在滚动体故障的滚动轴承，单独利用竖直方向信号分析得到的包络谱
图13为对于存在滚动体故障的滚动轴承，利用本申请所提到的一种矢量模极大值包络的滚动轴承故障信息融合方法得到的包络谱
图14为振动旋转矢量的模长序列实例图。
图15为首端线性插值及首端线性插值后对首端进行包络图。
图16为尾端线性插值及尾端线性插值后对尾端进行包络图。
图17为利用三次样条插值法对采样区间内的各极大值点以及步骤六中插值得到的左右端点处的极大值点进行插值，得到三次样条插值函数的图。
图18为矢量模极大值包络谱图。
具体实施方式
本实施例是一种矢量模极大值包络的滚动轴承故障诊断方法，具体过程是：
步骤1：构造10000个复数信号矢量
Figure BDA0003494664630000111
按常规方法测量得到滚动轴承的水平方向的振动加速度信号x(t)与竖直方向的振动加速度信号y(t)，并通过所述x(t)与y(t)构造复信号
Figure BDA0003494664630000112
在滚动轴承圆周表面布置两个相同型号的加速度传感器，并使这两个传感器之间的夹角为25°～155°。以采样频率为10kHz的等时间采样方法，在滚动轴承外环静止且内环转速为2400rpm时，得到在3s时段内，所述两个加速度传感器分别测得的振动加速度信号
Figure BDA0003494664630000113
与/>
Figure BDA0003494664630000114
单位为g，且g＝9.8m/s2。其中/>
Figure BDA0003494664630000115
Figure BDA0003494664630000116
中各时刻分别与0s,0.0001s,0.0002s…2.9998s,2.9999s,3.0000s一一对应。
选取3s时间段中的1.0001s时刻作为用于处理数据时的0.0001s时刻；选取3s时间段中的2s时刻作为用于处理数据时的1s时刻，以此排除接通传感器与关闭传感器时对数据采集的影响。故用于分析处理的采样数据包含的时长为0.9999s。
根据10kHz的采样频率，在所述0.9999s的时间段内，各加速度传感器能够测得在0.0001s,0.0002s,0.0003s…0.9998s,0.9999s,1.0000s共1000个时刻的振动加速度数据共计10000对。分别是：x(0.0001)，x(0.0002)，x(0.0003)…x(0.9998)，x(0.9999)，x(1.0000)，记为x(ti)；以及y(0.0001)，y(0.0002)，y(0.0003)…y(0.9998)，y(0.9999)，y(1.0000)，记为y(ti)。ti中i的取值为1，2，3…9998，9999，10000，所述ti对应的时刻分别为0.0001s,0.0002s,0.0003s…0.9998s,0.9999s,1.0000s时刻。
记滚动轴承周向所在平面内，原坐标系的坐标轴为x与y，方向分别是滚动轴承中心沿两个加速度传感器的方向。
记滚动轴承周向所在的同一平面内，存在一新坐标系，所述新坐标系的坐标轴为
Figure BDA0003494664630000117
Figure BDA0003494664630000118
所述新坐标系的坐标轴/>
Figure BDA0003494664630000119
的方向与原坐标系中y轴的方向一致，记为新坐标系中的竖直方向。所述新坐标系的坐标轴/>
Figure BDA00034946646300001110
的方向与新坐标系中/>
Figure BDA00034946646300001111
的方向垂直，记为新坐标系中的水平方向。
该新坐标系对应的复平面内的复数坐标系为将新坐标系的纵轴
Figure BDA00034946646300001112
乘以虚数单位j后，其余不变，所得的复数坐标系。
若布置的两个相同型号的加速度传感器的夹角为90°时，则将x(ti)直接记为
Figure BDA00034946646300001113
表示10000个新坐标系内的水平方向的振动加速度信号，并将y(ti)直接记为/>
Figure BDA00034946646300001114
表示10000个新坐标系内的竖直方向的振动加速度信号。
在实际的航空发动机的结构中，由于加工、装配的因素，两个相同型号的加速度传感器如果法正交布置，即二者之间的夹角θ≠90°。则需通过公式(1)与(2)进行正交变换：
Figure BDA0003494664630000121
Figure BDA0003494664630000122
其中，式(1)是将原坐标系中的沿y轴方向的10000个振动加速度信号y(ti)作为新坐标系中的沿
Figure BDA0003494664630000123
轴方向的10000个振动加速度信号/>
Figure BDA0003494664630000124
式(2)是将原坐标系中的沿x轴方向的10000个振动加速度信号x(ti)变换为新坐标系中的沿/>
Figure BDA0003494664630000125
轴方向的10000个振动加速度信号/>
Figure BDA0003494664630000126
将新坐标系中的两方向相互垂直的10000对振动加速度信号
Figure BDA0003494664630000127
与/>
Figure BDA0003494664630000128
通过
Figure BDA0003494664630000129
与虚数单位j相乘，再与/>
Figure BDA00034946646300001210
加和，构造成10000个复数坐标系内的振动加速度复数信号矢量/>
Figure BDA00034946646300001211
用以表征10000个变换到复数坐标系里的振动加速度复数信号矢量形式。
步骤2，确定与故障有关的旋转矢量
Figure BDA00034946646300001212
通过公式：
Figure BDA00034946646300001213
得到10000个/>
Figure BDA00034946646300001214
振动加速度信号的平均值/>
Figure BDA00034946646300001215
所述平均值为新坐标系中的沿/>
Figure BDA00034946646300001228
轴方向的振动加速度信号平均值。
通过公式
Figure BDA00034946646300001216
得到10000个/>
Figure BDA00034946646300001217
振动加速度信号的平均值/>
Figure BDA00034946646300001218
该平均值为新坐标系中的沿/>
Figure BDA00034946646300001219
轴方向的振动加速度信号平均值。
将
Figure BDA00034946646300001220
与虚数单位j相乘后与/>
Figure BDA00034946646300001221
相加，构造出复数坐标系内的振动加速度复数信号矢量/>
Figure BDA00034946646300001222
用以表示复数坐标系内的振动加速度信号平均值的矢量形式，代表故障无关的常矢量。在得到的10000个矢量/>
Figure BDA00034946646300001223
中分别有一个常矢量/>
Figure BDA00034946646300001224
从各矢量/>
Figure BDA00034946646300001225
中分别减去该常矢量，分别得到10000个复数坐标系内与故障有关的旋转矢量/>
Figure BDA00034946646300001226
Figure BDA00034946646300001227
步骤3，确定振动旋转矢量的模长序列：
对步骤2中得到的10000个与故障有关的旋转矢量
Figure BDA0003494664630000131
通过实部平方与虚部平方之和开二次方的运算，进行取模，得到的模为：
Figure BDA0003494664630000132
所述取的模为10000个，该10000个模分别对应ti各时刻。
各模长均为L(ti)，分别是L(t1),L(t2)……L(t9998),L(t9999),L(t10000)；各模长值分别代表相应时刻的振动加速度的幅值；该模长值单位为g，为标量。
将得到系列模长值按时间顺序排列在加速度振幅--时间图上：
t1＝0.0001s,t2＝0.0002s……t9998＝0.9998s，t9999＝0.9999s，t10000＝1.0000s
所述加速度幅值--时间图的横坐标是时间t，单位为s；纵坐标是加速度幅值A，亦即模长值，该加速度幅值的单位是g。得到一系列横坐标间隔为0.0001s的离散的点。将相邻的点用直线连接，得到连续函数曲线，如图14所示。
步骤4，确定连续函数曲线中的两端的用于线性插值的极大值点坐标：
所述用于线性插值的极大值点是分别位于所述得到连续函数曲线两端的各两个极大值，具体分别是该连续函数曲线中的按时间顺序从小到大排列的第一个极大值点、第二个极大值点，以及该按时间顺序排列的倒数第一个极大值点、倒数第二个极大值点。
得到所述连续函数曲线中的按时间顺序从小到大排列的第一个极大值点的坐标P1、第二个极大值点的坐标P2，以及倒数第一个极大值点的坐标Ps-1、倒数第二个极大值点的坐标Ps。分别为：
P1：(0.0005，4870.84)，
P2：(0.0009，3638.49)。
Ps-1：(0.9991，127.579)，
Ps：(0.9996，103.047)。
根据数学上对极大值点的定义得到s个极大值点：P1，P2，P3…Ps-2，Ps-1，Ps。所述s个极大值点的坐标分别为：P1(t1,A(t1))，P2(t2,A(t2)),P3(t3,A(t3))…Ps-2(ts-2,A(ts-2))，Ps-1(ts-1,A(ts-1))，Ps(ts,A(ts))。其中，0.0001s＜t1＜t2＜t3。。。＜ts-2＜ts-1＜ts＜1s。
步骤5，通过线性插值得出区间端点处的极大值点坐标：
在所述P1(0.0005，4870.84)与P2(0.0009，3638.49)之间构造直线，在加速度振幅--时间图上，A＝-3080875t+6411.2775。当t＝0.0001s时，代入A中得到A(0.0001)＝6103.19g。得到采样区间t＝0.0001s时刻的原图像中的极大值点，该极大值点的坐标为(0.0001，6103.19)。如图15所示。
在所述Ps-1(0.9991，127.579)与Ps(0.9996，103.047)之间构造直线，在加速度振幅--时间图上，为：A＝-49064t+49147.4214。当t＝1s时，代入A中得到A(1)＝83.4214g。得到采样区间t＝1.0000s时刻的原图像中的极大值点，该极大值点的坐标为(1，83.4214)。如图16所示。
步骤6，确定采样区间内的三次样条插值函数：
利用三次样条插值法，对采样区间内的各极大值点的坐标以及步骤5中插值得到的区间两端端点处的极大值点的坐标进行插值，得到三次样条插值函数。具有过程是：
I采样区间内的极大值点与区间端点处的极大值点均为三次样条插值所用的数据节点。本实施例中，采样区间内的极大值点个数为s，步骤5中通过插值得到的区间端点处的极大值点个数为2，故所述数据节点的个数为s+2，按横坐标时间顺序从小到大记为(t0，A0)，(t1，A1)，(t2，A2)……(ts+1，As+1)。各所述数据节点在横坐标方向上，将采样区间[0.0001s，1.0000s]分成了s+1个子区间：(t0，t1)，(t1，t2)，(t2，t3)...(ts，ts+1)。且三次样条插值函数的边界条件为自然边界条件，即区间两端端点的二阶导数均为0。
II计算步长hi＝ti+1-ti。其中，i＝0，1，2…s。
III将s+2个数据节点的坐标值，与自然边界条件：“m0＝0及ms+1＝0”入如下矩阵中：
Figure BDA0003494664630000141
并解出得到s+2个数据节点处的二阶导数的值：m0，m1，m2…ms+1。
IV计算样条曲线的系数：
ai＝Ai是各子区间中三次样条插值函数的常数项。
Figure BDA0003494664630000151
是各子区间[ti，ti+1]中三次样条插值函数的一次方项系数。i＝0，1，2...s。
Figure BDA0003494664630000152
是各子区间[ti，ti+1]s中，三次样条插值函数的平方项系数；i＝0，1，2…s。
Figure BDA0003494664630000153
是各子区间[ti，ti+1]中三次样条插值函数的立方项系数。i＝0，1，2…s。
V在各子区间[ti，ti+1]内构造三次样条插值曲线Ai(t)＝ai+bi(t-ti)+ci(t-ti)2+di(t-ti)3；
则各子区间[ti，ti+1]，i＝0，1，2…s内的共计s+1条三次样条插值曲线：A0(t)，A1(t)…As-1(t)，As(t)在采样区间内s个极大值点的坐标处是二阶连续的。将所述s+1条在各子区间[ti，ti+1]内构造的三次样条插值曲线：A0(t)，A1(t)…As-1(t)，As(t)，记为抽象函数A(t)，其中自变量t的取值范围是：[t0，ts+1]。当抽象函数A(t)的自变量t值，落在各子区间[ti，ti+1]内时，抽象函数A(t)所对应的具体函数，与在该区间内所构造的上述三次样条插值曲线一致。所述函数A(t)是通过三次样条插值方法得到的区间[t0，ts+1]内的三次样条插值函数。如图17。
步骤7，确定矢量模极大值包络谱：
对步骤6得到的三次样条插值函数A(t)进行傅里叶变换，得到矢量模极大值包络谱。具体是：
对区间[t0，ts+1]内，进行傅里叶变换：
Figure BDA0003494664630000154
其中，A表示按所述数据处理方法得到的滚动轴承的加速度幅值，单位是g。f表示按所述数据处理方法得到的滚动轴承的频率，单位是Hz。对步骤7得到的三次样条插值函数A(t)进行傅里叶变换后，将振动加速度振幅--时间曲线变换为振动加速度振幅--频率曲线。F[A(t)]表示对A(t)进行傅里叶变换。
得到矢量模极大值包络谱。如图18。
步骤8，结合矢量模极大值包络谱对滚动轴承故障进行诊断与分析。
根据得到的所述矢量模极大值包络谱，在滚动轴承外环故障特征频率131Hz及滚动轴承外环故障特征频率两侧的边频116Hz与146Hz处具有较大加速度值，在滚动轴承外环故障特征频率的二倍频261Hz及滚动轴承外环故障特征频率两侧的二倍频的边频246Hz与276Hz处，以及滚动轴承外环故障特征频率的三倍频392Hz及滚动轴承外环故障特征频率的三倍频两侧的边频377Hz与407Hz处，均有较大的加速度值。判定该滚动轴承外环处存在损伤故障。

Claims (9)

1.一种矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，具体过程是：
步骤1：构造10000个复数信号矢量
Figure QLYQS_1
滚动轴承周向所在的同一平面内，存在一新坐标系，所述新坐标系的坐标轴为
Figure QLYQS_2
；所述新坐标系的坐标轴/>
Figure QLYQS_3
的方向与原坐标系中y轴的方向一致，记为新坐标系中的竖直方向；所述新坐标系的坐标轴/>
Figure QLYQS_4
的方向与新坐标系中/>
Figure QLYQS_5
的方向垂直，记为新坐标系中的水平方向；
通过测量得到滚动轴承的水平方向的振动加速度信号x(t)与竖直方向的振动加速度信号y(t)，并通过所述x(t)与y(t)构造复信号
Figure QLYQS_6
该新坐标系对应的复平面内的复数坐标系为将新坐标系的纵轴
Figure QLYQS_7
乘以虚数单位j后，其余不变，所得的复数坐标系；
步骤2，确定与故障有关的旋转矢量
Figure QLYQS_8
通过公式：
Figure QLYQS_11
得到10000个/>
Figure QLYQS_13
振动加速度信号的平均值/>
Figure QLYQS_14
所述平均值为新坐标系中的沿/>
Figure QLYQS_10
轴方向的振动加速度信号平均值；通过公式
Figure QLYQS_12
得到10000个/>
Figure QLYQS_15
振动加速度信号的平均值/>
Figure QLYQS_16
该平均值为新坐标系中的沿/>
Figure QLYQS_9
轴方向的振动加速度信号平均值；
将
Figure QLYQS_17
与虚数单位j相乘后与/>
Figure QLYQS_18
相加，构造出复数坐标系内的振动加速度复数信号矢量/>
Figure QLYQS_19
用以表示复数坐标系内的振动加速度信号平均值的矢量形式，代表故障无关的常矢量；在得到的10000个矢量/>
Figure QLYQS_20
中分别有一个常矢量/>
Figure QLYQS_21
从各矢量/>
Figure QLYQS_22
中分别减去该常矢量，分别得到10000个复数坐标系内与故障有关的旋转矢量/>
Figure QLYQS_23
步骤3，确定振动旋转矢量的模长序列：
对得到的旋转矢量
Figure QLYQS_24
分别通过实部平方与虚部平方之和开二次方的运算，进行取模，得到的模为：
Figure QLYQS_25
所述取的模为10000个，该10000个模分别对应ti各时刻；
各模长均为L(ti)，分别是L(t1),L(t2)……L(t9998),L(t9999),L(t10000)；各模长值分别代表相应时刻的振动加速度的幅值；该模长值单位为g，为标量；
将得到系列模长值按时间顺序排列在加速度振幅--时间图上：
t1＝0.0001s，t2＝0.0002s......t9998＝0.9998s，t9999＝0.9999s，t10000＝1.0000s
所述加速度幅值--时间图的横坐标是时间t，单位为s；纵坐标是加速度幅值A，亦即模长值，该加速度幅值的单位是g；得到一系列横坐标间隔为0.0001s的离散的点；将相邻的点用直线连接，得到连续函数曲线；
步骤4，确定连续函数曲线中的两端的用于线性插值的极大值点坐标：
所述用于线性插值的极大值点是分别位于所述得到连续函数曲线两端的各两个极大值；
具体分别是该连续函数曲线中的按时间顺序从小到大排列的第一个极大值点、第二个极大值点，以及该按时间顺序排列的倒数第一个极大值点、倒数第二个极大值点；
步骤5，通过线性插值得出区间端点处的极大值点坐标：
在所述P1(0.0005，4870.84)与P2(0.0009，3638.49)之间构造直线，在加速度振幅--时间图上，A＝－3080875t+6411.2775；当t＝0.0001s时，代入A中得到A(0.0001)＝6103.19g；得到采样区间t＝0.0001s时刻的原图像中的极大值点，该极大值点的坐标为(0.0001，6103.19)；
在所述Ps-1(0.9991，127.579)与Ps(0.9996，103.047)之间构造直线，在加速度振幅--时间图上，为：A＝－49064t+49147.4214；当t＝1s时，代入A中得到A(1)＝83.4214g；得到采样区间t＝1.0000s时刻的原图像中的极大值点，该极大值点的坐标为(1，83.4214)；
步骤6，确定采样区间内的三次样条插值函数：
利用三次样条插值法，对采样区间内的各极大值点的坐标以及步骤5中插值得到的区间两端端点处的极大值点的坐标进行插值，得到三次样条插值函数A(t)；
步骤7，确定矢量模极大值包络谱：
对步骤6得到的三次样条插值函数A(t)进行傅里叶变换，得到矢量模极大值包络谱；具体是：
对区间[t0,ts+1]内，进行傅里叶变换：
Figure QLYQS_26
其中，A表示按所述数据处理方法得到的滚动轴承的加速度幅值，单位是g；f表示按所述数据处理方法得到的滚动轴承的频率，单位是Hz；对步骤7得到的三次样条插值函数A(t)进行傅里叶变换后，将振动加速度振幅--时间曲线变换为振动加速度振幅--频率曲线；F[A(t)]表示对A(T)进行傅里叶变换；得到矢量模极大值包络谱；
步骤8，结合矢量模极大值包络谱对滚动轴承故障进行诊断与分析；
根据得到的所述矢量模极大值包络谱，在滚动轴承外环故障特征频率131Hz及滚动轴承外环故障特征频率两侧的边频116Hz与146Hz处具有较大加速度值，在滚动轴承外环故障特征频率的二倍频261Hz及滚动轴承外环故障特征频率两侧的二倍频的边频246Hz与276Hz处，以及滚动轴承外环故障特征频率的三倍频392Hz及滚动轴承外环故障特征频率的三倍频两侧的边频377Hz与407Hz处，均有较大的加速度值；判定该滚动轴承外环处存在损伤故障。
2.如权利要求1所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，所述构造10000个复数信号矢量
Figure QLYQS_27
的具体过程是：
在滚动轴承圆周表面布置两个相同型号的加速度传感器，并使各传感器之间的夹角为25°～155°；以采样频率为10kHz的等时间采样方法，当滚动轴承外环静止且内环转速为2400rpm时，在3s时段内两个所述加速度传感器分别测得振动加速度信号
Figure QLYQS_28
与/>
Figure QLYQS_29
单位为g，且g＝9.8m/s2；其中/>
Figure QLYQS_30
Figure QLYQS_31
中各时刻分别与0s,0.0001s,0.0002s…2.9998s,2.9999s,3.0000s一一对应；选取3s时间段中的1.0001s时刻作为用于处理数据时的0.0001s时刻；选取3s时间段中的2s时刻作为用于处理数据时的1s时刻，以此排除接通传感器与关闭传感器时对数据采集的影响；故用于分析处理的采样数据包含的时长为0.9999s；
在所述0.9999s的时间段内，各加速度传感器能够测得在0.0001s,0.0002s,0.0003s…0.9998s,0.9999s,1.0000s共1000个时刻的振动加速度数据共计10000对；分别是x(0.0001)，x(0.0002)，x(0.0003)…x(0.9998)，x(0.9999)，x(1.0000)，记为x(ti)；以及y(0.0001)，y(0.0002)，y(0.0003)…y(0.9998)，y(0.9999)，y(1.0000)，记为y(ti)；所述ti中，i＝为1，2，3…9998，9999，10000；各ti时刻对应的时刻分别是0.0001s，0.0002s，0.0003s，……，0.9998s，0.9999s，1.0000s。
3.如权利要求2所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，所述两个传感器之间的夹角为25°～155°；各传感器的采样方法为等时间采样法；各传感器的采样频率均为10kHz。
4.如权利要求2所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，当布置的两个加速度传感器的夹角为90°时，则将x(ti)直接记为
Figure QLYQS_32
表示10000个新坐标系内的水平方向的振动加速度信号，并将y(ti)直接记为/>
Figure QLYQS_33
表示10000个新坐标系内的竖直方向的振动加速度信号；
当两个加速度传感器二者之间的夹角θ≠90°时，通过公式(1)与(2)进行正交变换：
Figure QLYQS_34
Figure QLYQS_35
其中，式(1)是将原坐标系中的沿y轴方向的10000个振动加速度信号y(ti)作为新坐标系中的沿
Figure QLYQS_36
轴方向的10000个振动加速度信号/>
Figure QLYQS_37
式(2)是将原坐标系中的沿x轴方向的10000个振动加速度信号x(ti)变换为新坐标系中的沿/>
Figure QLYQS_38
轴方向的10000个振动加速度信号/>
Figure QLYQS_39
将新坐标系中的两方向相互垂直的10000对振动加速度信号
Figure QLYQS_40
与/>
Figure QLYQS_41
通过/>
Figure QLYQS_42
与虚数单位j相乘，再与/>
Figure QLYQS_43
加和，构造成10000个复数坐标系内的振动加速度复数信号矢量/>
Figure QLYQS_44
用以表征10000个变换到复数坐标系里的振动加速度复数信号矢量形式。
5.如权利要求1所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，步骤4中得到所述连续函数曲线中的按时间顺序从小到大排列的第一个极大值点的坐标P1、第二个极大值点的坐标P2，以及倒数第一个极大值点的坐标Ps-1、倒数第二个极大值点的坐标Ps；分别为：
P1：(0.0005，4870.84)，
P2：(0.0009，3638.49)；
Ps-1：(0.9991，127.579)，
Ps：(0.9996，103.047)。
6.如权利要求5所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，对得到的连续函数曲线根据数学上对极大值点的定义得到s个极大值点，分别是P1，P2，P3...Ps-2，Ps-1，Ps；各所述极大值点的坐标分别为：P1(t1，A(t1))，P2(t2，A(t2))，P3(t3，A(t3))...Ps-2(ts-2，A(ts-2))，Ps-1(ts-1，A(ts-1))，Ps(ts，A(ts))；其中，0.0001s＜t1＜t2＜t3...＜ts-2＜ts-1＜ts＜1s。
7.如权利要求1所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，步骤6中所述确定采样区间内的三次样条插值函数的具体过程是：
I确定采样区间内的极大值点与区间端点处的极大值点均为三次样条插值所用的数据节点；
II计算步长hi＝ti+1-ti；其中，i＝0，1，2…s；
III将s+2个数据节点的坐标值，与自然边界条件：“m0＝0及ms+1＝0”入如下矩阵中：
Figure QLYQS_45
并解出得到s+2个数据节点处的二阶导数的值：m0，m1，m2……，ms+1；
IV计算样条曲线的系数；
V在各子区间[ti，ti+1]内构造曲线Ai(t)＝ai+bi(t-ti)+ci(t-ti)2+di(t-ti)3；则各子区间[ti，ti+1]，i＝0，1，2…s内的共计s+1条曲线：A0(t)，A1(t)...As-1(t)，As(t)在采样区间内s个极大值点的坐标处是二阶连续的；将所述s+1条在各子区间[ti，ti+1]内构造曲线记为函数A(t)，其中，t的取值范围：[t0，ts+1]；所述函数A(t)是通过三次样条插值方法得到的三次样条插值函数。
8.如权利要求7所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，所述采样区间内的极大值点个数为s，步骤5中通过插值得到的区间端点处的极大值点个数为2，故所述数据节点的个数为s+2，按横坐标时间顺序从小到大记为(t0，A0)，(t1，A1)，(t2，A2)……(ts+1，As+1)；各所述数据节点在横坐标方向上，将采样区间[0.0001s，1.0000s]分成了s+1个子区间：(t0，t1)，(t1，t2)，(t2，t3)...(ts，ts+1)；且三次样条插值函数的边界条件为自然边界条件，即区间两端端点的二阶导数均为0。
9.如权利要求7所述矢量模极大值包络的滚动轴承故障诊断方法，其特征在于，所述计算样条曲线的系数的具体过程是：
ai＝Ai是各子区间中三次样条插值函数的常数项；
Figure QLYQS_46
是各子区间[ti，ti+1]中三次样条插值函数的一次方项系数；i＝0,1,2…s；
Figure QLYQS_47
是各子区间[ti，ti+1]s中，三次样条插值函数的平方项系数；i＝0,1,2…s；
Figure QLYQS_48
是各子区间[ti，ti+1]中三次样条插值函数的立方项系数；i＝0,1,2…s。
CN202210108671.4A 2022-01-28 2022-01-28 一种矢量模极大值包络的滚动轴承故障诊断方法 Active CN114486252B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210108671.4A CN114486252B (zh) 2022-01-28 2022-01-28 一种矢量模极大值包络的滚动轴承故障诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210108671.4A CN114486252B (zh) 2022-01-28 2022-01-28 一种矢量模极大值包络的滚动轴承故障诊断方法

Publications (2)

Publication Number Publication Date
CN114486252A CN114486252A (zh) 2022-05-13
CN114486252B true CN114486252B (zh) 2023-06-30

Family

ID=81475840

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210108671.4A Active CN114486252B (zh) 2022-01-28 2022-01-28 一种矢量模极大值包络的滚动轴承故障诊断方法

Country Status (1)

Country Link
CN (1) CN114486252B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115270896B (zh) * 2022-09-28 2023-04-07 西华大学 一种用于识别航空发动机主轴承松动故障的智能诊断方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110806315A (zh) * 2019-11-20 2020-02-18 北京工业大学 一种基于倒位编辑的齿轮箱复合故障诊断方法
CN111307460A (zh) * 2020-03-14 2020-06-19 中国石化销售股份有限公司华南分公司 基于计算阶次跟踪与谱峭度的滚动轴承故障诊断方法
AU2020103681A4 (en) * 2020-11-26 2021-02-04 Anhui University Of Technology Rolling Bearing Fault Diagnosis Method Based on Fourier Decomposition and Multi-scale Arrangement Entropy Partial Mean Value
CN112577746A (zh) * 2020-12-07 2021-03-30 东南大学 一种转速波动下滚动轴承包络阶次谱故障特征的提取方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110806315A (zh) * 2019-11-20 2020-02-18 北京工业大学 一种基于倒位编辑的齿轮箱复合故障诊断方法
CN111307460A (zh) * 2020-03-14 2020-06-19 中国石化销售股份有限公司华南分公司 基于计算阶次跟踪与谱峭度的滚动轴承故障诊断方法
AU2020103681A4 (en) * 2020-11-26 2021-02-04 Anhui University Of Technology Rolling Bearing Fault Diagnosis Method Based on Fourier Decomposition and Multi-scale Arrangement Entropy Partial Mean Value
CN112577746A (zh) * 2020-12-07 2021-03-30 东南大学 一种转速波动下滚动轴承包络阶次谱故障特征的提取方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
降速工况下滚动轴承微弱故障特征信号提取新方法;栾孝驰;沙云东;;机械设计与制造(第03期);全文 *

Also Published As

Publication number Publication date
CN114486252A (zh) 2022-05-13

Similar Documents

Publication Publication Date Title
US9341511B2 (en) Timing analysis
Shi et al. Bearing fault diagnosis under variable rotational speed via the joint application of windowed fractal dimension transform and generalized demodulation: A method free from prefiltering and resampling
Wang et al. Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis
CN103471841B (zh) 一种旋转机械振动故障诊断方法
CN110987438B (zh) 水轮发电机变转速过程周期性振动冲击信号检测的方法
CN110567574A (zh) 一种旋转叶片叶端定时振动参数辨识方法与系统
CN109855874B (zh) 一种声音辅助振动微弱信号增强检测的随机共振滤波器
JP2003528292A (ja) 振動解析によるベアリングの状態ベースのモニタリング
Wang et al. Sparse and low-rank decomposition of the time–frequency representation for bearing fault diagnosis under variable speed conditions
CN104215456B (zh) 一种基于平面聚类和频域压缩感知重构的机械故障诊断方法
Lin et al. A review and strategy for the diagnosis of speed-varying machinery
CN112461934B (zh) 一种基于声发射的航空发动机叶片裂纹源定位方法
CN102721462B (zh) 旋转机械启停车过程波德图/奈奎斯特图的快速计算方法
CN114486252B (zh) 一种矢量模极大值包络的滚动轴承故障诊断方法
Wu et al. A modified tacho-less order tracking method for the surveillance and diagnosis of machine under sharp speed variation
Zhao et al. Rolling element bearing instantaneous rotational frequency estimation based on EMD soft-thresholding denoising and instantaneous fault characteristic frequency
Choudhury et al. An overview of fault diagnosis of industrial machines operating under variable speeds
CN114330489A (zh) 一种监测设备故障诊断方法及系统
CN112465068A (zh) 一种基于多传感器数据融合的旋转设备故障特征提取方法
Chen et al. A time-varying instantaneous frequency fault features extraction method of rolling bearing under variable speed
CN112345247B (zh) 一种滚动轴承的故障诊断方法及装置
Shang et al. Rolling bearing fault diagnosis method based on MOMEDA and IEWT
CN112781723B (zh) 一种基于频谱方差的谐波成分检测方法
Song et al. Multispectral Balanced Automatic Fault Diagnosis for Rolling Bearings under Variable Speed Conditions
CN115326396A (zh) 一种轴承故障的诊断方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant