基于广义线性模型的火工品可靠性数据分析

洪东跑,王英华,管 飞,温玉全

(1.中国运载火箭技术研究院,北京 100076; 2.北京理工大学爆炸科学与技术国家重点实验室,北京 100081)

1 引 言

在航天领域,火工品承担着一系列有关飞行器点火、姿态控制、分离等控制指令的执行,对飞行器飞行成败起到关键性的作用,是飞行器的关键元件。火工品可靠性要求高,需要有效的评价方法，作为高可靠性要求的敏感性产品,在国内外的工程应用中,一般利用感度试验数据来进行可靠性分析。即假设每个产品均存在一个临界刺激量,当外界施加的刺激量大于临界刺激量时,该产品“响应”; 否则“不响应”[1]。同时假定临界刺激量服从某分布,按照一定试验方法施加不同的刺激量进行感度试验,并利用试验数据进行可靠性分析[2-3]。由于火工品的临界刺激量不能直接测量,只能在若干刺激量下通过刺激试验来测定。每个火工品样本只能做一次试验,即使施加了刺激量不响应,也会破坏其结构[4]。故感度试验所获得的数据为施加的刺激量和相应的“响应”或“不响应”数,无法获得精确的临界刺激量。同时不同的方法获得的试验数据差异较大,数据分析方法也不完全一致。特别当样本量较小时,分析结果可能会有较大的偏差,影响了火工品可靠性分析的精度。而且利用感度试验数据进行参数估计时无法获得解析解,而利用数值算法求解时,其收敛性依赖于迭代初始值。由广义线性模型的相关研究[5-7]可知,该模型可有效地弥补数值算法中参数估计迭代不收敛或收敛慢的不足。

为此,本研究根据火工品可靠性数据分析的工程需求,探索了不同类型感度试验数据的共性特点和内在规律,采用广义线性模型描述了火工品可靠性与感度的关系,给出一种火工品可靠性数据分析方法,在小样本下实现对高可靠性要求的火工品可靠性分析。

2 感度试验数据特点分析

工程中常用的感度试验方法有升降法[1]、Langlie法[2]、E方法[8],D优化方法[4]等。由于升降法试验操作相对简单,是工程中最为常用的火工品感度试验方法[9]。为此,本研究以升降法试验为例,分析了火工品感度试验数据的特点。

升降法试验包括三个因素：试验样本量N、初始刺激量x0和步长d。x0和d选定后,用x0作第一次刺激-响应试验; 第二次及以后每次试验所用刺激量的取法如下：如前一次试探的反应结果为“响应”,则本次试探用刺激量为xi+1=xi-d; 如为“不响应”,则为xi+1=xi+d。如此循环试验,至完成预定试验样本量N为止。对升降法试验数据进行统计分析,按刺激量的升序排列,可表示成如下形式：

(1)

式中，k为刺激量个数,xi(i=1,2,…,k)为试验刺激量,mi为在xi试验的不响应数,ni为在xi试验的响应数。

本研究对Langlie法、E方法、D优化方法等试验数据也进行了分析,由分析结果可知,火工品的可靠性试验数据可统一表示成式(1)的形式。

3 基于广义线性模型的可靠性数据分析

3.1 感度分布函数

火工品的临界刺激量分布函数又称感度分布函数。工程中常用的火工品感度分布为正态分布、对数正态分布、Logistic分布和对数Logistic分布[3],其中对数正态分布和对数Logistic分布可通过对数变换,分别转换为正态分布和Logistic分布[10]。对于正态分布和Logistic分布,记火工品感度为X,可用位置-刻度模型[10]来统一表示：

X=μ+σε

(2)

式中，μ为位置参数,σ为刻度参数,ε为分布函数是G(·)的随机变量,其中G(·)与位置参数及刻度参数无关。

给定工作刺激量x,由式(2)可得火工品的可靠度函数：

(3)

3.2 分布参数极大似然估计

(4)

由式(4)的对数似然函数可以求得感度分布参数μ和σ的极大似然估计。由于没有解析解,一般利用解非线性方程组的数值方法来求解。然而数值方法中迭代算法的收敛性依赖于迭代初始值,如何确定适当的初始值以确保算法收敛一直是个难题。

为此,本研究利用广义线性模型[5-6]给出了一种有效的算法,用于求解感度分布参数的极大似然估计。

(5)

由广义线性模型的性质可知,式(5)对数似然函数可看成连接函数为G(β1+β2x)的二项分布变量的广义线性表达式[5]。对于广义线性模型,一般利用迭代加权方法来求得未知参数的极大似然估计。令β=(β1,β2)′,假设单个响应变量ni(i=1,2,…,k)的对数似然函数为li(ρi),其中ρi=β1+β2xi,则式(5)的似然函数可变换为:

(6)

(7)

(8)

(9)

(10)

(11)

式(11)可表示为矩阵形式:

I(β)=U′W(β)U

(12)

(13)

式中，依赖于初值β的部分已被抑制,通过迭代直到满足预定的精度要求,可得参数β的极大似然估计。

3.3 可靠性分析

(14)

在工程应用中,通常还需要获得可靠度置信下限。由于难以获得R(x)的分布,无法直接利用传统的区间估计方法。然而,相对而言一般比较容易获得近似区间估计,而且大多数情况它们没有显著差异,故通常利用近似区间估计来代替。

(15)

(16)

(17)

4 实验验证

例1某针刺雷管可靠性指标为：置信水平γ=0.95,可靠度R≥0.999,工作刺激量为：落锤重量(52±1) g,落高6 cm。

对该火工品进行三组升降法试验,试验结果如表1所示。

表1 某针刺雷管升降法试验数据
Table 1 Up-down test data of stab detonator

