CN105403407B

CN105403407B - 一种基于正常域估计的列车滚动轴承隐患辨识方法

Info

Publication number: CN105403407B
Application number: CN201510670716.7A
Authority: CN
Inventors: 张媛; 杜艳平; 李宏峰; 秦勇; 贾利民; 邢宗义
Original assignee: Beijing Institute of Graphic Communication
Current assignee: Beijing Institute of Graphic Communication
Priority date: 2015-10-13
Filing date: 2015-10-13
Publication date: 2017-09-22
Anticipated expiration: 2035-10-13
Also published as: CN105403407A

Abstract

本发明提供一种基于正常域估计的列车滚动轴承隐患辨识方法，通过采集列车滚动轴承在正常运行状态下的振动加速度数据、分段预处理所采集的数据、提取统计状态特征向量，并创新性的提出正常域的概念且凸包估计出正常域的边界，最后判断统计状态特征向量对应的实时状态点是否位于正常域边界内部来进行隐患辨识。本发明能够在仅有正常状态数据的情况下准确地辨识列车滚动轴承是否存在隐患，为滚动轴承的状态监测、隐患评估、故障预警提供了有效的方法。

Description

一种基于正常域估计的列车滚动轴承隐患辨识方法

技术领域

本发明属于及轨道交通安全技术领域。本发明涉及一种基于正常域估计的无故障数据环境下滚动轴承隐患辨识方法。

背景技术

据统计，火车是事故发生率最低的交通工具。可想而知，为保障安全运营，从一列新车上线直至服役期满下线，其绝大部分时间必然是处于正常服役状态的，尽管偶有关键部件发生故障，但相对来说故障发生的概率很小。尤其是随着机车车辆制造技术的不断发展，列车本身的质量越来越好，安全性也越来越高，其在线服役时关键部件发生异常的情况越来越少，且现代化的维修和养护手段也为列车运行时的高安全性提供了有效保障。因此，在某一关键装备的状态数据积累方面，面临的情况是：所采集的数据绝大部分是正常状态下的数据，相对而言有效的故障状态数据可能十分稀少甚至没有。在新车上线或旧装备换新时，这种情况更为普遍，甚至可能存在运行数月都没有故障或异常发生的情况。鉴于此，有效的无故障数据环境下隐患辨识方法和技术亟待研究和突破。

目前，国内外已有部分学者在隐患辨识相关的早期故障诊断和故障预测方面进行了研究。现有研究成果从研究手段上主要分为如下三大类。

(1)是基于全寿命状态数据和剩余寿命分析和早期故障检测，如：Chen Lin等针对齿轮箱的早期故障诊断问题，利用连续时间马尔科夫模型进行各种工况状态的模拟，并基于振动信号的向量回归建模和贝叶斯方法计算早期故障发生概率；任丽娜等研究了数控机床故障过程，提出了4参数非齐次泊松过程模型用于数控机场可靠性分析和早期故障预测。

(2)利用基于状态转移或过程模型的概率分析结果进行早期故障检测：JustynaPetke等提出了一种普适的基于多参数和多状态数据的组合交互式测试方法，可识别多种可能发生的早期故障，并将这些故障按发生概率进行排序；许丽佳等针对电子系统的故障预测问题，提出了线性辨别分析与马尔科夫链以及贝叶斯网络相结合的方法。

(3)基于实际的或仿真的早期故障数据，利用各种信号处理和模式识别的方法进行早期故障诊断，这类方法的研究成果最为丰富，如：Ilhan Aydin等针对旋转机械的早期故障诊断问题，采集了早期故障状态的振动数据，提出了小波和经验模式分解相结合的诊断方法；Henry David等针对飞机控制面伺服环路的早期故障，提出了基于模型的线性变参数的检测方法，并利用空客飞机的仿真数据进行了验证。

上述已有的隐患监测和早期故障识别相关研究中，几乎都采用了实验或仿真技术进行了故障建模或采集故障数据，无法直接应用于无故障数据环境。因此，针对轨道交通列车这种无法获取故障状态数据的特殊对象，采用常用的数据驱动及传统的模式识别方法无法十分有效地完成列车滚动轴承的定量化状态监测和隐患识别。鉴于此，本发明提出了基于正常域估计的无故障数据环境下滚动轴承的隐患辨识方法。

发明内容

本发明所采用的技术方案是：

本发明的目的是针对无故障数据环境下的列车滚动轴承隐患辨识问题，提出了基于正常域估计的隐患辨识方法。本发明创新性的提出正常域估计的概念，正常域是针对某一具体的研究对象(如列车上某些关键装备)，在研究对象的安全相关变量空间内，能够且仅包含正常运行状态数据的特征点的区域。直观来说，正常域由其边界确定，其边界可理解为能够包住所有正常运行状态点的最小的闭合的几何形状，即：在二维的安全相关变量空间内，正常域为能够包围所有正常状态点的最小闭合曲线；在三维的安全相关变量空间内，正常域为能够包住所有正常状态点的最小闭合曲面；在更高维的空间内，正常域则为能够包住所有正常状态点的最小超平面。附图1所示即为二维空间内的正常域示意图，其中能包围住所有正常状态点的最小闭合曲线即为正常域边界。

