石 慧， 康 辉， 任谦力， 曾建潮， 谷丰收

(1.太原科技大学 电子信息工程学院，太原 030024； 2.中北大学 大数据与视觉计算研究所，太原 030051； 3.哈德斯菲尔德大学 计算与工程学院， 英国 哈德斯菲尔德 HD1 3DH)

近年来，随着信息和传感技术的高速发展，设备健康监测变得越来越智能化。以机械设备中最重要的部件之一齿轮箱为例，其健康状况直接影响到整个设备的安全可靠性，当齿轮箱中的部件发生故障时，会加速其他部件的退化，导致机械设备损坏，不仅会产生巨大的经济损失，还会引发一系列安全隐患[1-2]。因此有必要针对系统的关键部件进行实时监测并预测其剩余寿命，以便对部件进行及时的更换或维修。

系统的剩余寿命预测是故障预测与健康管理(prognostics and health management，PHM)的核心。现阶段预测剩余寿命的方法有4种：基于知识的预测方法、基于物理模型的方法、基于数据驱动的方法以及混合预测的方法。基于知识的预测方法通常需要具备有关系统的专门知识[3]，同时在复杂机械设备中，故障的物理机理很难获得，从而导致基于物理模型的方法适用性不强[4-5]。而基于数据驱动的方法不依赖物理模型，利用接收到的数据进行预测建模，具有较高的预测精度和较强的泛化能力[6-7]。现有的一些基于统计及概率知识的数据驱动剩余寿命预测方法需假定退化过程并进行参数估计，然而假设的模型可能与实际的模型不符，并且参数寻优过程容易陷入局部最优。核密度估计方法不附加任何假设，是一种非参数估计方法，它从数据本身出发来研究样本的分布特征。同时与深度学习等方法相比，核密度估计方法能够清楚地表征系统退化特征的变化。因此，在预测与健康管理领域具有较高的关注度[8]。

Hu等[9]提出非参数核密度估计方法来估计风速的概率密度函数，验证了该方法的准确性和灵活性。Wang等[10]将核密度估计方法引入鲁棒随机配置网络，以减少噪声数据或离群点对最终学习模型的影响。Zhang等[11]提出基于核密度估计的非线性系统异常检测方法来模糊正常样品之间的差异，强化潜在异常样品的异常。Sidibé等[12]利用核密度估计系统可靠性函数，研究了核平滑参数的可变性对成本模型的影响。虽然核密度估计方法已应用于系统的故障预测与健康管理，但大多数文献都是针对系统中的单个部件进行可靠性建模研究。

现代的设备系统越来越复杂，系统中有许多重要的子系统和部件协同工作，部件之间退化过程中存在随机相关性影响[13]。例如风电机组齿轮箱的齿轮和轴承退化相互影响。构建多部件系统的PHM系统就要充分考虑部件间的随机相关性[14]。

前期研究多部件系统随机相关性的文献主要集中于由部件故障触发的相互作用。一些文献研究部件随机相关性时仅考虑部件故障会对其他部件的退化产生影响，建立基于时间的维修决策模型[15-16]。而在实际运行过程中，获取故障率需要大量破坏性试验，研究具有随机相关的多部件故障预测与健康管理问题更关注的是系统故障前各部件的连续随机退化特征与实时剩余寿命的关联性。Rasmekomen等[17]指出部件间的随机相关性不仅可以由部件故障触发，还可能由部件的退化触发，并将此定义为退化相互作用，使用回归模型来表征部件退化状态对其他部件退化速率的影响。Bian等研究一个部件的退化程度如何影响系统中其他部件的退化率，并将这种由随机相关引起的部件间相互作用称为退化率相互作用(degradation-rate-interaction，DRI)，建立了离散DRI模型来描述系统中相互依赖的部件退化过程之间的相互作用。Do等[18]针对具有随机相关性和经济相关性的两部件系统，提出了基于状态的维修模型。Shahraki等[19]针对由随机相关部件组成的复杂系统进行系统建模并制定选择性维修策略。Xi等[20]将多部件系统中隐藏的退化相关性视为噪声影响建立状态空间模型，提出了在线剩余寿命预测方法。上述文献大多假设部件退化模型，并进行参数估计。宋仁旺等[21]利用Copula函数对齿轮箱振动和噪声之间的相关性进行建模，求得齿轮箱的剩余寿命。Li等[22]研究了考虑随机退化的多部件系统的视情维修问题，采用Lévy Copula函数对系统的随机相关性进行建模。杨志远等[23]基于伽马过程建立系统退化模型，使用Copula函数描述退化过程的相关性，进行系统可靠度分析和剩余寿命预测。Copula函数可以灵活构造多变量分布，但通过Copula函数构建多部件系统中部件间的联合分布函数，不能明显地反映部件间的随机相关性特征，不同Copula函数的选择也会对剩余寿命的预测结果产生较大偏差。

