基于廣義線性模型的火工品可靠性數(shù)據(jù)分析

洪東跑,王英華,管 飛,溫玉全

(1.中國(guó)運(yùn)載火箭技術(shù)研究院,北京 100076; 2.北京理工大學(xué)爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100081)

1 引 言

在航天領(lǐng)域,火工品承擔(dān)著一系列有關(guān)飛行器點(diǎn)火、姿態(tài)控制、分離等控制指令的執(zhí)行,對(duì)飛行器飛行成敗起到關(guān)鍵性的作用,是飛行器的關(guān)鍵元件?；鸸て房煽啃砸蟾?需要有效的評(píng)價(jià)方法，作為高可靠性要求的敏感性產(chǎn)品,在國(guó)內(nèi)外的工程應(yīng)用中,一般利用感度試驗(yàn)數(shù)據(jù)來(lái)進(jìn)行可靠性分析。即假設(shè)每個(gè)產(chǎn)品均存在一個(gè)臨界刺激量,當(dāng)外界施加的刺激量大于臨界刺激量時(shí),該產(chǎn)品“響應(yīng)”; 否則“不響應(yīng)”[1]。同時(shí)假定臨界刺激量服從某分布,按照一定試驗(yàn)方法施加不同的刺激量進(jìn)行感度試驗(yàn),并利用試驗(yàn)數(shù)據(jù)進(jìn)行可靠性分析[2-3]。由于火工品的臨界刺激量不能直接測(cè)量,只能在若干刺激量下通過(guò)刺激試驗(yàn)來(lái)測(cè)定。每個(gè)火工品樣本只能做一次試驗(yàn),即使施加了刺激量不響應(yīng),也會(huì)破壞其結(jié)構(gòu)[4]。故感度試驗(yàn)所獲得的數(shù)據(jù)為施加的刺激量和相應(yīng)的“響應(yīng)”或“不響應(yīng)”數(shù),無(wú)法獲得精確的臨界刺激量。同時(shí)不同的方法獲得的試驗(yàn)數(shù)據(jù)差異較大,數(shù)據(jù)分析方法也不完全一致。特別當(dāng)樣本量較小時(shí),分析結(jié)果可能會(huì)有較大的偏差,影響了火工品可靠性分析的精度。而且利用感度試驗(yàn)數(shù)據(jù)進(jìn)行參數(shù)估計(jì)時(shí)無(wú)法獲得解析解,而利用數(shù)值算法求解時(shí),其收斂性依賴于迭代初始值。由廣義線性模型的相關(guān)研究[5-7]可知,該模型可有效地彌補(bǔ)數(shù)值算法中參數(shù)估計(jì)迭代不收斂或收斂慢的不足。

為此,本研究根據(jù)火工品可靠性數(shù)據(jù)分析的工程需求,探索了不同類型感度試驗(yàn)數(shù)據(jù)的共性特點(diǎn)和內(nèi)在規(guī)律,采用廣義線性模型描述了火工品可靠性與感度的關(guān)系,給出一種火工品可靠性數(shù)據(jù)分析方法,在小樣本下實(shí)現(xiàn)對(duì)高可靠性要求的火工品可靠性分析。

2 感度試驗(yàn)數(shù)據(jù)特點(diǎn)分析

工程中常用的感度試驗(yàn)方法有升降法[1]、Langlie法[2]、E方法[8],D優(yōu)化方法[4]等。由于升降法試驗(yàn)操作相對(duì)簡(jiǎn)單,是工程中最為常用的火工品感度試驗(yàn)方法[9]。為此,本研究以升降法試驗(yàn)為例,分析了火工品感度試驗(yàn)數(shù)據(jù)的特點(diǎn)。

升降法試驗(yàn)包括三個(gè)因素：試驗(yàn)樣本量N、初始刺激量x0和步長(zhǎng)d。x0和d選定后,用x0作第一次刺激-響應(yīng)試驗(yàn); 第二次及以后每次試驗(yàn)所用刺激量的取法如下：如前一次試探的反應(yīng)結(jié)果為“響應(yīng)”,則本次試探用刺激量為xi+1=xi-d; 如為“不響應(yīng)”,則為xi+1=xi+d。如此循環(huán)試驗(yàn),至完成預(yù)定試驗(yàn)樣本量N為止。對(duì)升降法試驗(yàn)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,按刺激量的升序排列,可表示成如下形式：

(1)