本发明可以为缺乏故障数据情况下的列车滚动轴承隐患辨识提供一种新方法，为解决新车或无故障历史车辆的状态监测提供技术支持和指导。具体采用如下技术方案：该方法包括如下步骤：

1)采集列车滚动轴承在正常运行状态下的振动加速度数据；

2)预处理所采集的数据：按时间间隔t_I对采集到的振动加速度数据进行分段，对划分完成的每段数据进行局部均值处理，获得每段数据的乘积函数矩阵，所述时间间隔t_I可在开区间0～T_c内任意取值，其中T_c为振动加速度数据采集时间；

3)提取滚动轴承正常运行状态的统计状态特征向量：对每个乘积函数矩阵利用主元分析方法进行处理，计算其T²和SPE统计量的控制限，两个控制限的值构成了每段数据的二维统计状态特征向量，每个向量对应一个状态点；将获得的T²和SPE统计量控制限数据分别进行归一化处理；

4)估计正常域边界：将正常域用凸包进行形式化描述，利用凸包估计中的Jarvis算法求解能够包围3)中统计状态特征点集的最小凸包，此凸包即为估计出的正常域边界；正常域，是指研究对象的安全相关变量空间内，能够且仅包含正常运行状态数据的特征点的区域，正常域边界为能包围住所有正常状态点的最小闭合曲线；

5)进行滚动轴承的隐患辨识：将实时采集的振动加速度数据经步骤2)～3)的处理可获得统计状态特征向量，判断统计状态特征向量对应的实时状态点是否位于步骤4)正常域边界内部，若位于正常域内，则表示当前状态正常，若位于正常域外，则表示出现隐患。

进一步地，步骤4)中利用凸包估计正常域边界的具体步骤如下：

①设所有正常状态点组成点集P，找出点集P中y轴坐标最小的点，若存在多个这样的点，则取最左边的点记为p₀，该点必为凸包的顶点；

②从点p₀向右引一条平行于x轴的射线，即为l₀；

③沿着点p₀逆时针旋转射线l₀，直到直线l₀再次与点集P内的点相交，若这样的交点有2个或2个以上，则按y轴坐标升序排列这些点，记为p₁,p₂,,...,p_m，p₁,p₂,,...,p_m必为凸包上的顶点；记过p₀,p₁,p₂,,...,p_m的直线为p₀p_m；

④以p_m为中心点，重复步骤③直到重新回到点p₀为止，即可求得点集P的凸包的全部顶点。

进一步地，步骤2)中获得乘积函数矩阵具体方法如下：将每段数据进行基于三次样条函数的局部均值分解，分解过程如下步骤①～⑧，得到每段数据的乘积函数分量后，由各个乘积函数分量构成乘积函数矩阵；

①设原始信号为x(t)，找出其所有局部极值点，将所有的局部极大值点和局部极小值点分别用三次样条曲线连接起来，得到x(t)的上包络线E_max(t)和下包络线E_min(t)；

②按如下两式计算局部均值函数m₁₁(t)和包络估计函数c₁₁(t)；

③将局部均值函数m₁₁(t)从原始信号x(t)中分离，得

h₁₁(t)＝x(t)-m₁₁(t)

④用h₁₁(t)除以包络估计函数c₁₁(t)，以对h₁₁(t)进行解调，得

⑤按照①中方法求出s₁₁(t)所对应的包络估计函数c₁₂(t)，若c₁₂(t)＝1，则s₁₁(t)为纯调频函数，若c₁₂(t)≠1，则需重复上述①～④迭代过程，直至s_1n(t)的包络估计函数c_1(n+1)(t)＝1；故有

理论上，迭代终止条件为

⑥把迭代过程中产生的所有包络估计函数相乘可得包络信号

⑦将包络信号c₁(t)＝1和纯调频信号s_1n(t)相乘，即得x(t)的第一个乘积函数分量PF₁

PF₁(t)＝c₁(t)·s_1n(t)

⑧将PF₁从x(t)中分离，得到一个新的信号r₁(t)，将此信号作为待分解信号重复上述①～⑦步骤k次，直至r_k(t)为一单调函数为止，即

则如下所示，原始信号x(t)被分解为k个PF分量和一个r_k(t)之和：

本发明具有如下有益的技术效果：创新地提出正常域估计的概念，能够在仅有正常状态数据的情况下准确地辨识列车滚动轴承是否存在隐患，为滚动轴承的状态监测、隐患评估、故障预警提供了有效的方法。