本文在对部件间的退化随机相关性特征及对退化状态的影响进行系统分析的基础上，研究多部件系统基于非参数核密度估计的实时剩余寿命预测方法。首先针对多部件系统中的复杂双向随机相关性进行建模；引入窗宽因子构建基于自适应窗宽的核密度估计退化模型以解决高密度区域过度平滑或低密度区域噪声干扰的问题；然后通过首达时间的概念，求解部件的剩余寿命概率密度函数；最后在斜齿轮箱试验台上进行试验，利用调制信号双谱分析对传感器数据进行特征提取，在识别非线性影响的同时考虑幅值和相位的影响，从而有效地抑制随机噪声的干扰，验证了所提模型的准确性和适用性。

1 核密度估计模型的构建

1.1 随机相关性分析

从多部件系统部件间存在的复杂随机相关性的特征以及对部件连续退化状态的不同影响角度出发，可将部件间的随机相关性分为三类：单一单向随机相关性、单一多向随机相关性和双向随机相关性。随机相关性示意图，如图1所示。

(a) 单一单向随机相关性

(b) 单一多向随机相关性

(c) 双向随机相关性图1 随机相关性示意图Fig.1 Schematic diagram of stochastic dependence

(1) 单一单向随机相关性：是指某个部件的退化只会对系统中的单个部件的退化产生单向影响；部件a自身的退化只对部件b的退化产生单向影响，部件a不受其他任何部件退化的影响。

(2) 单一多向随机相关性：是指某个部件的退化会对系统中多个部件的退化产生单向影响或者是某个部件的退化会受到系统中其他多个部件退化的单向影响；部件a的退化会对部件b、部件c等多个部件的退化产生影响，其自身退化不受其他任何部件退化的影响。

(3) 双向随机相关性：是指部件的退化会受到系统中其他部件的影响，同时自身的退化也会对其他部件的退化产生影响。部件a自身的退化会对部件b和部件c等多个部件的退化产生影响，同时，这些受影响的部件退化到一定程度也会对部件a的退化产生影响，这种影响是双向存在的。

本文考虑更为复杂的部件间具有双向随机相关性的多部件系统进行研究。假设一个系统中的部件i和部件j具有双向随机相关性，即一个部件的退化会加剧另一个部件的退化，反之亦然。随着现代传感器技术的发展，可以通过传感器监测系统中各部件的历史退化状态数据及实时运行数据，利用这些样本数据进行随机相关性建模并建立剩余寿命预测模型。

1.2 双向随机相关性的模型构建

(1)

(2)

(3)

(4)

(5)

(6)

将式(5)代入式(6)，可以求出

(7)

式中，σk-1为k-1个已知样本的标准差。

(8)

(9)

(10)

1.3 自适应窗宽下核密度估计退化分布模型的实时更新

(11)

部件i在tk时刻基于自适应窗宽的核密度估计为

(12)

2 部件退化分布的计算

(13)

(14)

tk-1+m时刻获得m个样本数据后，概率密度函数为

(15)