式中，k為刺激量個(gè)數(shù),xi(i=1,2,…,k)為試驗(yàn)刺激量,mi為在xi試驗(yàn)的不響應(yīng)數(shù),ni為在xi試驗(yàn)的響應(yīng)數(shù)。

本研究對(duì)Langlie法、E方法、D優(yōu)化方法等試驗(yàn)數(shù)據(jù)也進(jìn)行了分析,由分析結(jié)果可知,火工品的可靠性試驗(yàn)數(shù)據(jù)可統(tǒng)一表示成式(1)的形式。

3 基于廣義線性模型的可靠性數(shù)據(jù)分析

3.1 感度分布函數(shù)

火工品的臨界刺激量分布函數(shù)又稱感度分布函數(shù)。工程中常用的火工品感度分布為正態(tài)分布、對(duì)數(shù)正態(tài)分布、Logistic分布和對(duì)數(shù)Logistic分布[3],其中對(duì)數(shù)正態(tài)分布和對(duì)數(shù)Logistic分布可通過(guò)對(duì)數(shù)變換,分別轉(zhuǎn)換為正態(tài)分布和Logistic分布[10]。對(duì)于正態(tài)分布和Logistic分布,記火工品感度為X,可用位置-刻度模型[10]來(lái)統(tǒng)一表示：

X=μ+σε

(2)

式中，μ為位置參數(shù),σ為刻度參數(shù),ε為分布函數(shù)是G(·)的隨機(jī)變量,其中G(·)與位置參數(shù)及刻度參數(shù)無(wú)關(guān)。

給定工作刺激量x,由式(2)可得火工品的可靠度函數(shù)：

(3)

3.2 分布參數(shù)極大似然估計(jì)

(4)

由式(4)的對(duì)數(shù)似然函數(shù)可以求得感度分布參數(shù)μ和σ的極大似然估計(jì)。由于沒(méi)有解析解,一般利用解非線性方程組的數(shù)值方法來(lái)求解。然而數(shù)值方法中迭代算法的收斂性依賴于迭代初始值,如何確定適當(dāng)?shù)某跏贾狄源_保算法收斂一直是個(gè)難題。

為此,本研究利用廣義線性模型[5-6]給出了一種有效的算法,用于求解感度分布參數(shù)的極大似然估計(jì)。

(5)

由廣義線性模型的性質(zhì)可知,式(5)對(duì)數(shù)似然函數(shù)可看成連接函數(shù)為G(β1+β2x)的二項(xiàng)分布變量的廣義線性表達(dá)式[5]。對(duì)于廣義線性模型,一般利用迭代加權(quán)方法來(lái)求得未知參數(shù)的極大似然估計(jì)。令β=(β1,β2)′,假設(shè)單個(gè)響應(yīng)變量ni(i=1,2,…,k)的對(duì)數(shù)似然函數(shù)為li(ρi),其中ρi=β1+β2xi,則式(5)的似然函數(shù)可變換為:

(6)

(7)

(8)

(9)

(10)

(11)

式(11)可表示為矩陣形式:

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

(12)

(13)

式中，依賴于初值β的部分已被抑制,通過(guò)迭代直到滿足預(yù)定的精度要求,可得參數(shù)β的極大似然估計(jì)。

3.3 可靠性分析

(14)

在工程應(yīng)用中,通常還需要獲得可靠度置信下限。由于難以獲得R(x)的分布,無(wú)法直接利用傳統(tǒng)的區(qū)間估計(jì)方法。然而,相對(duì)而言一般比較容易獲得近似區(qū)間估計(jì),而且大多數(shù)情況它們沒(méi)有顯著差異,故通常利用近似區(qū)間估計(jì)來(lái)代替。

(15)

(16)

(17)

4 實(shí)驗(yàn)驗(yàn)證

例1某針刺雷管可靠性指標(biāo)為：置信水平γ=0.95,可靠度R≥0.999,工作刺激量為：落錘重量(52±1) g,落高6 cm。

對(duì)該火工品進(jìn)行三組升降法試驗(yàn),試驗(yàn)結(jié)果如表1所示。

表1 某針刺雷管升降法試驗(yàn)數(shù)據(jù)
Table 1 Up-down test data of stab detonator

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

由于步進(jìn)法試驗(yàn)樣本量較大,其參數(shù)估計(jì)較為穩(wěn)定,將本方法的分析結(jié)果與利用步進(jìn)法試驗(yàn)數(shù)據(jù)的分析結(jié)果進(jìn)行對(duì)比,以驗(yàn)證本方法的合理性,該火工品的步進(jìn)法試驗(yàn)數(shù)據(jù)如表2所示。

