CN105241680A - 一种基于概率密度函数的旋转机械健康状态评估方法 - Google Patents

Abstract

本发明涉及一种基于核密度估计和K-L散度等统计学方法来实现旋转机械的健康状态评估方法。包括如下几个步骤：步骤1：采集被监测对象的原始振动数据；步骤2：从原始振动数据中提取时、频域特征；步骤3：对时、频域特征进行降维处理获得敏感特征；步骤4：使用窗宽为k的移动滑动窗来动态地选定样本集，计算样本集合的统计特性；步骤5：计算各个样本集在每个敏感特征下的概率密度函数；步骤6：计算相同敏感特征下相邻的两个样本集概率密度函数的K-L散度值；步骤7：计算集成K-L散度作为监测对象的健康评估指标。本发明的有益效果在于：充分考虑到样本的统计不确定性，提高了健康状态评估模型的准确性和泛化性。

Description

一种基于概率密度函数的旋转机械健康状态评估方法

技术领域

本发明属于机械装备的性能退化评估领域，具体涉及一种基于核密度估计和K-L散度等统计学方法以实现旋转机械的健康状态评估。

背景技术

旋转机械广泛应用于大型制造系统和重要技术装备中，例如数控机床、风力发电机、航空航天发动机等。在这些制造系统和装备中，滚动轴承、齿轮等旋转机械中的关键部件发挥着重要作用，但是旋转机械零件在运转中需要承受交变的机械应力和偶然冲击，加上本身固有的制造误差，经常会产生一些早期缺陷(如：轻度磨损、点蚀等)。如果这些缺陷不及时诊断发现，就会随着时间的推移不断恶化，最终导致系统失效，带来巨大的财产损失。对于现代化大型复杂旋转机械设备，尽管可以通过改善设计、制造工艺来提高零部件的质量，但仍然不能完全避免故障的发生。健康状态评估和故障诊断已经成为保障旋转机械装备正常运转和预防故障发生的重要技术手段。通过这两种技术的结合，可以尽早地发现设备运行过程中的早期故障，从而通过对早期故障加以隔离，避免设备的突发性失效和意外停机。健康状态评估的首要工作是利用先进的传感和监测技术，对关键零部件及系统进行状态监测。为了监测旋转机械设备的健康状况，需要对该设备的运行状态进行准确地评估，然后根据评估结果制定合理的维修计划。使用敏感的健康指标探测机械设备异常状态和正常状态的分界点或分界面，不但有助于预防机械产品发生严重的功能失效，而且可以避免灾难性失效引起的巨额财产损失。在旋转机械健康状态评估中，选择合适的健康评估指标和建立健康状态评估模型尤为重要。

基于数据驱动的健康状态评估和故障诊断方法是近年来逐渐兴起的一种新技术，计算机技术的快速发展使得大数据并行高速计算变得非常容易，推动了依靠大量数据分析的健康评估技术的发展。从应用的角度看，基于数据驱动的健康状态评估和故障诊断方法比基于物理模型的方法更为切实可行，这是由于数据采集通常要比精确建立物理模型更加容易。除此之外，基于数据驱动的健康评估方法还有两个明显的优点：一是该类方法更容易实现自动化评估，这与现代工业的智能化发展是切合的；二是该类方法不需要太多参数设置和专家经验知识，所以具有一定的稳健性。一般来说，一种数据驱动的健康评估方法应包括数据获取、特征提取、特征选择或特征变换、选择评估指标及设计评估模型、评估结果输出五个步骤。其中，评估指标的选择及评估模型的设计是该类方法的关键。

目前，现有的基于数据驱动的健康评估方法大多都是建立在数据融合的基础上，然后通过一些比较成熟的机器学习等智能技术进行回归预测，当选取的健康评估指标达到预先设定的阈值时，则认为所监测的对象发生故障或者达到失效状态。例如，基于支持向量机(SupportVectorMachine,SVM)的健康状态评估方法，基于人工神经网络(ArtificialNeuralNetwork,ANN)的健康状态评估方法等。然而由于噪声测量误差等对有效振动信号的污染，使得回归及分类问题存在一定的不确定性，进而导致了结果误差较大，即很难找到一个恰当的超平面使得所有样本落在允许的较小误差带范围内。