同理，在任意时刻获得新的样本数据后，都可获得相应的累积退化概率密度函数。

3 剩余寿命预测模型的建立

(16)

由式(15)可知

(17)

将式(17)代入式(16)，可得

(18)

则tk-1时刻部件i的剩余寿命概率密度函数为

(19)

4 试验验证与分析

试验中将监测到的部件振动信号数据利用调制信号双谱(modulation signal bispectrum，MSB)变换进行特征提取，得到可以表征部件退化过程的退化状态曲线，同时确定部件的阈值。将特征提取后的退化数据代入文中的剩余寿命预测模型，便可通过式(16)～式(19)求得部件的剩余寿命概率密度函数，求得部件的剩余寿命。剩余寿命预测试验验证流程如图2所示。

图2 剩余寿命预测试验流程图Fig.2 Experimental flow chart of remaining useful life prediction

4.1 数据采集

本文采用英国哈德斯菲尔德大学效率与性能工程中心的斜齿轮箱试验台进行试验，试验台的基本结构，如图3所示。两个齿轮箱(GB1和GB2)采用背靠背的结构进行安装，它们分别与驱动电机和负载电机相连。在GB1和GB2的箱体表面分别安装传感器，用于采集信号。采用16通道数据采集系统(data acquisition system，DAS)采集振动信号、声学信号和油温信号，电机额定转速为1 500 r/min，振动信号的采样频率为96 kHz。

图3 斜齿轮箱试验台Fig.3 Helical gearbox test bench

本文使用的两个斜齿轮箱的几何参数,如表1所示。

表1 齿轮箱规格Tab.1 Gearbox specification

图4显示了安装在两个齿轮箱中的齿轮分布情况，试验采用齿轮正反交错搭接的啮合方式。

图4 齿轮箱示意图Fig.4 Schematic diagram of gearboxes

轴频和啮合频率可以通过表2所示的公式计算。

表2 轴频和啮合频率Tab.2 Shaft frequency and mesh frequency

为模拟工业设备在实际运行过程中环境和负载的随机变化情况，在试验期间，首先控制驱动电机以正弦曲线变化的转速运转30 min，然后在5种不同负载条件下进行匀速运转(全速的50%)，最后在同样的5种负载下以不同的恒定速度(全速的70%)运行。详细的测试条件,如图5所示。

图5 循环运行条件Fig.5 Cycle operating conditions

试验持续了838 h，监测到GB2的振动出现明显峰值后停止试验。打开齿轮箱后，观察到GB2在低速阶段发生磨损，齿根部位磨损较为严重，如图6所示。

(a)

(b)图6 838 h的齿轮磨损情况Fig.6 Gear wear of 838 hours

4.2 数据处理

试验运行了838 h，齿轮在运行初期样本特征波动较大，为齿轮啮合阶段，使用这个阶段的数据来分析部件的磨损退化过程是不恰当的，因此选择运行300 h之后的数据来分析，对振动数据进行MSB变换，进行特征提取并预测其剩余寿命。

传统的频域特征提取方法有频谱分析、功率谱分析、常规双谱分析等，频谱和功率谱分析在对信号经过傅里叶变换后仅对幅值信息进行分析，忽略了相位信息。而常规的双谱B(fx,fc)通过傅里叶变换后在频域内可表示为

B(fx,fc)=E[X(fc)X(fx)X*(fc+fx)]

(20)

式中：X(f)为信号x(t)的离散傅立叶变换；fx为调制频率；fc为载波频率；X*为X的复共轭。

对于振动信号x(t)其MSB变换在频域中表示为

BMS(fx,fc)=

E[X(fc+fx)X(fc-fx)X*(fc)X*(fc)]

(21)

MSB对常规双谱改进后在进行退化特征提取时，同时考虑信号幅值和相位的变化，能够更好地说明调制信号和载波信号之间的非线性关系，充分表示信号的调制特性，要比常规的双谱分析方法更准确[26-27]。

以幅值和相位来表示，可将式(21)写为

