基于威布尔分布的婴儿培养箱关键部件的剩余寿命预测

2021-10-16缪吉昌纪晓宏王兆源陈宏文张斌

中国医学物理学杂志 2021年9期

缪吉昌，纪晓宏，王兆源，陈宏文，张斌

1.南方医科大学南方医院设备器材科，广东广州510515；2.南方医科大学南方医院医疗质量管理科，广东广州510515；3.南方医科大学南方医院总务处，广东广州510515

前言

近年来，国家愈发重视医疗质量。医疗设备的质量控制标准逐步完善，这对医疗设备的可靠性和安全性提出了严苛的要求。医疗设备在使用过程中，由于疲劳、磨损、腐蚀、冲击等原因，不可避免的会存在性能退化，并最终演化成故障。风机是婴儿培养箱中的关键部件，风机的性能决定婴儿培养箱中噪声、温度、湿度等质量控制指标是否满足要求。因此，为了保证婴儿培养箱风机的可靠性，对其维护策略的研究具有重要意义。

为了做好婴儿培养箱的维护，提高婴儿培养箱的运行可靠性。文献［1-5］采用响应式维护（Reactive maintenance）或称为事后维护策略，在设备故障后进行维修，是一种传统的被动维护方式。文献［6］采用预防性维护（Preventive maintenance）策略，针对传统模式不能满足高可靠性设备的可靠性要求，转而频繁更换可以继续使用的零部件，这种维护方法基于时间周期，对高可靠性的设备进行过度维护，维护成本高昂。文献［7-8］采用状态监测维护（Condition-based maintenance）策略，根据从设备采集的物理量，发现早期的故障现象从而避免严重的故障发生，该方式不单纯依赖时间，但是状态信号采集的难度大。预测性维护（Predictive maintenance）是一种精准、最优化的维护策略［9-11］，通过预测和避免非预期的停机事件，将“计划外停机”转变为越来越少的“计划停机”，从而减少过度维护，最大化零部件使用寿命，降低成本，提高设备的可靠性，保障医疗设备的诊疗质量。在医疗设备维护领域，预测性维护鲜有应用，主要原因是用于预测性维护策略的数据获取难度大。

本文采用基于失效数据统计分布拟合的方法，对婴儿培养箱风机进行剩余寿命预测，从而指导进行预防性维护。首先，对获取的婴儿培养箱风机的失效数据进行预处理；其次，运用先验知识对婴儿培养箱风机的失效数据进行建模；再次，对建立好的模型进行参数估计；最后，对婴儿培养箱的风机进行剩余寿命预测。

1 方法

在本节中，首先简要介绍威布尔（Weibull）分布，然后对获取的数据进行预处理，建立可靠性模型，并对模型的参数进行估计，最后进行剩余寿命预测。

1.1 威布尔分布

威布尔［12］分布的概率密度函数的表达式为：

其中，t是自由变量，η>0 是尺度参数，β>0 是形状参数，γ≥0是位置参数。

威布尔分布函数的表达式为:

当位置参数γ=0 时，三参数威布尔分布转换为含η,β的二参数威布尔分布。

在威布尔概率密度函数表达式中，参数的不同可以表示不同的分布类型，如指数分布、瑞丽分布、正态分布等。尺度参数η决定其尺度变化［13］。图1a为不同η值时概率密度函数f(t)的曲线。从图中可以看出，概率密度曲线形状相似，曲线在坐标轴上的尺度（比例）不同。形状参数β决定概率密度函数的曲线呈单峰形状或者单调下降形状［13］。图1b为不同β值时概率密度函数f(t)的曲线。从图中可以看出，当β=1 时，曲线近似指数分布；当β=2 时，曲线近似瑞丽分布；当β∈[ ]3,4 时，曲线近似正态分布。威布尔分布的形状参数β表征函数曲线的形状，威布尔分布的尺度参数η表征函数曲线在横坐标轴上的伸缩。位置参数γ在不改变概率密度函数曲线形状的情况下，小幅度修正曲线，提升拟合效果。

图1 威布尔分布函数概率密度曲线Fig.1 Probability density curve of Weibull distribution

1.2 数据预处理