由以上常用方法得到的分类及回归预测结果误差较大的原因与它们的训练过程的原理有关。现有的健康状态评估方法往往忽略了样本间的统计信息和关联信息，而统计信息对于随机信号处理是极为关键的，即对于正确分类及回归预测很有帮助。目前，国内外从样本统计性角度开展的智能故障诊断方法的研究还非常有限。

发明内容

本发明的目的是为了从原始样本中提取更加全面、有效的统计学信息，从而提高回归预测的准确率和推广能力，提出一种基于概率密度函数的旋转机械健康状态评估方法。

本发明利用核密度估计(KernelDensityEstimation，KDE)和K-L散度(Kullback-LeiblerDivergence)两种统计学工具实现健康状态的评估。因而本发明一种基于概率密度函数的旋转机械健康状态评估方法，包括如下步骤：

步骤1：采集被监测对象的原始振动数据；

步骤2：从步骤1中得到的原始振动数据中提取多个的时域和频域特征，获得时域特征集和频域特征集；

步骤3：对步骤2中得到的时域和频域特征集进行降维处理获得敏感特征，并且计算这些敏感特征的贡献率，然后获得贡献率向量；

步骤4：采集被监测对象每一个时刻的原始振动数据为一个样本，使用一个窗宽为k的移动滑动窗来动态地选定样本集，即落入滑动窗的k个样本构成一个样本集合，计算获得样本集合的统计特性；

步骤5：利用核密度估计方法，计算各个样本集关于步骤3中提取的敏感特征下的概率密度函数；

步骤6：计算出相同敏感特征下相邻两个样本集的概率密度函数的K-L散度值；

步骤7：计算集成K-L散度，以集成K-L散度值作为监测对象的健康评估指标，集成K-L散度组成新的时间序列用来判断监测对象健康状态退化的严重程度。集成K-L散度越小，相邻的两个样本集的概率密度函数越相似；反之，当集成K-L散度越大，相邻的两个样本集的概率密度函数差别越大；换句话说，集成K-L散度值越小表示相邻两个时刻的退化程度越不明显，而当集成化K-L散度越大表示相邻两个时刻的退化程度越明显。

进一步地，所述步骤2中得到的通过总体平均经验模态分解方法和希尔伯特变换等信号处理方法获取多个时域和频域特征。

进一步地，所述步骤3中使用主成分分析进行特征变换；具体计算过程如下：

步骤3.1：由所述步骤2中提取的时域和频域特征组成的样本观测矩阵：

$X=(x_1,x_2,...,x_p)=\left[\begin{array}{cccc}{x}_{11}& x_{12}& ...& x_{1p}\\ x_{21}& x_{22}& ...& x_{2p}\\ \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}\\ x_{n1}& x_{n2}& ...& x_{np}\end{array}\right]---\left(1\right)$

其中，n维向量x_i(i＝1,2,…p)代表某一个具体的时域或频域的统计特征，p代表特征总维数，n代表总的观测样本数；

步骤3.2：对观样本观测矩阵进行数据标准化处理，减少不同特征单位不同引起的影响，得到新的标准化处理之后的矩阵

步骤3.3：计算矩阵X^nor的协方差矩阵 $C=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}(x_i^{nor}-\mathrm{\μ}){(x_i^{nor}-\mathrm{\μ})}^T,$ 其中， $\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i^{nor};$

步骤3.4：计算协方差矩阵C的特征值以及特征向量：

Cν_i＝λ_iν_i(2)

其中，λ_i表示协方差矩阵C的第i个特征值，ν_i是λ_i对应的特征向量；