BMS(fx,fc)=E〈|X(fc+fx)|X(fc-fx)||X*(fc)|·

|X*(fc)|exp[jφMS(fx,fc)]〉

(22)

式(22)中MSB的总相位可通过以下方式计算

φMS(fx,fc)=φ(fc+fx)+φ(fc-fx)-

φ(fc)-φ(fc)

(23)

φ(fc)和φ(fx)耦合时，相位关系可以表示为

φ(fc+fx)=φ(fc)+φ(fx)

(24)

φ(fc-fx)=φ(fc)-φ(fx)

(25)

将式(24)和式(25)代入式(23)，可以得到MSB的总相位为零，它的幅值由其4个分量幅值的乘积所确定，因此，如果(fc+fx)和(fc-fx)来源于fx和fc的非线性影响，在双频BMS(fx,fc)处会出现明显的双谱峰。这样对MSB的表述将更为准确。与此相反，如果fx和fc无非线性作用，或者它们是随机噪声，其MSB的相位随时间变化，最终平均结果为零，由此MSB具有很好地抑制噪声作用。图7给出了采用MSB对本文试验台数据特征提取的调制双谱分析图。

(a) MSB幅值谱

(b) MSB相干谱图7 MSB双谱分析Fig.7 MSB bispectrum analysis

如图7所示，通过MSB双谱分析先从幅值谱发现具有峰值的信号幅值，同时，该坐标点对应的相位谱中也有较大的相干系数，则说明该幅值可以反映部件的退化特征，若相位谱中对应着非常小的相干系数，说明该幅值不是通过调制作用生成，可能是由噪声产生的，从而将其剔除。其在识别非线性影响的同时考虑幅值和相位的影响，从而有效地抑制随机噪声的干扰。提取能够表征部件退化的振动信号特征数据出来，便可得到部件的退化状态曲线图。

图4中GB2中的齿轮5和GB1中的齿轮3通过主轴进行连接，二者之间具有较强的相关性，故将齿轮5和齿轮3作为部件i和部件j进行研究。采用MSB特征提取方法从MSB频谱中将能够表征部件退化的特征幅值提取出来，便可得到部件i的退化状态曲线如图8所示，GB2在838 h监测到磨损严重，发生故障。其故障阈值为0.229 mm/s2。

图8 部件i的退化状态Fig.8 Degradation state of component i

4.3 结果分析

为了说明选择不同的窗宽对核密度估计结果的影响，图9给出了部件i在450 h自适应窗宽与固定窗宽下的核密度估计结果。

(a) 样本分布情况

(b) 核密度估计结果图9 两种窗宽下核密度估计结果Fig.9 Results of kernel density estimation under two bandwidths

由图9可知，基于自适应窗宽的核密度估计对样本的拟合程度更高，与固定窗宽下的核密度估计结果相比，更能够反映样本的分布情况。

将MSB特征提取后的退化数据代入模型，可根据式(16)～式(19)求得部件i在当前监测时刻的剩余寿命概率密度函数。图10给出本次试验在7个监测时间下的剩余寿命预测结果。可以看出，随着监测时间的增加，样本数据增大，本文提出方法的剩余寿命预测结果越来越接近真实值。说明本文提出的基于自适应核窗宽的核密度估计剩余寿命预测方法可以很好的对部件i剩余寿命的概率密度函数进行估计。

图10 核密度估计方法在不同监测时间下剩余寿命预测结果Fig.10 Remaining useful life prediction results of kernel density estimation method under different monitoring time

Do等考虑具有随机相关性和经济相关性的两部件系统，使用参数估计进行建模计算。为了对本文所用的非参数方法进行进一步验证，使用相同的样本数据，基于参数估计的方法进行剩余寿命预测，与本文方法进行比较。