附图说明

图1是正常域示意图。

图2是凸包示意图。

图3正常域边界估计结果。

图4存在隐患时的辨识结果。

具体实施方式

下面结合附图对本发明的实施例具体说明。一种基于正常域估计的无故障数据环境下列车滚动轴承隐患辨识方法包括以下具体实施步骤：

1)采集正常运行状态数据：采集列车滚动轴承在正常运行状态时的振动加速度数据；

2)振动加速度数据的分段处理：将所采集到的振动加速度数据按固定的时间间隔t_I进行分段，划分得到的每个数据段将对应一个状态特征向量，即状态点；

3)获得乘积函数矩阵：将每段数据进行基于三次样条函数的局部均值分解，分解过程如下步骤①～⑧，得到每段数据的乘积函数分量后，由各个乘积函数分量构成乘积函数矩阵；

③将局部均值函数m₁₁(t)从原始信号x(t)中分离，得

h₁₁(t)＝x(t)-m₁₁(t) (3)

⑤按照①中方法求出s₁₁(t)所对应的包络估计函数c₁₂(t)，若c₁₂(t)＝1，则s₁₁(t)为纯调频函数，若c₁₂(t)≠1，则需重复上述①～④迭代过程，直至s_1n(t)的包络估计函数c_1(n+1)(t)＝1。故有

理论上，迭代终止条件为

⑥把迭代过程中产生的所有包络估计函数相乘可得包络信号

⑦将包络信号c₁(t)＝1和纯调频信号s_1n(t)相乘，即得x(t)的第一个乘积函数(product function,PF)分量PF₁

PF₁(t)＝c₁(t)·s_1n(t) (8)

则如下所示，原始信号x(t)被分解为k个PF分量和一个r_k(t)之和。

4)提取统计状态特征向量：在获得正常状态下的滚动轴承振动数据的乘积函数矩阵后，基于主元分析的方法，选取最常用的T²和SPE两统计量，分别计算此两个统计量的值及其控制限，并提取两个统计量的控制限的值作为统计状态特征量。T²和SPE统计量的控制限计算按如下步骤①～④进行：

①建立多元统计模型

将某段数据的乘积函数矩阵表示为数据集Y，Y＝[c₁ c₂ … c_n]，则对数据集Y_a×b(a为样本个数，b为变量个数)中每一时刻的数据向量(表示实数域)按下式进行标准化

式中，E(y)＝[μ₁,μ₂,…,μ_b]^T为y对应的均值向量，D_σ＝diag(σ₁,σ₂,…,σ_b)为方差矩阵，为第j个变量的标准差，j＝1,2…b。

记标准化后的数据集为对的相关系数矩阵作奇异值分解

R＝UD_λU^T (12)

式中，为一酉矩阵，D_λ＝diag(λ₁,λ₂,…,λ_b)为一对角阵。在新的坐标系U的各个方向上的方差满足λ₁＞λ₂＞…＞λ_b。称U的前d(d<b)维线性无关向量P＝[u₁,u₂,…,u_d]构成的子空间为主元空间后b-d维向量P'＝[u_d+1,u_d+2,…,u_b]构成的子空间为残差空间主元个数d通常采用方差累计贡献率法确定。则数据向量可分解为

式中，分别为在和上的投影。

②在主元空间中建立T²统计量并计算其控制限

T²统计量的定义为

式中，为D_λ的前d个对角元素组成的矩阵，为主元打分向量。对于样本个数为a，主元个数为d的数据向量T²服从自由度为d和a-d的F分布，即

式中，F(d,a-d)为自由度d和a-d的中心F分布。则置信度为α的T²统计量控制限T² _CL为

式中，F_α(d,a-d)为自由度d和a-d的中心F分布的上100α百分位点，其值可由F分布表查得。本实施例取常用置信度水平α＝0.95。

③在残差空间中建立SPE统计量并计算其控制限

SPE统计量定义为

当检验水平为α时，SPE的控制限SPE_CL为

式中，θ₁＝λ_d+1+λ_d+2+…+λ_b，C_α为标准正态分布的100α百分位点，(其中与②中相同，α＝0.95

④保存T² _CL和SPE_CL值所确定的数据样本点为正常状态点

将T² _CL和SPE_CL变量分别作为横坐标和纵坐标，构成一个二维平面空间，每段数据获得一个T² _CL值和一个SPE_CL值，从而可在二维平面上构成一个数据点(滚动轴承的正常状态点)，分别保存各段数据对应的此数据点以用于后续的正常域估计和隐患辨识。

5)完成正常域的形式化描述：利用凸包的定义对二维正常域边界进行形式化描述，如下定义1～定义3中的描述：

定义1：设集合若对于和都有αx₁+(1-α)x₂∈S,则称S是凸集。