步骤3.5：按照降序对所有的特征值进行重新排序，得到新的特征值序列：Λ＝[λ₁,λ₂,…,λ_p]，其中，λ₁＞λ₂＞…＞λ_p，对应的特征向量序列为V＝[ν₁,ν₂,…,ν_p]；

步骤3.6：按照预先规定的累积贡献率选择出m个主成分：

$\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_i/\underset{i=1}{\overset{p}{\Σ}}{\mathrm{\λ}}_i\≥\mathrm{\θ}---\left(3\right)$

步骤3.7：由m个主成分组成的特征变换矩阵：

V＝(v₁,v₂,…,v_m)(4)

得到新的特征矩阵：

Y＝V^T(x₁,x₂,…,x_m)^T＝(y₁,y₂,…,y_m)^T(5)

步骤3.8：分别计算m个敏感特征的分类贡献率F_j；定义：

$F_j=\frac{{\mathrm{\λ}}_j}{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i},(j=1,2,...,m)---\left(6\right)$

构造贡献率向量F＝(F₁,F₂,…,F_m)。

进一步地，步骤3.6中累计贡献率θ取值为85％～90％。

进一步地，所述步骤4中提到的移动滑窗，窗宽为8～15。

进一步地，所述步骤6中提到的K-L散度也称为相对熵或者信息增益。在概率论和信息论中，K-L散度用来计算两个分布的分对称性或者差异性。K-L散度值越小，表明两个分布越相似。本方法中定义K-L散度为：

$D_{KL}(P,Q)=\frac{1}{2}\left(D_{KL}\right(P\left|\right|Q)+D_{KL}(Q\left|\right|P\left)\right)---\left(7\right)$

其中，P和Q分别代表两个离散的概率分布，D_KL(P||Q)表示分布P和分布Q的对数差期望，具体的表示形式为由步骤5得到各类样本集的概率密度函数之后，就可以根据式(7)计算出相邻两个样本集对应的K-L散度值。

进一步地，所述步骤7中集成K-L散度值IKL_i(i＝1,2,…,n-k)为：

${\mathrm{IKL}}_i=\underset{j=1}{\overset{m}{\Σ}}{F}_j\×{\mathrm{KL}}_i^j---\left(8\right)$

其中，F_j为步骤3得到的第j个敏感特征的贡献率，表示在第j个敏感特征描述下的第i个K-L散度值。式(8)中通过将第j个敏感特征描述下的K-L散度值与第j个敏感特征的贡献率F_j加权，得到各类样本集的集成K-L散度值。

本发明的有益效果：由于使用了核密度估计(KernelDensityEstimation,KDE)和K-L散度(Kullback-LeiblerDivergence)两种统计学工具，所以充分考虑到样本的统计不确定性，进而提高了健康状态评估模型的准确性和泛化性。

附图说明

图1为本发明的基于概率密度函数的健康状态评估方法的流程示意图；

图2为本发明实例所针对的来自美国智能维修系统中心的滚动轴承加速寿命试验平台的整体布局图；

图3为本发明的基于核密度估计和K-L散度的健康评估模型的原理图；

图4为样本集的动态选择过程示意图；

图5为第二组试验1号轴承的原始特征经过PCA特征变换后主成分的重要度谱线图；

图6为第二组试验1号轴承的健康评估结果；

图7为EMD分解和Hilbert变换得到第二组试验1号轴承第680个样本的四个IMF分量及对应的包络谱。

具体实施方式

下面将结合附图和具体的实施例对本发明作进一步的阐述。

如图1所示，其步骤包括：

步骤1：采集被检测对象的原始振动数据。

本发明以美国智能维修系统中心(IntelligentMaintenanceSystem,IMS)的滚动轴承加速寿命试验平台为例，具体的试验参数如下：