产品丧失其规定的功能称为失效［14］。从产品工作开始到产品失效的时间称为产品的失效时间或者使用寿命［15］。本研究的数据是婴儿培养箱风机的失效时间，来源于医院及某婴儿培养箱厂家的售后工程师。由于目前婴儿培养箱风机的故障记录主要靠维修工单，存在产品启用时间不准确的问题；此外，在统计过程中存在笔误及重复记录的问题。为了提高数据统计的准确性，需要对原始故障记录从准确性和有效性两方面进行预处理，完成数据审核、筛选工作。除了常规审核之外，还需要识别并剔除系统性故障（或称为非随机性故障）［16］，如控制主板故障引发的风机失效。然后将筛选的失效时间数据从小到大排序，按时间分段，统计在每一时间段内的频数，作出相应的直方图，并计算失效数据的描述性统计量。根据失效数据的直方图形状，运用先验知识初步判断其是否服从威布尔分布。

1.3 可靠性模型建立

可靠性分析是剩余寿命预测的基础［10］。可靠性是指设备在规定条件和规定时间间隔内执行指定功能的能力［17］。它的特征量有可靠度、失效率等。从工程实践看，风机的轴磨损是风机失效的主要原因。在风机的耗损失效期［18］，机械磨损导致性能下降，预期的剩余使用寿命减少。本研究中，假设1：风机的性能退化从最佳状态开始，直到发生故障，期间不存在维修；假设2：风机性能退化至失效的时间为随机变量；假设3：风机的性能退化量随时间的增长单调递增。

威布尔模型所回归出的性能曲线数学规律较为明显，对未来性能曲线的走势拟合更为方便［19］。威布尔分布既能描述失效率上升的产品寿命，又能描述失效率下降的产品寿命［15］。本研究中，产品的使用寿命时间是指风机从正式使用开始到风机失效的时间长度，产品的失效时间是随机的，则风机在t时刻失效的概率可表示为：

其中，F(t)为失效（寿命）分布函数，表示在t时刻前发生失效的概率。f(t)为F(t)的概率密度函数。在工程实践中，参数通过试验评估确定。

可靠度函数，表示产品的正常使用时间长于t的概率，用R(t)表示为：

R(t)与F(t)为互补函数，即R(t)+F(t)= 1。在风机工作的初始状态，R( 0 )= 1，F( 0 )= 0。随着使用时间的增长，失效数增加，可靠性降低。

失效（或称故障）率，指产品在t时刻正常，在下一单位时间内发生故障的概率。失效率函数为：

那么，可靠度函数和失效率函数的关系为：

风机的平均无故障时间为随机变量T的数学期望：

其中，Г是伽马（Gamma）函数。平均无故障时间表示产品平均工作时间的可靠性指标。人们更关注可靠度下降到一定程度时的时间，即可靠寿命。当产品的可靠度下降到r时，产品的寿命为tr，则R(tr)=r。称t0.5为产品的中位寿命，te-1为产品的特征寿命。

1.4 参数估计

文献［20］对威布尔参数估计问题进行分析，比较了图形法、矩估计、极大似然估计、最小二乘估计以及贝叶斯估计的优缺点。结果表明，极大似然法的风险最小，精度最高。本研究采用极大似然法进行参数估计［21］，首先构造似然函数，然后对似然函数求极值，可求出函数中的参数。

设总体分布的密度函数为f(t;η,β)，其中(η,β)为待估参数。获取的婴儿培养箱风机失效时间组(t1,t2,…,tn)的概率为取其概率的最大值，即可以求得参数的估计值。此处，构建似然函数为：

对上式两边取对数，化简得：

1.5 剩余寿命预测

设T是风机的失效时间，假设风机已经运行时间t，那么剩余使用寿命时间为条件随机变量Tt=T-t，其中T>t。记剩余使用寿命Tt的分布函数为Ft(x)，则：

由于可靠度函数和分布函数为互补函数，则风机的剩余使用寿命Tt的分布函数可转化为：

风机的剩余寿命Tt的可靠度函数为：

风机的剩余寿命Tt的概率密度函数为：

2 结果

