基于混合Gamma分布的机载产品可靠性建模

李军亮, 祝华远, 王 正, 王利明, 张鑫磊

(海军航空大学青岛校区, 山东 青岛 266041)

0 引 言

舰载机在其寿命周期内,要遂行多种作战任务,需在沿海、近海、远海等环境中服役,机载设备受到温度、湿度、气压、盐雾、振动、加速度等多种自然环境应力和平台环境应力耦合作用,故障率高且失效模式复杂,严重影响舰载机作战效能的发挥。因此建立科学的评价方法,并准确评价舰载机机载产品的可靠性是亟待解决的问题。

机载产品主要包括机械类、电子类以及机电耦合类产品,具有小批量、多批次、高可靠、长寿命等特点,研制生产阶段的可靠性试验只能收集在特定应力下的产品可靠性信息,不能完全模拟真实使用环境,且GJB-899A主要对寿命分布为指数分布的设备级产品的可靠性试验进行规范,未能考虑非指数型分布或具有性能退化特征的产品特征,适用范围有限[1]。工程中处理寿命分布的方法有3种：一是根据先验信息(主/客观信息)选择适用的分布;二是根据经典分布模型进行寿命数据拟合、拟合优度检验以及参数估计[2-3];三是采用贝叶斯统计方法综合上述两种方法。常用的寿命分布模型有指数、正态、对数正态、威布尔和Gamma分布等[4],但是不同的分布模型都有各自的使用范围,工程经验表明，电子产品的失效多为指数分布[5],机械部件的疲劳过程则服从威布尔分布或对数正态分布[6-7]。机载产品结构复杂,一般是机械、电子和液压等器件的组合,在实际使用过程中故障行为复杂,往往存在故障传播、故障耦合以及共因失效等特征[8-9],因此单一分布的拟合能力有限。Gamma分布既可以拟合产品的退化数据,也可以拟合其寿命数据,黄卓等[10]证明了有限Gamma分布的稠密性,并采用混合Gamma分布构建了一种通用的寿命数据拟合方法[11],基于最大期望(expectation maximun， EM)算法求解分布参数,并设计了混合分布的分支控制策略,其拟合结果与单一分布拟合相比更为精确。在此,拟根据Gamma分布的自身特性,构建一种更为简洁和有效的混合分布模型和优化求解算法,并结合实际算例对构建方法进行验证。

1 三维混合Gamma 分布的寿命分布拟合

有限混合 Gamma 分布在[0, +∞]上的全体概率分布函数中稠密[10],即可以逼近任意寿命分布函数。假设任意系统故障时间分布函数可由多个Gamma分布构成,不同产品的故障密度概率可表示为

(1)

分析式(1)可知,要想明确系统的寿命分布曲线,首先要确定式(1)中的各项参数集{al,αl,λl}以及M,文献[11]在其推理过程中并未明确给出M的确定方法,在算例分析中假设M≤5,和分布参数一起求解,求解过程比较复杂。分析Gamma分布的特点可知,当形状参数αl∈(0,+∞)时,其故障率特点可以分为3种趋势:单调增加、单调递减和恒定不变。当形状参数0<α<1,α=1和α>1时(λ=1),分别代表故障率递减、恒定或者递增,如图1所示。

图1 不同参数的Gamma分布形状 Fig1 Shape of Gamma distribution with different parameters

又因为有限维混合Gamma分布在实数概率空间中稠密[10-11],理论上不同形状的分布参数可以通过智能优化算法求解,从而获得最优解。因此,拟对M的混合策略进行简化,假设任意寿命分布曲线可由3条Gamma曲线组成,即maxM≤3,然后根据优化模型的优化目标精度控制拟合精度,从而控制各支分布的参数,其概率密度函数为

(2)

其故障分布函数为

(3)

定理1由3支Gamma 分布组成的混合Gamma 分布在[0, +∞) 上的全体概率分布函数中稠密。

证明假设g(x)为任意的概率密度函数,其数学期望为E(X)=θ,其方差为D(X)=σ2,欲证明f(x)在g(x)中稠密,只需证明存在度量函数ρ和ε→0,使得ρ≤ε,则原假设成立。