如附图2所示，该试验台在同一根转轴上安装有四个测试轴承。转轴通过一台交流电机通过带传动驱动。转轴转速保持为2000r/min。转轴径向加有2700Kg的载荷，通过弹簧机构施加。所有被测试轴承均采取强制润滑措施。在润滑供油系统的反馈油管加入一个磁性柱塞，用来收集润滑油液中的轴承脱落颗粒，当收集到的颗粒达到一个特定水平，整个试验结束。试验中使用的轴承为双排圆柱棍子轴承，型号为RexnordZA-2115。该类轴承每排有16个滚子，轴承节径为71.5mm，滚动体直径为8.4mm，接触角为15.17°。选择型号为PCB353B33的高灵敏度加速度传感器进行振动信号采集，该传感器布置于轴承座上。试验过程中，使用美国NI公司的数据采集仪器每隔10min进行一次数据采集，采样频率恒定为20kHz，采样时长为1s，即一次采样得到20480个数据点。通过三次全寿命试验得到对应的三组试验数据。第一组试验数据中，3号轴承出现了内圈故障，而四号轴承出现了滚动体故障；第二组试验结束时，1号轴承出现了外圈故障；第三组试验结束时，3号轴承出现了外圈故障。

为了验证本发明所提方法的可行性和有效性，但考虑到第一、三组试验的数据采集时间过长，所以我们在本发明中只使用了第二组试验数据进行方法验证。第二组试验结束时，观察到1号轴承发生了外圈故障，此时总共收集到984组振动数据。轴承的具体参数如下表所示：

表1试验轴承的参数表

由以上轴承的参数值可以得到试验中所用到的滚动轴承外圈故障特征频率(BallPassFrequencyatOuterRace,BPFO)f_BPFO：

步骤2：利用EEMD和希尔伯特变换等信号处理方法从原始振动数据中提取制定时域频域特征。

利用EEMD对各个数据集对应的振动信号样本进行分解，白噪声初始值为0.3，集成数为100，于是可以得到每个样本的前四阶固有模态分量(IntrinsicModeFunction,IMF)。接下来，分别从各个数据集对应的原始振动信号和前四阶IMF中提取指定的9个时域特征和10个频域特征，从而得到每类样本集所对应的95个统计特征。所提取的时域特征和频域特征种类如表2和表3所示。

表2时域特征表

表3频域特征表

步骤3：利用PCA特征变换方法，把从步骤2中选择出10个敏感特征，并且计算这些敏感特征的贡献率，以供后续步骤使用。应用PCA的特征变换方法得到第二组试验1号轴承的原始特征PCA变换后主成分的重要度谱线图如图3所示。由图3可以看出前10个主成分的贡献率已经可以达到90％，因此，前10个主成分被选作敏感特征，即达到了特征降维的目的。

步骤4：使用窗宽k＝10的移动滑窗来动态选定样本集，即使落入滑动窗的10个监测样本构成一个样本集合。如图3所示，1-10个样本被选进第一个样本集，接着滑动窗往后滑动一个，第1个样本出窗，第11个样本进窗，即第2-11个样本构成新的样本集，直到遍历完所有样本。更直观地，如图4中所示例子，窗宽为10的绿色样本框选中第一个样本集，一般作为正常样本集；而紫色框选中的样本集即为退化过程样本集，紫色框不断向后推移；蓝色框选中的样本集是故障样本集。从绿色框到紫色框再到蓝色框，每个框所包含的样本集的统计特性都是逐渐变化的，差异不断积累，宏观地表现就是发生故障。本发明中使用的试验数据长度是984，最终可以动态选择出975个样本集。

步骤5：利用核密度估计计算每一个样本集关于步骤3中提取的敏感特征的概率密度函数。

一种估计统计数据概率密度函数的非参数方法，设x₁,x₂,…,x_k是取自一元连续总体X的样本，在任意点x处的总体密度函数f(x)的核密度估计定义为：

$\widehat{f_h}\left(x\right)=\frac{1}{kh}\underset{i=1}{\overset{k}{\Σ}}K\left(\frac{x-x_i}{h}\right)---\left(7\right)$

其中，K(·)称为核函数，h为窗宽。上式的窗宽h是核密度估计方法中唯一需要优化的参数，根据最小化集成平方误差的方法，针对任何一个样本集，样本集的大小是10，带入下面的式子有：