风机失效时间的原始数据在时间轴上的分布，如图2a 所示，由图可知，风机的失效时间t分布在区间［7 473,15 766］，表示婴儿培养箱风机失效主要发生在耗损失效期。将风机的失效时间数据按从小到大排序，然后将失效时间数据点投影到时间轴。从时间轴对数据进行分段，统计每一个失效时间段内的失效数据的频次，结果如图2b所示，本文中的分段数为16。

图2 数据预处理结果Fig.2 Data preprocessing results

用二参数威布尔概率密度函数对图2b频次图进行拟合，构建极大似然函数，利用Matlab（Version 9.7,The MathWorks,Inc.,美国）求解极大似然函数，求得形状参数β=8，尺度参数η=12 817，则婴儿培养箱风机使用寿命的概率密度函数为：

风机使用寿命的分布函数为：

图3a 是采用威布尔概率密度函数对失效数据频次拟合的结果，可以看出，威布尔分布能较好地拟合风机的失效数据。通过计算，风机失效数据的均值为12 044 h，标准差为1 831 h，表示风机在运行到12044 h 时，故障高发，应予以关注。图3b 是风机的故障率函数曲线，故障率随时间增长不断提高，这与工程实际相符。图3c 是风机运行的可靠度函数，随着风机运行时间的增长，风机的可靠度逐渐降低。图3d是风机运行的瞬时失效率函数。失效率函数单调递增，表示风机失效发生在耗损失效期，这间接证明了数据来源的真实性。

图3 风机失效数据建模结果Fig.3 Modeling results of fan failure data

2.1 模型评估

Q-Q 图用于比较已知样本的分布和猜测分布的接近程度［22］。Q-Q图是一种散点图，横坐标为某一样本的分位数，纵坐标为另一样本的分位数，横坐标与纵坐标组成的散点图代表同一个累积概率所对应的分位数。如果散点分布在直线y=x附近，则这两个样本是同等分布。假设风机的寿命分布与威布尔分布非常接近，则由分位数构成的点应大致成一条直线。图4是风机失效数据与威布尔分布的Q-Q 图，由图可知风机的失效数据大部分落在y=x直线上，表示风机的失效数据近似为上述威布尔分布。而个别数据偏离y=x直线，可能是风机的故障发生在随机失效期［18,23］，也有可能是数据统计过程中的错误记录。这就要求医疗机构应重视医疗设备的状态监测，重视设备运行过程中的数据收集和分析。

图4 失效数据-威布尔分布的Q-Q图Fig.4 Q-Q graph of failure data-Weibull distribution

2.2 预测

上述结果表明，婴儿培养箱风机的寿命分布符合威布尔分布，则风机的剩余寿命概率分布函数为：

对于待进行预测性维护的婴儿培养箱A、B、C、D，假设设备A 的风机运行时长为10 213 h，设备B 的风机运行时长为12 044 h，设备C 的风机运行时长为13 875 h，设备D 的风机运行时长为15 706 h，可以得出A、B、C、D 4台设备的风机的剩余寿命可靠度。图5a为4台设备剩余寿命的可靠度曲线。由图可知，运行时长最长的D 设备，其可靠度曲线斜率变化最快；运行时长最短的A 设备，其可靠度曲线斜率变化最慢。随着风机的运行时长增加，同一水平可靠度的剩余寿命逐步减少。

图5b 为4 台设备剩余寿命的概率密度曲线。由图可知，运行时长最短的A 设备，其剩余寿命的概率密度曲线趋于平坦，对于估计剩余寿命置信区间的范围宽；随着风机的运行时长增加，其概率密度曲线的最大值点左移，曲线越来越陡峭。运行时长最长的D 设备，概率密度曲线非常陡峭，对于估计剩余寿命置信区间的范围非常窄。这种形态的变化，体现为概率密度函数的方差逐渐变小，表明剩余寿命预测的结果趋于准确。在待预测的4 台设备中，D（15 706 h）设备对应的剩余寿命的概率密度函数最窄，预测的精度最高，其中位寿命为253 h。A（10 213 h）设备对应的剩余寿命的概率密度函数最宽，预测的精度最低，其中位寿命为2 357 h。在对婴儿培养箱进行预测性维护时，优先对设备D进行维护。