两边取平方可得

证毕

根据定理1可知,由3支Gamma分布组成的混合分布在[0, +∞)的故障分布函数中稠密,即式(3)可以拟合出任意故障分布函数,且参数集中的未知参数集为{al,αl,λl},而文献[11]的算例中则含有15个未知参数,因此本文的混合模型未知参数更少,求解过程相对简单。特别的,当al=1时,则混合分布为Gamma分布[13-14]。

2 三维混合Gamma分布优化模型构建

定理1虽然证明了混合Gamma分布的稠密性,但是在实际应用中分布参数集都是未知的,只能通过产品在服役环境下产生的故障样本得出其样本分布的一些统计特性的

2.1 小样本数据的优化模型

当样本量n≤20时,一般采用海森公式、数学期望公式或者近似中位秩公式来直接代替经验分布函数值,在此选用中位秩公式来替代经验分布,如下所示:

(4)

式中:n为样本容量;i为样本失效时间由小到大排列的顺序,即t1≤t2≤…≤tn。

采用均方误差MSE作为度量测度,即

(5)

当样本已知时,式(5)中含有未知参数集{al,αl,λl},优化模型如下:

(6)

(7)

当样本已知时,可以确定分布参数置信上下限。

2.2 大样本数据的优化模型

Person检验是一种比较通用的拟合优度检验方法,Person统计量适用于一般分布,而不针对特定分布。当样本量较大时,将Person统计量作为拟合模型和真实分布量的测度。

外场数据采集的时效性,不可能像在实验室环境下一样精确,尤其是整个机群中同类产品数据的采集可能会存在时间误差和随机误差,为了减少这种误差的影响,可以根据样本情况对采集的时间区间进行划分。在此,根据样本量n将样本按时间取值范围分为k个时间区间,则ni为第i个区间内的失效样本数,总体X落入第i个区间的概率pi为

pi=F(t(i))-F(t(i-1)),i=1,2,…,k

(8)

则Person检验统计量可表示为

(9)

将式(3)代入式(9)可得

χ2=

(10)

构建以下优化模型:

minχ2=

(11)

(12)

3 基于改进粒子群算法的优化算法设计

单一分布模型中比较常见的分布参数评估方法有最小二乘法、极大似然估计(maximum likelihood estimate, MLE)、EM算法等。文献[11]在EM算法的基础上设计了控制分支M的一种迭代算法,用来求解有限混合Gamma分布的参数求解,其求解方法和过程都比较复杂。近年来,遗传算法、粒子群算法等群智能算法被用来求解可靠性参数优化问题,并且取得了不错的效果[14-15]。其中,粒子群算法收敛速度快,且模型构建过程比较简单,在此采用自适应粒子群算法对混合分布的参数进行求解,算法的基本流程设计如图2所示。

图2 基于自适应粒子群优化算法的参数求解流程Fig.2 Parameter solving process based on adaptive particle swarm optimization algorithm

如图2所示,在使用粒子群算法的时候,首先要确定参数集P(N,c1,c2,xLIMT,vLIMT,wLIMT,M,D),其中N为群粒子数目,c1和c2为加速系数,D为种群维度即目标函数中自变量个数,M为迭代次数,xLIMT和vLIMT分别为根据约束条件设定的粒子自变量取值范围和粒子迭代移动幅度参数,vLIMT为加权系数阈值。则在第t时刻,各个粒子位置参数为Xi(t)={xi,1(t),xi,2(t),…,xi,n(t)},速度公式为Vi(t)={vi,1(t),vi,2(t),…,vi,n(t)},个体最优解为pBesti(t)={pBesti,1(t),pBesti,2(t),…,pBesti,n(t)}，全局最优解为gBesti(t)={gBesti,1(t),gBesti,2(t),…,gBesti,n(t)},位置和速度的更新方程公式分别为

vi,j(t+1)=wvi,j(t)+c1r1[pi,j-xi,j(t)]+c2r2[pg,j-xi,j(t)]