$h={\left(\frac{4{\widehat{\mathrm{\σ}}}^5}{3\×10}\right)}^{\frac{1}{5}}\≈1.06\widehat{\mathrm{\σ}}\×{10}^{-\frac{1}{5}}$

则可以计算出核密度估计的最佳带宽h。

基于上述方法，可得到975个样本集的关于第j个特征的概率密度函数KD_i ^j(j＝1,2,…,10；i＝1,2,…,975)。

步骤6：在选定的特征描述下，求出各个样本集的概率密度函数。

由步骤5得到各个样本集的概率密度函数，根据对称化K-L散度定义，有：

$D_{KL}(P,Q)=\frac{1}{2}\left(D_{KL}\right(P\left|\right|Q)+D_{KL}(Q\left|\right|P\left)\right)$

求得在步骤3中选择的10个敏感特征描述下，相邻两个样本集的概率密度函数的K-L散度值KL_i ^j(j＝1,2,…,10；i＝1,2,…,974)。

步骤7：计算集成K-L散度，并通过集成K-L散度的大小来评估监测对象的退化严重程度。

由步骤6可知，我们通过计算得到了10个敏感特征描述下的K-L散度值，根据式(11)中集成K-L散度值IKL_i的定义，有：

${\mathrm{IKL}}_i=\underset{j=1}{\overset{10}{\Σ}}{F}_j\×{\mathrm{KL}}_i^j(i=1,2,...,974)$

式中，F_j为步骤3得到的10个敏感特征中的第i个敏感特征的分类贡献率。

集成K-L散度的大小反映了监测对象的退化严重程度，进而来评估监测对象的健康状况。本发明所用到的第二组试验1号轴承的健康评估曲线如图6所示。由图6中可以看出，第二组试验的1号轴承的健康状态明显出现了两个阶段，分别是轻度损伤状态和严重损伤状态。在健康状态曲线取得第一个峰值之后，健康状态有好转趋势，这是因为轴承进入了新的非健康稳定状态，直到稳定被打破，损伤快速加剧。从轻度损伤状态到严重损伤状态，大约经历200个样本(70小时)，维修时间同样充足。从健康状态曲线图上可以看出第680个样本有明显上升，我们有理由认为监测对象此时已经发生轻度损伤，而本发明提出的这一时刻要比同类健康评估方法的预测时间提前三个多小时。

为了验证第680个样本是否真的出现故障，将第680个样本提出，并利用经验模式分解的方法对其进行4层分解，进一步得到4层IMF对应的Hilbert包络谱，如图7所示。图7(b)第一个IMF对应的包络谱中可以清晰地发现236.3Hz的频率成分及其2～4倍频成分，而236.3Hz这一频率成分与之前由公式(12)计算的该滚动轴承的外圈故障特征频率极其接近，这一点准确地证实了所监测轴承发生了外圈故障。综上所述，本发明提出的健康状态评估方法能更早的判断出监测对象的故障发生时刻并且有很高的准确性。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (7)

1.本发明属于机械装备的性能退化评估领域，具体涉及一种基于核密度估计和K-L散度等统计学方法以实现旋转机械的健康状态评估，包括如下步骤：

步骤1：采集被监测对象的原始振动数据；

步骤2：从步骤1中得到的原始振动数据中提取多个的时域和频域特征，获得时域特征集和频域特征集；

步骤3：对步骤2中得到的时域和频域特征集进行降维处理获得敏感特征，并且计算这些敏感特征的贡献率，然后获得贡献率向量；

步骤4：采集被监测对象每一个时刻的原始振动数据为一个样本，使用一个窗宽为k的移动滑动窗来动态地选定样本集，即落入滑动窗的k个样本构成一个样本集合，计算获得样本集合的统计特性；

步骤5：利用核密度估计方法，计算各个样本集关于步骤3中提取的敏感特征下的概率密度函数；

步骤6：计算出相同敏感特征下相邻两个样本集的概率密度函数的K-L散度值；

