1
引言
当前, 工业生产体系伴随着数字化技术与移动互联网的蓬勃发展, 掀起了一场万物互联智慧化的新兴革命, 信息技术正在与工业生产中的基础设施和管理系统相融合, 以将传统工业体系提升到更高水平. 德国“工业4.0”[1 ] 、“中国制造2025”[2 ] , 标志着工业生产制造从自动化时代全面转向信息化与智慧化时代, “智慧”体现了信息化所采取的方法和工具[3 ] . 物联网 (Internet of Things, IoT) 是工业信息化发展的基础技术, 依托无线网络、移动设备、SOC、传感器等多种技术的进步, 在集成度、灵敏性以及成本控制等方面愈发成熟[4 ] . 侯一鸣等研究了基于物联网和工业云的选矿设备状态监控系统[5 ] , 利用多传感器融合技术实现了视频数据与设备过程参数融合监控. 相互连接的设备定期收集、分析数据, 借助云与AI技术, 为复杂重型装备的健康监测和工业生产的规划、管理、决策提供了智慧辅助支持.
而设备健康状态的准确评估对于复杂重型装备有着重要意义. 目前, 设备故障状态评估主要根据设备的监测数据和实际生产系统的信息, 得出设备的健康状况, 确定设备是否继续工作或进行维护. 评估方法大多集中在系统建模与数据驱动建模两个方面. 系统建模主要有故障树分析法、神经网络分析、隐马尔可夫模型等[6 ] . 随着工业物联网与大数据分析的发展, 结合人工智能技术, 相继出现了支持向量机[7 ] 、神经网络[8 ] 、深度学习[9 ] 等健康评估方法, 这些方法可以较好地检测设备故障, 也存在着一些缺陷: ① 各部件传感器收集的特征参数较多, 直接输入到评估模型中会增加计算复杂度, 导致评估结果不准确; ② 大多是针对已存在故障的状态分类识别, 不能直观量化当前设备偏离正常状态的程度, 不能全面描述设备的退化过程[10 ] . 为了解决该问题, 马田系统 (Mahalanobis-Taguchi System, MTS)被应用于滚动轴承初始故障检测和状态监测方法研究, 可以准确地检测出轴承的初始故障和判断轴承的退化状态, 但由于故障的复杂性和多变性, 该方法还需进一步改进. 针对上述问题, 本文提出了一种结合马田系统与设备状态指数(Device Status Index, DSI)的状态评估方法. 通过筛选多个设备特征参数, 构建马田系统的基准空间, 利用故障敏感度计算特征参数对应的加权马氏距离, 之后结合设备状态指数, 构建评估模型, 利用设备状态指数的变化来判断设备的健康状态, 最后基于轴承标准数据集进行模拟分析, 验证方法的有效性.
2
设备健康状态评估基本流程
2.1
设备故障监测问题分析
设备从购买、安装、开始工作到故障报废的整个生命周期中, 其内在部件的状态会随着时间推移发生一系列变化. 因此, 设备故障的变化趋势是从轻微故障征兆开始, 逐渐发展到整个设备丧失工作能力的过程, 如图1 所示. 其中, 设备一开始的状态处于稳定阶段, A点是即将出现轻微故障的征兆点; B点表示出现故障征兆点, 在A点与B点之间的阶段, 设备处于亚稳定阶段, 可以通过监测数据发现异常, 不影响设备正常工作; C点表示设备出现明显故障, 该故障可以直接通过观察或外在特征体现出来.
1
设备生命周期曲线
在进行设备健康状况评估时, 需要提取已采集的特征数据中与设备退化相关的特征数据集, 构建评估矩阵. 假设一台设备有m 个特征值与设备退化相关, 在单位时间段内共有n 个数据, 可定义评估矩阵:
\begin{document}${{{{E}}}} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{{{x}}_{11}}}&{{{{x}}_{12}}} \end{array}} \\ {\begin{array}{*{20}{c}} {{{{x}}_{21}}}&{{{{x}}_{22}}} \end{array}} \\ {\begin{array}{*{20}{c}} {...}&{...} \end{array}} \\ {\begin{array}{*{20}{c}} {{{{x}}_{n1}}}&{{{{x}}_{n2}}} \end{array}} \end{array}}&{\begin{array}{*{20}{c}} {...} \\ {...} \\ {...} \\ {...} \end{array}}&{\begin{array}{*{20}{c}} {{{{x}}_{1 m}}} \\ {{{{x}}_{2 m}}} \\ {....} \\ {{{{x}}_{nm}}} \end{array}} \end{array}} \right]_{n\times m}}$\end{document}
其中, x ij 表示设备第i 次采集时第j 个特征值的数据.
距离是体现不同样本之间差异性的常用工具, 其值的大小与样本间的差异性成正比. 因此, 可以通过合适的距离度量设备正常状态与故障状态的差异度. 马氏距离 (Mahalanobis Distance, MD ) 是由印度统计学家马哈拉诺比斯 (P. C. Mahalanobis)提出, 用于表示数据的协方差距离, 有效计算两个样本集之间的相似度[11 ] , 有
\begin{document}$\mu =({\mu }_{1}, {\mu }_{2}, \mathrm{\cdots}, {\mu }_{{p}}{)}^{{\rm{T}}}$\end{document}
, 协方差矩阵为Z 的多变量向量
\begin{document}${{ x = }}{({x_{\rm{1}}},{x_{\rm{2}}},\cdots,{x_{{p}}})^{\rm{T}}}$\end{document}
, 其马氏距离定义为:
1
\begin{document}$ {MD}(x) = \sqrt {{{(x - \mu )}^{\rm T}}{{Z}}_{}^{ - 1} {x - \mu } } $ \end{document}
相比于欧式距离, 马氏距离在考虑特征值的前提下, 排除了特征相关性的干扰, 在多维特征尺度下可以较好量化设备状态. 从图1 可知, 运行时间增加, 设备运行性能会逐渐下降, 而用于量化设备退化状态的MD 会加速增长. 可以通过映射函数阈值直观地评价设备的健康状态, 将MD 映射到指定范围内. 映射过程可表述如式(2):
2
\begin{document}$ \left\{ {\begin{array}{*{20}{l}} {F:[0,x) \to [a,b],x \in (0, + \infty )} \\ {DSI = F(M{D_t}) \in [0,1]} \end{array}} \right. $ \end{document}
其中, x 可表达为特征对应的MD 范围, 映射函数F 的函数值称为设备状态指数.
2.2
评估基本流程
从设备运行的历史数据中提取多种特征参数构建稳定基准空间, 将特征按照故障敏感度进行筛选, 求得对应的加权马氏距离. 然后利用Box-Cox变换获得设备状态指数阈值, 构建健康状态模型, 实现评估分析. 如图2 所示.
2
设备状态评估流程
根据图2 可知评估基本流程主要有以下4个过程:
(1) 构建特征集. 对设备运行状态初始数据进行预处理, 提取其中相关特征参数, 构建初始特征集
\begin{document}${X}= $\end{document}
\begin{document}$ ({{x}}_{ij}{)}_{n\times m}$\end{document}
;
(2) 计算加权马氏距离. 将X 作为输入, 通过特征参数的故障敏感性筛选中评估重要特征, 并计算出相应的加权马氏距离(Weighted Mahalanobis Distance, WMD);
(3) 构建设备状态指数模型. 根据WMD构建DSI模型, 利用Box-Cox变换和3
\begin{document}$\sigma $\end{document}
准则确定DSI模型中的阈值;
(4) 评估设备健康状态. 根据确定的阈值与特征参数的DSI, 对设备的健康状态进行评估.
3
马田系统算法设计
田口玄一博士提出的马田系统将马氏距离与田口方法有效集成, 广泛应用于疾病诊断、数据分类、模式识别以及样本的诊断、预测分析[12 ] . 利用MTS进行设备状态识别, 需要收集正常的样本数据集, 并将其定义为基准空间; 随后以基准空间为基点, 求得样本参数对应的MD 作为度量尺度, 测度未知样本与基准空间的距离; 最后通过设定阈值进行状态识别.
3.1
基准空间的构建
通过收集设备特征数据构建正常样本集, 并计算相应MD , 详细步骤如下:
(1) 构建初始特征集: 识别设备的m 个重要特征参数
\begin{document}${{{x}}_j} \;(j = 1,2,3,\cdots,m)$\end{document}
以此构建初始特征集
\begin{document}${X}=({{x}}_{{i}}{}_{{j}}{)}_{{n\times}m}$\end{document}
.
(2) 剔除异常点, 构建稳定基准空间: 第一步构建的初始特征集中, 可能存在部分异常数据, 会导致最后计算的样本MD 不稳定. 根据Rousseeuw提出的改进最大分类器差异(Maximum Classifier Discrepancy, MCD )算法[13 ] , 该方法大致步骤为: 先找到一个样本量为h 的子集, 使得在所有大小为h 的子集中, 该子集的协方差矩阵的行列式是最小的. 根据MCD 计算协方差估计量, 计算公式如下, 获取均值和协方差估计量后, 最后可以通过计算得到每个样本与中心之间的马氏距离, 如果马氏距离大于某个临界值, 则该点视为离群点.
3
\begin{document}$ \sum {MCD = \frac{{{k_{MCD}}(h,n,p)}}{{h - 1}}} \sum\limits_{i \in {I_{MCD}}} {({x_i} - {{\hat{\mu } }_{MCD}})({x_i} - {{\hat{\mu } }_{MCD}})'} $ \end{document}
该方法是一种高鲁棒性的估计方法, 它的目标就是找出协方差矩阵具有最低行列式的观测值, 改进的MCD 算法, 其基本思想是一个包含顺序统计和行列式的不等式, 以及我们称之为 “选择性迭代” 和 “嵌套扩展” 的技术, 该方法对数据中的异常值进行有效识别, 剔除初始特征集中的异常点, 构建稳定的基准空间.
(3) 正常特征参数标准化:
4
\begin{document}$ {{{\hat{X} }}_{{{ij}}}} = \frac{{{{{x}}_{ij}} - {{\bar x}_j}}}{{{\delta _j}}} $ \end{document}
5
\begin{document}$ {{{\bar X}}_{{j}}} = \frac{{\displaystyle\sum\nolimits_{{{i}} = {\rm{1}}}^{{n}} {{{{x}}_{ij}}} }}{{{n}}},{\delta _j} = \frac{{\displaystyle\sum\nolimits_{{{i}} = {\rm{1}}}^{{n}} {{\rm{(}}{{{x}}_{ij}}} - {{{{\bar x}}}_{{j}}}{)^2}}}{{{{n - 1}}}} $ \end{document}
标准化后的稳定基准空间表示为:
6
\begin{document}$ {{{\hat{X} }}_{{{ij}}}} = {[{{{\hat{X} }}_{\rm{1}}},\;{{{\hat{X} }}_{\rm{2}}}{\rm{,}}\;{{{\hat{X} }}_{\rm{3}}},\cdots,{{{\hat{X} }}_{{n}}}]^{\rm{T}}} $ \end{document}
(4) 计算MD : 在MTS中, MD 平方后被用作标准度量尺度 (即: 将得到的MD 进行平方计算再赋值给MD ). 通过式 (6) 处理后, 正常样本的MD 期望值在1附近分布, 利于区分异常样本[14 ] . 其计算公式被定义为:
7
\begin{document}$ {{M}}{{{D}}_{{i}}} = \frac{1}{{{m}}}{{{\hat{X} }}_{{i}}}{{{S}}^{{\rm{ - }}1}}{{\hat{X} }}_{{i}}^{{\rm{T}}} $ \end{document}
8
\begin{document}$ {{S}} = \frac{{\displaystyle\sum\nolimits_{{{i}} = 1}^{{n}} {{{{{\hat{X} }}}_{{i}}}{{\hat{X} }}_{{i}}^{{\rm{T}}}} }}{{{{n - 1}}}} $ \end{document}
式 (7) 是正常样本中m 个特征的相关系数矩阵. 为避免出现式 (4) 中S 矩阵不可逆问题, 采用田口玄一博士提出的GSP方法来计算MD , 计算公式为:
9
\begin{document}$ {{M}}{{{D}}_{{i}}} = \frac{1}{{{m}}}\sum\limits_{j = 1}^m {\frac{{u_{ij}^2}}{{\zeta _j^2}}} $ \end{document}
其中,
\begin{document}${U_i} = {[{u_{i1}},{u_{i2}}, \cdots ,{u_{ij}}]^{\rm{T}}}$\end{document}
为
\begin{document}${{{\bar X}}_{{i}}}$\end{document}
通过GSP方法得到的正交向量,
\begin{document}${\zeta _{{j}}}$\end{document}
为
\begin{document}${{{U}}_{{i}}}$\end{document}
的特征标准差.
3.2
基准空间有效性验证
在基准空间优化前需要对其进行有效性验证. 使用正常样本的期望和标准差对异常样本进行标准化处理, 获取异常样本的MD . 如果正常样本的MD 小于异常样本, 证明度量尺度良好, 建立的基准空间是有效的, 反之, 则需要重新选择特征变量定义基准空间.
3.3
基准空间优化
设备状态评估的过程中, 并非所有的特征变量都有助于提高计算精度, 有些特征变量可能对最终数据没有影响, 而有些甚至存在干扰. 因此, 有必要对经验证后的基准空间进行特征优化, 选择对构建MTS基准空间有正收益的特征变量. 在MTS中, 使用正交表OA和信噪比SNR 相结合来筛选有效特征. 根据特征参数的个数设计正交表, 假设样本空间有p 个特征, 安排在正交表的前p 列上, 用两个标量分别表示该特征是否参与构建基准空间. 对于每次试验(正交表的行数), 使用被选择特征计算异常样本的马氏距离
\begin{document}${{MD}}_{{i}}, \;{i}= $\end{document}
\begin{document}$ 1,2,\mathrm{\cdots},m$\end{document}
并计算SNR :
10
\begin{document}$ {{SNR}} = {\rm{ - }}10{\rm{lg}}\left( {\frac{1}{m}\sum\limits_{i = 1}^m {\frac{1}{{M{D_i}}}} } \right) $ \end{document}
对于每一个特征, 分别用
\begin{document}${{{t}}_1}$\end{document}
和
\begin{document}${{{t}}_2}$\end{document}
表示特征
\begin{document}${{{x}}_j}$\end{document}
参与实验的SNR 均值与未参与的SNR 均值.
11
\begin{document}$ \Delta = {{{t}}_1} - {{{t}}_2} $ \end{document}
如果
\begin{document}$\Delta > 0$\end{document}
, 则表明该特征可以保留, 反之, 删除该特征.
3.4
基样本空间的加权MD
优化后的特征变量可以作为评估的重要特征参数, 之后可以根据特征的重要性赋予特征参数不同的权重, 以体现其贡献程度[15 ] . 经过GSP处理后的加权马氏距离可定义如下:
12
\begin{document}$ {{WM}}{{{D}}_{{i}}} = \frac{1}{{{m}}}\sum\limits_{j = 1}^m {\frac{{u_{ij}^2}}{{\zeta _j^2}}} {W_j} $ \end{document}
由式(12)可知权重会影响到基准空间的有效性. 特征的权重应当与该特征对故障的敏感性有关. 特征对异常样本越敏感, 标明该特征包含的变化信息越多, 更利于分类和预测. 在设备运行生命周期中, 不同的特征参数在经过A点到B点的时间间隔中, 数据表现是不同的. 根据不同特征对故障的敏感性不同, 本文通过线性变化函数将所有的正常样本归一化, 计算其对应的敏感性. 归一化计算公式如下:
13
\begin{document}$ {{\bar X}} = \frac{{{{{x}}_{{i}}} - {x_{\min }}}}{{{{{x}}_{{\rm{max}}}} - {x_{\min }}}},\; i = 1,2,3,\cdots,n $ \end{document}
上述由于故障的敏感性主要由设备运行生命周期的AB段体现, 故对特征参数
\begin{document}${{\bar X}}$\end{document}
的AB段进行提取, 获得峰值p 与时间t , 通过式(14)求得其敏感性:
14
\begin{document}$ {{S}} = \frac{{{p}}}{{{t}}} $ \end{document}
算法特征对故障的敏感性来确定的权重被定义为:
15
\begin{document}$ {{{W}}_{{k}}} = \frac{{{{{s}}_{{k}}}}}{{\displaystyle\sum\nolimits_{{{k}} = 1}^{{m}} {{{{s}}_{{k}}}} }} $ \end{document}
至此, 通过基准空间中的特征参数获得了WMD, 以此确定故障的可能性.
4
设备状态指数模型
4.1
状态指数模型构建
由基准空间求得的WMD与设备的运行状态密切相关, 其关系可通过映射函数体现. 当函数值接近范围上限b 时, 意味着设备工作正常, 生命周期处于稳定阶段; 当设备退化到一定程度时, 函数值处于(a , b )范围内, 设备工作存在隐患, 生命周期可能处于亚稳定阶段; 当函数值接近范围下限a 时, 设备故障可能性极高, 甚至停工. 加权后的马氏距离, 其期望值分布在[0, 1]附近, 利用
\begin{document}${{{\rm e}}^x}$\end{document}
数可以在保持原函数单调性的前提下增加敏感性的特点, 构造DSI函数:
16
\begin{document}$ {{DSI}} = \frac{2}{{1 + {{\rm{e}}^{\alpha \cdot WMD}}}} $ \end{document}
其中,
\begin{document}$\alpha $\end{document}
是调节参数, 保证严重故障时DSI趋于范围下限, 其值可以通过设备稳定状态的WMD期望中
\begin{document}$\overline {{{WM}}{{{D}}_{{h}}}} $\end{document}
和设备健康置信度
\begin{document}${{Be}}{{{f}}_{{h}}}$\end{document}
确定:
17
\begin{document}$ \alpha = {\rm{ - }}\dfrac{{{\rm{ln}}\left( {\dfrac{{Be{f_h}}}{{2 - Be{f_h}}}} \right)}}{{\overline {{{WM}}{{{D}}_{{h}}}} }} $ \end{document}
设备健康置信度可由同类设备的历史运行数据统计得出. 从式(17)可知, 设备健康运行, WMD在1附近波动, DSI值略低于1; 设备轻微故障时, WMD会超过某阈值, DSI会逐渐减小, 与WMD成反比; 设备严重故障时, WMD远大于1, DSI急速减小并趋于0.
4.2
DSI阈值确定
DSI阈值可以准确区分设备生命周期的正常和异常状态, 小于DSI阈值表明正常; 超过DSI阈值表明异常. 利用设备特征参数的MD 计算DSI阈值. 考虑到特征参数经过剔除异常点后仍可能存在错误, 影响判断的准确率. Kumar等[16 ] 提出一种Box-Cox变换方法, Box-Cox变换的一般形式如下.
18
\begin{document}$ y(\lambda ) = \left\{ \begin{array}{l} \dfrac{{{y^\lambda } - 1}}{\lambda },\lambda \ne 0 \\ \ln y,\lambda = 0 \\ \end{array} \right. $ \end{document}
其中,
\begin{document}$\lambda $\end{document}
是一个待定变换参数, Box-Cox变换是一族变换, 这是因为不同的
\begin{document}$\lambda $\end{document}
, 对其进行的变换也不尽相同, 它包括了平方根变换, 对数变换, 倒数变换等常用变换, 在确定变换参数后, 使得
\begin{document}${y^{(\lambda )}}$\end{document}
满足:
19
\begin{document}$ {y^{(\lambda )}} = X\beta + e,e \sim N(0,{\sigma ^2}I) $ \end{document}
即要求通过因变量的变换, 使得变换过的向量
\begin{document}${y^{(\lambda )}}$\end{document}
回归自变量具有线性相依关系误差也服从正态分布, 误差各分量是等方差且相互独立, 故Box-Cox变换是通过参数的适当选择. 达到对原来数据的“综合治理”, 使其满足一个正态线性回归模型的所有假设条件.
故而可以将不服从正态分布的WMD值转化为近似正态分布的数据, 具体转换参照式(20).
20
\begin{document}$ {{M}}{{{D}}_{{i}}}(\gamma ) = \left\{ {\begin{array}{l} {\dfrac{1}{\gamma }({{WMD}}_{{i}}^\gamma - 1),\gamma \ne 0;}\\ {{\rm{ln}}({{WMD}}_{{i}}^{}),\gamma = 0.} \end{array}} \right. $ \end{document}
其中,
\begin{document}${{WM}}{{{D}}_{{i}}}$\end{document}
第i 个样本的加权马氏距离,
\begin{document}$ {{MD}}_{{i}}(\gamma )$\end{document}
为相应的变换值,
\begin{document}$\gamma $\end{document}
值由最大似然求得:
21
\begin{document}$\begin{split} {{f}}(MD,\gamma ) =& - \frac{n}{2}\ln \left[ {\sum\nolimits_{i = 1}^n {\frac{{{{(M{D_i}(\gamma ) - \overline {MD} (\gamma ))}^2}}}{n}} } \right]\\ &+ (\gamma - 1)\sum\nolimits_{i = 1}^n {\ln (M{D_i})} \end{split} $\end{document}
其中,
22
\begin{document}$ \overline{MD}(\gamma )=\frac{1}{{n}}{\displaystyle \sum _{{i}=1}^{{n}}{{MD}}_{{i}}(\gamma )}$\end{document}
然后基于3
\begin{document}$\sigma $\end{document}
准则 (拉伊达准则) 确定DSI 阈值, 满足正态分布或近似正态分布的样本数据处理, 如表1 所示.
1
3
\begin{document}$\sigma $\end{document}
数值分布
数值分布
占比
\begin{document}$ (\mu \rm{-}\sigma , \mu +\sigma )$\end{document}
0.6827
\begin{document}$ (\mu \rm{-2}\sigma , \mu +2\sigma )$\end{document}
0.9545
\begin{document}$ (\mu \rm{-3}\sigma , \mu +3\sigma )$\end{document}
0.9973
利用3
\begin{document}$\sigma$\end{document}
准则对经过变换后的
\begin{document}$ {{MD}}_{{i}}(\gamma )$\end{document}
确定阈值
\begin{document}${{{T}}_{2\sigma }}$\end{document}
与
\begin{document}${{{T}}_{3\sigma }}$\end{document}
:
23
\begin{document}$ \left\{ {\begin{array}{*{20}{c}} {{{{T}}_{2\sigma }} = \mu MD(\gamma ) + 2\sigma MD(\gamma )} \\ {{{{T}}_{3\sigma }} = \mu MD(\gamma ) + 3\sigma MD(\gamma )} \end{array}} \right. $ \end{document}
根据确定的
\begin{document}${{{T}}_{2\sigma }}$\end{document}
、
\begin{document}${{{T}}_{3\sigma }}$\end{document}
与
\begin{document}$\gamma $\end{document}
得到变换前正常特征样本的WMD阈值, 通过式 (16) 计算DSI的阈值,
\begin{document}$2\sigma $\end{document}
阈值对应的是设备生命周期的故障征兆点,
\begin{document}$3\sigma $\end{document}
阈值对应故障点.
24
\begin{document}$ \left\{ {\begin{array}{c}{{DSI}}_{{\text{征兆}}}=\dfrac{2}{1+{\rm{e}}^{\alpha \cdot M{D}_{2\sigma }}}\\ {{DSI}}_{{\text{故障}}}=\dfrac{2}{1+{\rm{e}}^{\alpha \cdot M{D}_{3\sigma }}}\end{array}} \right.$\end{document}
利用此方法确定评估模型的阈值, 在不同情况下DSI值映射范围内的置信度达到99%以上, 说明构建的设备健康状态识别模型具备良好的准确率.
5
实验分析
滚动轴承是复杂重型装备的关键部件, 其健康状态关乎复杂重型装备运行情况. 本文选择采用滚动轴承作为实验研究对象来验证方法有效性, 其结果对复杂重型装备具有较大的适用性. 使用由西安交大与昇阳科技联合实验室发布轴承加速寿命实验数据XJTU-SY Bearing Dataset[17 ] , 实验平台如图3 所示. 数据采集由固定在测试轴承的水平和竖直方向上的两个单向加速度传感器获得. 采样频率为25.6 kHz, 采样间隔为1 min, 每次采样时长为1.28 s, 每份样本包含32 768个数据点, 实验选用1号工况数据集.
3
轴承加速退化测试平台
(1) 数据预处理
振动信号在采集的过程中会因为外界诸多因素的干扰导致其内部存在噪声, 这些噪声会使得真实数据出现非平滑, 非线性等特点, 从而导致后续分析存在误差甚至错误[18 ] . 因此, 信号在进行分析前进行预处理是非常必要的. 通过对Storm中的模型Bolt进行自定义, 可以在信号数据被后续数据分析服务消费前进行预处理. 本文采用小波变换的方法[19 ] 对振动信号进行预处理.
25
\begin{document}$ WT(a,\tau ) = \frac{1}{{\sqrt a }}\smallint _{ - \infty }^\infty f(t)\times\varphi \left( {\frac{{t - \tau }}{a}} \right){{d}}t$\end{document}
式(25)为信号
\begin{document}${{f}}(t)$\end{document}
的连续小波变换公式,
\begin{document}${{a}}$\end{document}
控制其伸缩量,
\begin{document}$\tau $\end{document}
控制其平移量. 通过改变其值, 达到调整频窗的效果, 从而自适应的表达信号中的低频与高频. 当信号被小波函数分解后, 对分解出的参数进行阈值处理, 最后将处理后的参数进行逆变换, 将信号重构, 保证去噪后仍旧保留信号特征. 含有噪声的以为信号模型表达式可定义为:
26
\begin{document}${{F}}(n{\rm{)}} = f(n) + \varepsilon \times \zeta (n),\;n \in [0,1,2,\cdots,k - 1]$\end{document}
其中,
\begin{document}${{F}}(n)$\end{document}
为包含噪声信号,
\begin{document}$f(n)$\end{document}
为原始信号,
\begin{document}$\zeta (n)$\end{document}
为噪声系数的标准偏差. 通常信号的有用部分表现在低频部分或表达平稳, 因此, 可以按照如图4 所示进行分解, 一般情况下, 噪声位于
\begin{document}${{{b}}_1}$\end{document}
,
\begin{document}${{{b}}_2}$\end{document}
,
\begin{document}${{{b}}_3}$\end{document}
处. 只需对其进行小波处理, 将得到的信号重构即可得到目标结果. 对于振动数据, 选用db4小波函数分解4层后, 效果如图5 所示, 变换后的信号更平滑.
4
信号分解流程
(2) 特征选择
在滚动轴承的振动信号分析中, 常在时域与频域上选取特征参数. 设共有m 份样本数据, 每份数据的长度为n , 第i 份样本数据的第j 个数据点用
\begin{document}${{{Z}}_{{{ij}}}}$\end{document}
表示. 本文选用的特征参数表如表2 所示[20 ] . 其中, 频域特征参数中的k , d , D , a 分别代表轴承中滚珠个数、滚珠直径、轴承中径与接触角, 其幅值根据样本时域信号的快速傅立叶变换计算得到. 本文选择数据集中一组轴承从完好到失效的全生命周期的数据构建初始特征集, 共123个样本.
由于实验数据是轴承的全生命周期数据, 因此选取样本数据的前50个作为正常样本, 后50个样本作为异常样本, 构建训练集. 首先对初始特征集中的正常样本数据进行异常点筛选, 过滤异常点后将剩余的样本数据构建稳健基准空间. 之后分别计算两份样本的WMD, 通过比较发现所建立的空间是有效的, 之后通过正交表与SNR 进行特征筛选, 最终筛选出
\begin{document}${{{P}}_3}$\end{document}
、
\begin{document}${{{P}}_4}$\end{document}
、
\begin{document}${{{P}}_5}$\end{document}
、
\begin{document}${{{P}}_7}$\end{document}
、
\begin{document}${{{P}}_{12}}$\end{document}
、
\begin{document}${{{P}}_{13}}$\end{document}
、
\begin{document}${{{P}}_{14}}$\end{document}
、
\begin{document}${{{P}}_{15}}$\end{document}
构建优化后的基准空间.
(3) 结果分析
通过对正常样本的 WMD值进行Box-Cox变换, 得到近似正态分布, 从而得到最优变换参数
\begin{document}$\gamma = 0.1876$\end{document}
. 利用式(23)计算得出MD 的征兆阈值与故障阈值, 其中,
\begin{document}$ {{MD}}_{{\text{征兆}}}=2.2175$\end{document}
,
\begin{document}$ {{MD}}_{{\text{故障}}}=5.2366$\end{document}
, 正常样本与异常样本的WMD分布如图6 所示.
图6 中, 红色线为故障阈值点, 绿色线为征兆阈值点, 由图6 可知, 正常样本的WMD 值均小于故障阈值, 将近98.6%的样本值小于征兆预警阈值. 同时, 异常样本的WMD值有98.6%高于征兆阈值, 说明本文方法确定的两个阈值有较高的置信度. 因此, 可以以此构建相应的DSI模型. 以98.6%的置信度作为标准, 计算得到
\begin{document}${{\textit{DSI}}}_{{\text{征兆}}}$\end{document}
与
\begin{document}${{\textit{DSI}}}_{{\text{故障}}}$\end{document}
, 状态指数如图7 所示.
5
处理前后的振动信号对比
从图7 可以看到, 随着时间增加, 在80 min附近状态指数开始超过故障阈值, 由此说明阈值模型可以匹配实际信号, 体现了设备状态的变化趋势. 此外, 通过阈值可以对设备的运行状态提供标准, 准确识别异常状态, 判断设备何时出现故障征兆, 为设备的维修与管理提供数据支持.
2
特征参数表
域
特征参数
公式
时域
均值
\begin{document}${P_1} = \dfrac{1}{n}\displaystyle\sum\limits_{j = 1}^n { { {\textit{z} }_{ij} } }$\end{document}
标准差
\begin{document}${P_2} = \sqrt {\dfrac{1}{n}\displaystyle\sum\limits_{j = 1}^n { { {({{z_{ij} } } - {P_1})}^2} } }$\end{document}
整流平均值
\begin{document}${P_3} = \dfrac{1}{n}\displaystyle\sum\limits_{j = 1}^n {\left| {{{\textit{z}}_{ij}}} \right|} $\end{document}
均方根
\begin{document}${P_4} = \sqrt {\dfrac{1}{n}\displaystyle\sum\limits_{j = 1}^n {{{\textit{z}}_{ij}}^2} } $\end{document}
峰值
\begin{document}${P_{\rm{5}}} = \max \left| {{{\textit{z}}_{ij}}} \right|$\end{document}
偏度
\begin{document}${P_6} = \dfrac{{\displaystyle\sum\nolimits_{j = 1}^n {{{({{\textit{z}}_{ij}} - {P_1})}^3}} }}{{(n - 1){P_2}^3}}$\end{document}
峰度
\begin{document}${P_{\rm{7}}} = \dfrac{{\displaystyle\sum\nolimits_{j = 1}^n {{{({{\textit{z}}_{ij}} - {P_1})}^4}} }}{{(n - 1){P_2}^4}}$\end{document}
峰度因子
\begin{document}${P_8} = \dfrac{{{P_5}}}{{{P_4}}}$\end{document}
脉冲因子
\begin{document}$P{}_{\rm{9}} = \dfrac{{{P_{\rm{5}}}}}{{{P_{\rm{3}}}}}$\end{document}
波形因子
\begin{document}${P_{10}} = \dfrac{{{P_4}}}{{{P_3}}}$\end{document}
裕度因子
\begin{document}${P_{11}} = \dfrac{{{P_5}}}{{\dfrac{1}{n}\displaystyle\sum\nolimits_{j = 1}^n {{{\textit{z}}_{ij}}^2} }}$\end{document}
频域
外圈特征
\begin{document}${P_{12}} = \dfrac{{ - k{f_i}\left( {1 - \dfrac{d}{D}\cos \alpha } \right)}}{2}$\end{document}
内圈特征
\begin{document}${P_{13}} = \dfrac{{ - k{f_i}\left( {1 + \dfrac{d}{D}\cos \alpha } \right)}}{2}$\end{document}
滚动体特征
\begin{document}${P_{14}} = - \dfrac{{{f_i}}}{2}\dfrac{D}{d}\left( {1 - {{\left( {\dfrac{d}{D}\cos \alpha } \right)}^2}} \right)$\end{document}
保持架特征
\begin{document}${P_{15}} = \dfrac{{{f_i}\left( {1 - \dfrac{d}{D}\cos \alpha } \right)}}{2}$\end{document}