給定置信水平為γ=0.95和工作刺激量水平x=6 cm,對(duì)表2中的試驗(yàn)數(shù)據(jù)進(jìn)行分析可得可靠性下限R2L=0.9998,故該火工品的可靠性達(dá)到了指標(biāo)要求。

通過(guò)對(duì)兩種方法給出的針刺雷管可靠性下限進(jìn)行對(duì)比分析可知,本方法給出的可靠性下限為0.9996,略低于大樣本方法給出的0.9998。經(jīng)理論分析可知,當(dāng)試驗(yàn)樣本來(lái)源于同一總體時(shí),試驗(yàn)樣本量越大,參數(shù)置信區(qū)間越窄,故本研究給出的結(jié)果略低于大樣本方法是合理的。

表2 某針刺雷管步進(jìn)法試驗(yàn)數(shù)據(jù)
Table 2 Run-down test results of stab detonator

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

例2某電雷管可靠度指標(biāo)為：置信水平γ=0.95,可靠度R≥0.9999,工作刺激量為：700 mA電流。

對(duì)該火工品進(jìn)行三組升降法試驗(yàn),試驗(yàn)結(jié)果如表3所示。

表3 某電雷管升降法試驗(yàn)數(shù)據(jù)
Table 3 Up-down test results of electronic detonator

stimulus/cmgroupAresponsenumbernonresponsenumbergroupBresponsenumbernonresponsenumbergroupCresponsenumbernonresponsenumber300010101320171101834071210981136012494114380414141400101010

同樣地將本方法的分析結(jié)果與利用步進(jìn)法試驗(yàn)數(shù)據(jù)的分析結(jié)果進(jìn)行對(duì)比,收集了該火工品試驗(yàn)樣本量為1800發(fā)的步進(jìn)法試驗(yàn)數(shù)據(jù)。給定置信水平為γ=0.95和工作刺激量水平x=700 mA,利用步進(jìn)法試驗(yàn)數(shù)據(jù)可得可靠性下限R4L=0.99999,故該火工品的可靠性達(dá)到了指標(biāo)要求。

通過(guò)對(duì)兩種方法給出的電雷管可靠性下限進(jìn)行對(duì)比分析可知,本文方法給出的可靠性下限為0.99998,略低于大樣本方法給出的0.99999,其原因與針刺雷管可靠性分析結(jié)果一致。

結(jié)合上述兩種火工品的分析結(jié)果,對(duì)比分析可知本文方法是合理可行的,可以利用約150個(gè)樣本實(shí)現(xiàn)對(duì)可靠性要求為0.999以上的火工品進(jìn)行有效地可靠性評(píng)價(jià)。

5 結(jié) 論

根據(jù)不同類型感度試驗(yàn)數(shù)據(jù)的共性特點(diǎn)和內(nèi)在規(guī)律,利用廣義線性模型來(lái)描述火工品可靠性與感度的關(guān)系,給出了適用于火工品感度試驗(yàn)數(shù)據(jù)的參數(shù)估計(jì)方法及相應(yīng)算法,有效地改善了參數(shù)估計(jì)算法的收斂性過(guò)于依賴初始值的缺點(diǎn)。同時(shí),結(jié)合參數(shù)的漸近正態(tài)性質(zhì),給出了火工品可靠度置信下限,提高了可靠性分析精度。分別以某針刺雷管和電雷管的為例,綜合利用其多組升降法試驗(yàn)數(shù)據(jù)進(jìn)行可靠性分析,并與利用大樣本步進(jìn)法試驗(yàn)數(shù)據(jù)分析結(jié)果進(jìn)行對(duì)比,結(jié)果表明本方法合理可行,可以利用約150個(gè)樣本實(shí)現(xiàn)對(duì)可靠性要求為0.999以上的火工品進(jìn)行有效地可靠性評(píng)價(jià)。

參考文獻(xiàn):

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

[2] 洪東跑,趙宇,溫玉全.基于序約束的火工品可靠性試驗(yàn)數(shù)據(jù)分析[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] 蔡瑞嬌,翟志強(qiáng),董海平,等.火工品可靠性評(píng)估試驗(yàn)信息熵等值方法[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] 洪東跑,馬小兵,趙宇.利用變環(huán)境數(shù)據(jù)的Weibull分布可靠性綜合評(píng)估[J].北京航空航天大學(xué)學(xué)報(bào),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.