(13)

xi, j(t+1)=xi, j(t)+vi, j(t+1),j=1,2,…,d

(14)

式中:r1和r2分别为随机数介于0～1。加权系数的大小会影响算法的收敛性能,自适应权重可以有效避免陷入局部最优解,从而求得全局最优解,自适应权重更新公式为

(15)

式中:fmin和favg分别为最小和平均目标值；wmin和wmax为最小和平均权重。

4 算例分析

不同类型的机载产品的失效机理不同,在不同环境表现出的故障行为也不一致。某舰载直升机机群首次大修间隔内故障统计数据表明,机载电子产品和机电产品的故障率相对较高,本节分别以某型舰载直升机多功能显示器和微动开关的故障样本为例来验证算法的正确性。

4.1 算例1

多功能显示器故障主要是因为其内部的印制电路板(printed circuit board, PCB)失效引起,在其服役过程中受到随机振动、酸性气体、盐雾以及温湿度的作用,该设备的密封结构、防腐涂层失效进一步引起PCB上的叉指电极、通孔等发生腐蚀等故障行为。其样本故障如表1所示,表中第5列Fn(ti)为采用近似中位秩公式计算的寿命分布近似值。采用本文设计方法,将计算结果Fn(ti)代入式(6)并展开可得优化目标函数。

表1 某设备PCB板故障数据样本

(16)

分别采用经典粒子群算法和自适应粒子群算法进行求解,其中粒子群算法参数为算法的初始参数P(50,3,3,0.8,100,9),自适应粒子群算法的参数设置为P(50,3,3,0.8,0.6,100,9),其计算结果如图3所示。如图3所示,两种算法收敛速度较快,且适应度函数均接近0.003 771 096 562 428,其中自适应用粒子群算法的参数集的优化结果为(0.9, 12.19, 221.12, 0.09, 1, 221.12, 0.01, 0.2, 221.12),将该参数代入式(3)可得产品的寿命分布函数,如下所示:

F(t)=0.9Γ(t,12.194 1,221.12)+0.09Γ(t,1,1.03)+ 0.01Γ(t,20,0.99)

(17)

图3 不同算法的收敛性及最小适应值Fig.3 Convergence of different algorithms and fitness value

为了更直观地显示寿命分布的拟合曲线,分别将混合分布模型与近似中位秩法以及残存比率法进行比较,结果如图4所示,系统的可靠度如图5所示。如图4所示,近似中位法计算的故障率最高,残存比例最低,混合分布的故障函数值介于近似中位秩和残存比例法中间,因为近似中位秩法不考虑删失样本数,所以数值偏大,而残存比率法需要考虑删失样本,所以数值偏小。实际上,近似中位秩法更适用于完全样本的故障率近似,但是在实际使用过程中很难获得产品的完全样本,尤其是航空产品具有高可靠、长寿命的特征,在一定时间内很难获得产品的完全样本,其样本一般为截尾数据,残存比率法在工程中经常被用于计算具有删失特性的截尾样本数据,但其计算结果偏于保守,因此混合Gamma分布拟合方法可以有效改良两种方法的缺陷,接近产品的真实故障分布特性。

图4 3种算法的比较Fig.4 Comparison of three algorithms

图5 PCB板的可靠度函数Fig.5 Reliability function of PCB board

在装备保障的实际过程中,可结合式(17)和图5合理设定产品的可靠性阈值或者剩余寿命,从而有针对性地制定预防性维修策略。

4.2 算例2

以机载微动开关为研究对象,在舰载服役环境中,微动开关易受到酸性盐雾及温湿度的影响,开关内部触点发生腐蚀从而引起开关接触电阻变大,直至功能失效,因此其失效过程为退化失效和突发失效的综合作用。其失效前后的照片如图6所示,图6(a)为装机前的新样品,图6(b)为失效后拆卸下来的样本。

图6 DK1-2微动开关的失效前后照片Fig.6 Photos of DK1-2 microswitch before and after the failure

假设机群中每架飞机的使用环境和使用频次基本一致,采集其近5年的故障数据,得到80组故障数据,如表2所示。