定义2：设x₁,x₂,...,x_k∈Rⁿ，如果存在满足且α_i≥0使得则称x是x₁,x₂,...,x_k的一个凸组合，其中α₁,α₂,...,α_k是相应的凸组合系数。

定义3：设集合且S＝{x₁,x₂,...,x_k}，则将包含S的所有的凸集的交集成为S的凸包，记为co(S)。S的凸包co(S)是包含S的最小凸集，且可以由集合S内所有点的凸组合构造而成，即

平面点集的凸包是指包含平面点集内所有点并且顶点属于平面点集的最小简单凸多边形，可形象地将其想象为一条刚好包围所有点的橡皮圈，如附图2所示。二维正常域边界即为状态点在二维平面上时平面点集的凸包。

6)进行正常域边界估计：基于Jarvis算法完成正常域边界估计，凸包计算按如下步骤①～④进行；

②从点p₀向右引一条平行于x轴的射线，即为l₀；

③沿着点p₀逆时针旋转射线l₀，直到直线l₀再次与点集P内的点相交，若这样的交点有2个或2个以上，则按y轴坐标升序排列这些点，记为p₁,p₂,,...,p_m，p₁,p₂,,...,p_m必为凸包上的顶点。记过p₀,p₁,p₂,,...,p_m的直线为p₀p_m；

7)进行滚动轴承的隐患辨识：实时采集列车滚动轴承的振动加速度数据，经实施步骤2)～4)的处理可获得统计状态特征向量(即状态点)，将实时状态点与正常域边界进行对比，若状态点位于正常域内，则表示当前状态正常，若位于正常域外，则表示出现隐患。

以下是本发明的一个实施例。

实施例中所需数据是由Dr.Kenneth A.Loparo提供的滚动轴承实验数据，轴承型号为205-2RS JEM SKF型深沟球轴承，电机负载3马力，转速1730r/min(约28.8r/s)，振动加速度数据采集点事负载端，采样频率12k Hz，共采样4次，每次采样时间为10s；采样时间10s。

本实施例将4次正常状态下的采样数据按照前述具体实施步骤提取T²控制限和SPE控制限的二维统计特征，形成二维平面上的状态点集；利用Jarviss算法进行凸包计算并绘图。试验结果如下所示，附图3为所获正常状态点集的凸包计算结果，即正常域边界估计结果。

为验证正常域边界估计的有效性，本实施例中利用滚动轴承存在隐患时(滚子表面存在一个直径为0.1778mm深度0.2794mm的凹坑)的振动加速度数据，仍按照前述具体实施步骤中的2)～4)计算了存在隐患时的状态点，由附图4可见，该状态点落了在正常域之外，此结果表明本发明所提出的基于正常域估计的无故障数据环境下列车滚动轴承隐患辨识方法是有效的。

Claims

1.一种基于正常域估计的列车滚动轴承隐患辨识方法，其特征在于，包括如下步骤：

1)采集列车滚动轴承在正常运行状态下的振动加速度数据；

2.根据权利要求1所述的方法，其特征在于：步骤4)中利用凸包估计正常域边界的具体步骤如下：

①设所有正常状态点组成点集P，找出点集P中y轴坐标最小的点，若存在多个这样的点，则取最左边的点记为p₀，该点设为凸包的顶点；

②从点p₀向右引一条平行于x轴的射线，记为l₀；

③沿着点p₀逆时针旋转射线l₀，直到直线l₀再次与点集P内的点相交，若这样的交点有2个或2个以上，则按y轴坐标升序排列这些点，记为p₁,p₂,,...,p_m，p₁,p₂,,...,p_m设为凸包上的顶点；记过p₀,p₁,p₂,,...,p_m的直线为p₀p_m；

④以p_m为中心点，重复步骤③直到重新回到点p₀为止，求出点集P的凸包的全部顶点，最后将全部顶点连成封闭的曲线，该曲线为正常域边界。

3.根据权利要求1所述的方法，其特征在于：步骤2)中获得乘积函数矩阵具体方法如下：将每段数据进行基于三次样条函数的局部均值分解，分解过程如下步骤①～⑧，得到每段数据的乘积函数分量后，由各个乘积函数分量构成乘积函数矩阵；

③将局部均值函数m₁₁(t)从原始信号x(t)中分离，得

h₁₁(t)＝x(t)-m₁₁(t)

理论上，迭代终止条件为

⑥把迭代过程中产生的所有包络估计函数相乘可得包络信号

<mrow> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>c</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msub> <mi>c</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>...</mo> <msub> <mi>c</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>c</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow>

PF₁(t)＝c₁(t)·s_1n(t)

则如下所示，原始信号x(t)被分解为k个乘积函数分量和一个r_k(t)之和：

<mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>v</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <msub> <mi>PF</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow> 2