基于参数估计的剩余寿命预测方法结果，如图11所示。图11中：虚线为基于参数估计方法计算的剩余寿命预测结果；实线代表本文利用自适应核窗宽核密度估计方法得到的剩余寿命预测结果；“*”号为部件i的真实剩余寿命结果。从图11可以看出，随着监测时间的增加，样本数据量增大，采用非参数核密度估计方法进行剩余寿命预测更加接近真实的剩余寿命，说明本文提出的非参数核密度估计方法与参数估计相比，有一定的优越性。

图11 参数估计与核密度估计方法比较结果Fig.11 Comparison between parameter estimation and kernel density estimation

表3给出了参数估计和核密度估计方法对部件剩余寿命预测结果的误差分析。从表3可以看出，随着监测时间的推移，获得的样本数据越多，两种方法对部件剩余寿命的预测就越精确。在同一监测时间下，本文提出的核密度估计的方法要比参数估计方法更精确，进一步验证了非参数核密度估计方法的有效性和精确性。

表3 参数估计和核密度估计剩余寿命预测误差分析Tab.3 Error analysis of residual life prediction of parameter estimation and kernel estimation

图12用三维网格曲面图和二维等高线图绘制了部件i的平均剩余寿命和部件i与部件j退化的关系图，以此来说明部件j对部件i的随机相关性影响。图12(a)中，在部件i自身退化量一定的情况下，部件i的平均剩余寿命随着部件j退化量的增加而减少；在部件j退化量一定的情况下，部件i的剩余寿命随着部件i自身退化量的增加而减少。图12(b)中，从浅色区域到深色区域，部件i的平均剩余寿命逐渐减小，同样反映出部件i的平均剩余寿命不仅与自身退化有关，还会受到部件j退化量的影响。说明部件j对部件i具有随机相关性影响。

(a) 三维网格图

(b) 二维等高线图图12 部件j对部件i的随机相关性影响Fig.12 The stochastic dependence effect of component j on component i

同理，图13用三维网格曲面图和二维等高线图来说明部件i对部件j的随机相关性影响。图13(a)中，在部件j自身退化量一定的情况下，部件j的平均剩余寿命随着部件i退化量的增加而减少；在部件i退化量一定的情况下，部件j的剩余寿命随着部件j自身退化量的增加而减少。图13(b)中，从浅色区域到深色区域，部件j的平均剩余寿命逐渐减小，同样反映出部件j的平均剩余寿命不仅与自身退化有关，还会受到部件i退化量的影响。说明部件i对部件j具有随机相关性影响。

(a) 三维网格图

(b) 二维等高线图图13 部件i对部件j的随机相关性影响Fig.13 The stochastic dependence effect of component i on component j

综上可知，部件i和部件j具有双向随机相关性影响，进一步说明了本文考虑多部件系统部件间的双向随机相关进行建模的必要性。

为了更好的说明所提模型的通用性，采用Xi等考虑多部件随机相关性试验所用的数据，再次对文中模型进行验证分析。

部件i的退化状态,如图14所示。

图14 部件i的退化状态Fig.14 Degradation state of component i

将退化数据代入本文核密度估计模型，并与参数估计模型进行比较，得到两种方法下部件i在不同监测时刻的剩余寿命预测结果,如图15所示。

图15 参数估计与核密度估计方法比较结果Fig.15 Comparison between parameter estimation and kernel density estimation

从图15可以看出，随着监测时间的增加，两种方法的预测结果越来越接近真实值，而本文提出的核密度方法比参数估计的方法更加准确，再一次证明了文中提出的自适应核密度估计模型的适用性。

5 结 论

针对多部件系统部件连续退化过程中存在的双向随机相关性，基于自适应窗宽的核密度估计方法建立了实时剩余寿命预测模型。首先针对部件间存在的双向随机相关性影响，引入条件核密度估计来进行建模，然后采用自适应窗宽的核密度估计方法求得相应的概率密度函数，最后建立实时剩余寿命预测模型，得到部件的剩余寿命。利用斜齿轮箱试验台进行验证分析，结果表明该方法可以很好的进行寿命预测，通过和基于参数估计的剩余寿命预测方法进行比较，进一步验证了该方法的准确性和有效性。