表2 微动开关故障数据样本

如表2所示,样本数据中最小故障时间为218 h,最大故障时间为266 h,样本容量为80,取时间间隔Δt=10 h,可得样本的故障频率直方图如图7所示,即可获得的ni。观察图7可知，样本数据非正态、指数分布等单一分布,因此按照式(10)构建样本的优化目标函数为

minχ2=

(18)

图7 微动开关故障样本频率分布直方图Fig.7 Histogram of frequency distribution of micro switch fault samples

采用文中设计方法求解可得,参数设置为(500, 3, 3, 0.8, 0.5, 200, 9),算法收敛过程如图8所示,最佳适应值为0.001 057 3,优化结果的参数集为(0.95, 260, 1, 0.02, 1, 1.03, 0.03, 30, 0.99)。

图8 算例2的收敛过程Fig.8 Convergence process of Example 2

将所求的参数代入式(3)可得产品的寿命分布模型为

F(t)=0.95Γ(t,260,1)+0.02Γ(t,1,1.03)+ 0.03Γ(t,30,0.99)

(19)

可靠度曲线如图9所示。工程实践中,一般以威布尔分布或者对数正态分布拟合机电产品的故障曲线[13-14],图9中“△”的曲线为以算例2的样本按照99.5%置信度拟合出的服从对数正态分布的产品可靠度曲线;“*”为采用本文设计的算法计算的产品可靠度;“+”为采用文献[11]设计的算法计算的产品的可靠度。其中,“△”曲线在使用时间小于219 h之前,认为产品可靠度一直大于99%,在使用时间接近220 h产品可靠度迅速衰减,接近0,与表2中的故障样本数据反映的特征有所差别,不能有效反映出产品的真实可靠性和故障特性。本文设计算法和文献[11]中的算法计算的可靠度曲线形状比较相似,但文献[11]中的算法计算产品可靠度略低于本文算法,本文算法的可靠度值基本介于对数正态分布和文献[11]算法之间,且与表2中的故障数据比较吻合,可以较为有效表现产品的可靠性特征。

图9 微动开关可靠度Fig.9 Product reliability of micro switch

如图9所示,微动开关的可靠性在不同时间内呈现出不同的特征,在0～40 h产品可靠性较高,但处于下降阶段;在40～227 h之间可靠度趋于平稳,其最小值大于0.94;在270 h之后可靠度迅速下降;在300 h之后产品的可靠度趋近于零。因此,可结合式(19)在不同的寿命阶段合理设定产品的可靠性阈值以及预防性维修策略,将比目前定期检查策略更为有效[16-18]。

5 结 论

本文分析了在不需要借助工程经验判断产品寿命分布模型类型的前提下,如何基于故障数据构建产品的混合Gamma分布模型,并从理论上证明了该方法的可行性,设计了求解模型的粒子群算法,并通过算例进行了验证和分析。算例1结果显示,当样本量小于20时,本文方法与样本真实误差小于0.000 3,且分别与近似中位秩法和残存比例法进行了比较,证明了本文方法更为精确。算例2说明在大样本情况下,本文方法同样适用,最佳适应值为0.001 057 3,通过算法和误差控制可以获得分布参数精确数值解,拟合出合理的产品寿命分布模型。

本文方法适用于在复杂环境中服役的复杂产品,如机电组合设备,因为其失效机理较为复杂,受到多种因素的综合作用,且电子产品和机械设备的寿命分布模型往往不一致,单一分布模型在拟合其寿命分布模型的精度有限,而本文方法通过混合分布可以更加真实地反映产品的故障分布统计特征。

文中只讨论了基于故障数据即寿命数据的可靠性评估方法,对于系统的退化数据并未进行讨论,而基于退化数据和寿命数据相融合的方法更能准确的评估系统的可靠性。另外,舰载产品的服役环境非常复杂,目前尚未有有效的舰载机机载产品在服役过程中适用于环境因子折算的环境剖面划分方法以及混合分布模型的环境因子折算方法,开展此类研究将非常有意义。