stimulus/cmgroupAresponsenumbernonresponsenumbergroupBresponsenumbernonresponsenumbergroupCresponsenumbernonresponsenumber1．00101011．517110172．071410117142．51431131433．0303030

由于步进法试验样本量较大,其参数估计较为稳定,将本方法的分析结果与利用步进法试验数据的分析结果进行对比,以验证本方法的合理性,该火工品的步进法试验数据如表2所示。

给定置信水平为γ=0.95和工作刺激量水平x=6 cm,对表2中的试验数据进行分析可得可靠性下限R2L=0.9998,故该火工品的可靠性达到了指标要求。

通过对两种方法给出的针刺雷管可靠性下限进行对比分析可知,本方法给出的可靠性下限为0.9996,略低于大样本方法给出的0.9998。经理论分析可知,当试验样本来源于同一总体时,试验样本量越大,参数置信区间越窄,故本研究给出的结果略低于大样本方法是合理的。

表2 某针刺雷管步进法试验数据
Table 2 Run-down test results of stab detonator

stimulus/cm11．522．533．544．5testnumber400200200200200200200200responsenumber02489138184191199200

例2某电雷管可靠度指标为：置信水平γ=0.95,可靠度R≥0.9999,工作刺激量为：700 mA电流。

对该火工品进行三组升降法试验,试验结果如表3所示。

表3 某电雷管升降法试验数据
Table 3 Up-down test results of electronic detonator

stimulus/cmgroupAresponsenumbernonresponsenumbergroupBresponsenumbernonresponsenumbergroupCresponsenumbernonresponsenumber300010101320171101834071210981136012494114380414141400101010

同样地将本方法的分析结果与利用步进法试验数据的分析结果进行对比,收集了该火工品试验样本量为1800发的步进法试验数据。给定置信水平为γ=0.95和工作刺激量水平x=700 mA,利用步进法试验数据可得可靠性下限R4L=0.99999,故该火工品的可靠性达到了指标要求。

通过对两种方法给出的电雷管可靠性下限进行对比分析可知,本文方法给出的可靠性下限为0.99998,略低于大样本方法给出的0.99999,其原因与针刺雷管可靠性分析结果一致。

结合上述两种火工品的分析结果,对比分析可知本文方法是合理可行的,可以利用约150个样本实现对可靠性要求为0.999以上的火工品进行有效地可靠性评价。

5 结 论

根据不同类型感度试验数据的共性特点和内在规律,利用广义线性模型来描述火工品可靠性与感度的关系,给出了适用于火工品感度试验数据的参数估计方法及相应算法,有效地改善了参数估计算法的收敛性过于依赖初始值的缺点。同时,结合参数的渐近正态性质,给出了火工品可靠度置信下限,提高了可靠性分析精度。分别以某针刺雷管和电雷管的为例,综合利用其多组升降法试验数据进行可靠性分析,并与利用大样本步进法试验数据分析结果进行对比,结果表明本方法合理可行,可以利用约150个样本实现对可靠性要求为0.999以上的火工品进行有效地可靠性评价。

参考文献:

[1] Dixon W J,Mood H M.A method for obtaining and analyzing sensitivity data[J].JournaloftheAmericanStatisticalAssociation,1948,43(241)：109-126.

[2] 洪东跑,赵宇,温玉全.基于序约束的火工品可靠性试验数据分析[J].含能材料,2008,16(5)：556-559.

HONG Dong-pao,ZHAO Yu,WEN Yu-quan.Order restricted analysis of reliability tests for explosive initiator[J].ChineseJournalofEnergeticMaterials(HannengCailiao),2008,16(5)：556-559.

[3] 蔡瑞娇,翟志强,董海平,等.火工品可靠性评估试验信息熵等值方法[J].含能材料,2007,15(1)：79-82.

CAI Rui-jiao,ZHAI Zhi-qiang,DONG Hai-ping,et al.Assessment method for reliability of initiating devices based on test information entropy equivalence[J].ChineseJournalofEnergeticMaterials(HannengCailiao),2007,15(1)：79-82.

[4] Barry T Neyer.A D-optimality-based sensitivity test[J].Technometrics,1994,1：61-70.

[5] Uusipaikka E.Confidence intervals in generalized regression models[M].London：Chapman&Hall/CRC,2009：176-178.

[6] McCullagh P,Nelder J A.Generalized linear models[M].Landon：Chapman and Hall,1989：130-135.

[7] Wendai W,Dimitri B K.Fitting the weibull log-linear model to accelerated life-test data[J].IEEETransactiononReliability,2000,49(2)：217-223.

[8] Jeff Wu.Efficient sequential designs with binary data[J].JournalofAmericaStatisticalAssociation,1985 (8)：974-984.

[9] Chao M T,Fuh C D.Bootstrap method for the up and down test on pyrotechnology sensitivity analysis[J].StatisticaSinica,2001 (11)：1-21.

[10] Balakrishnan N,Nevzorov V B.A Primer on statistical distributions [M].Hoboken：John Wiley & Sons,2003.

[11] 洪东跑,马小兵,赵宇.利用变环境数据的Weibull分布可靠性综合评估[J].北京航空航天大学学报,2012,38(11)：1485-1491.

HONG Dong-pao,MA Xiao-bing,ZHAO Yu.Integrated reliability assessment for weibull distribution using varied environment data[J].JournalofBeijingUniversityofAeronauticsandAstronautics,2012,38(11):1485-1491

[12] Bagdonavicius V,Nikulin M.Accelerated life models modeling and statistical analysis[M].London：Chapman&Hall/CRC,2002：120-125.

[13] Kalbfleisch J D,Prentice R L.The statistical analysis of failure time data[M].Hoboken：John Wiley & Sons,2002：105-108.