步骤7：计算集成K-L散度，以集成K-L散度值作为监测对象的健康评估指标，集成K-L散度组成新的时间序列用来判断监测对象健康状态退化的严重程度；集成K-L散度越小，相邻的两个样本集的概率密度函数越相似；反之，当集成K-L散度越大，相邻的两个样本集的概率密度函数差别越大；换句话说，集成K-L散度值越小表示相邻两个时刻的退化程度越不明显，而当集成化K-L散度越大表示相邻两个时刻的退化程度越明显。

2.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤2中得到的通过总体平均经验模态分解方法和希尔伯特变换等信号处理方法获取多个时域和频域特征。

3.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤3中使用主成分分析进行特征变换；具体计算过程如下：

第一步：由所述步骤2中提取的时域和频域特征组成的样本观测矩阵：

其中，n维向量x_i(i＝1,2,…p)代表某一个具体的时域或频域的统计特征，p代表特征总维数，n代表总的观测样本数。

第二步：对观样本观测矩阵进行数据标准化处理，减少不同特征单位不同引起的影响，得到新的标准化处理之后的矩阵

第三步：计算矩阵X^nor的协方差矩阵其中，

第四步：计算协方差矩阵C的特征值以及特征向量：

Cν_i＝λ_iν_i(2)

其中，λ_i表示协方差矩阵C的第i个特征值，ν_i是λ_i对应的特征向量。

第五步：按照降序对所有的特征值进行重新排序，得到新的特征值序列：Λ＝[λ₁,λ₂,…,λ_p]，其中，λ₁＞λ₂＞…＞λ_p，对应的特征向量序列为V＝[ν₁,ν₂,…,ν_p]。

第六步：按照预先规定的累积贡献率选择出m个主成分：

其中，累计贡献率θ一般取值85％～90％。本方法中θ取值为90％，表示特征变换以后，新特征代表了原始数据总方差90％的比例。

第七步：由m个主成分组成的特征变换矩阵：

V＝(v₁,v₂,…,v_m)(4)

得到新的特征矩阵：

Y＝V^T(x₁,x₂,…,x_m)^T＝(y₁,y₂,…,y_m)^T(5)

第八步：分别计算m个敏感特征的分类贡献率F_j；定义：

构造贡献率向量F＝(F₁,F₂,…,F_m)。

4.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤4中提到的移动滑窗，窗宽分别选取8～15中的每一个数，通过结果对比发现窗宽选定为10最合适。

5.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤5中提到的核密度估计是概率论中一种估计统计数据概率密度函数的非参数方法，设x₁,x₂,…,x_k是取自一元连续总体X的样本，在任意点x处的总体密度函数f(x)的核密度估计定义为：

其中，K(·)称为核函数，h为窗宽。上式的窗宽h是核密度估计方法中唯一需要优化的参数，根据最小化集成平方误差的方法，得到下式：

其中，是样本标准差，k为各类样本集所包含的样本数。

6.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤6中提到的K-L散度也称为相对熵或者信息增益。在概率论和信息论中，K-L散度用来计算两个分布的分对称性或者差异性。K-L散度值越小，表明两个分布越相似。本方法中定义K-L散度为：

其中，P和Q分别代表两个离散的概率分布，D_KL(P‖Q)表示分布P和分布Q的对数差期望，具体的表示形式为由步骤5得到各类样本集的概率密度函数之后，就可以根据式(7)计算出相邻两个样本集对应的K-L散度值。

7.如权利要求1所述的基于核密度估计和K-L散度的旋转机械健康评估方法，其特征在于：所述步骤7中集成K-L散度值IKL_i(i＝1,2,…,n-k)为：

其中，F_j为步骤3得到的第j个敏感特征的贡献率，表示在第j个敏感特征描述下的第i个K-L散度值。式(8)中通过将第j个敏感特征描述下的K-L散度值与第j个敏感特征的贡献率F_j加权，得到各类样本集的集成K-L散度值。