火工品發(fā)火可靠性小樣本鑒定試驗(yàn)方法

田玉斌，王典鵬

(北京理工大學(xué) 理學(xué)院，北京100081)

火工品是武器裝備的關(guān)鍵元件，其發(fā)火可靠性指標(biāo)往往很高，如可靠度下限p =0.999，置信水平0.95.為了鑒定火工品的發(fā)火可靠性，目前使用的方法有計(jì)數(shù)法［1］、升降法［2］和計(jì)數(shù)--計(jì)量綜合方法［3］。

計(jì)數(shù)法直接在工作水平處進(jìn)行試驗(yàn)，基于二項(xiàng)分布模型給出發(fā)火可靠度的置信下限。計(jì)數(shù)法簡單，所需樣本量很大，如為了鑒定火工品的發(fā)火可靠度不低于0.999(置信水平0.95)，需試驗(yàn)2996 發(fā)火工品，且無一失效［1］。

升降法做如下假設(shè): 1)每一產(chǎn)品存在一個(gè)固有的臨界刺激量，當(dāng)外界施加的刺激水平大于該臨界刺激量時(shí)，產(chǎn)品發(fā)火; 否則不發(fā)火。2)產(chǎn)品臨界刺激量X 的分布F(x)稱為感度分布，常取

式中:G(·)為已知函數(shù); μ 為位置參數(shù);σ 為刻度參數(shù)。在升降法中，按照升降試驗(yàn)規(guī)則，獲得發(fā)火或不發(fā)火數(shù)據(jù);通過計(jì)算參數(shù)(μ，σ)的最大似然估計(jì)，求出火工品發(fā)火可靠度的點(diǎn)估計(jì)和近似置信下限。由于升降法使用了感度分布模型，樣本量大大減小。但是，升降法的估計(jì)精度不高，往往高估火工品的發(fā)火可靠性。

在計(jì)數(shù)—計(jì)量綜合方法中，首先利用升降法估計(jì)參數(shù)(μ，σ)，擬合感度分布。然后，基于該擬合分布，按照一定原理，尋找與工作水平等價(jià)的試驗(yàn)水平，該水平對應(yīng)的發(fā)火率較低。最后，在該水平處進(jìn)行試驗(yàn)，通過考察該水平處的發(fā)火率，等價(jià)地鑒定工作水平處的發(fā)火可靠性。為了鑒定火工品發(fā)火可靠度不低于0.999(置信水平0.95)［3］，計(jì)數(shù)-計(jì)量綜合方法所需樣本量在300 發(fā)左右。

遵循目前國際相關(guān)科學(xué)研究的前沿方向［4-5］，研究的熱點(diǎn)不僅僅局限于鑒定火工品的發(fā)火可靠性，更多關(guān)注的是: 在中小樣本情形，推斷火工品以預(yù)定概率p 發(fā)火的刺激水平xp(也稱為火工品100% p 發(fā)火點(diǎn))，它滿足

推斷包括:給出xp的點(diǎn)估計(jì)和置信上限。這樣的方法，不僅能夠鑒定火工品的發(fā)火可靠性，而且能夠定量給出火工品的作用裕度。

本文將基于估計(jì)xp的優(yōu)化試驗(yàn)設(shè)計(jì)，給出xp的點(diǎn)估計(jì)和置信上限、作用裕度的區(qū)間估計(jì)、以及鑒定火工品發(fā)火可靠性的方法。最后，應(yīng)用蒙特卡羅方法模擬研究上述方法的精確性。

1 鑒定火工品發(fā)火可靠性的優(yōu)化試驗(yàn)設(shè)計(jì)

1.1 數(shù)學(xué)模型

假設(shè)火工品的感度分布如(1)式所示，對一產(chǎn)品施加刺激水平x，記試驗(yàn)結(jié)果:yx=1 為發(fā)火;yx=0 為不發(fā)火。在水平x 處進(jìn)行試驗(yàn)，產(chǎn)品發(fā)火的概率

此時(shí)，xp=μ+σG-1(p).為了較精確地推斷xp，應(yīng)采用序貫試驗(yàn)設(shè)計(jì)［5］。假設(shè)(x1，y1)…(xn，yn)是對應(yīng)的試驗(yàn)數(shù)據(jù)，其中xi為第i 次試驗(yàn)的水平，yi為試驗(yàn)結(jié)果，n 為樣本容量。基于該觀測數(shù)據(jù)，似然函數(shù)［6］為

式中:pi=G((xi-μ)/σ.對(4)式關(guān)于μ，σ 求導(dǎo)數(shù)，并求期望

構(gòu)成Fisher 信息矩陣。該矩陣是試驗(yàn)數(shù)據(jù)包含(μ，σ)的信息的一種度量［7］。

1.2 推斷xp的優(yōu)化試驗(yàn)設(shè)計(jì)

根據(jù)文獻(xiàn)［7］，在合理的試驗(yàn)水平范圍內(nèi)，優(yōu)化分散試驗(yàn)水平，可以精確地估計(jì)參數(shù)(μ，σ).根據(jù)文獻(xiàn)［4］，將試驗(yàn)水平逐步設(shè)置在xp附近，將提高xp估計(jì)的精度。文獻(xiàn)［8］綜合這2 種思想，給出3 段式優(yōu)化試驗(yàn)設(shè)計(jì)。本文應(yīng)用該優(yōu)化試驗(yàn)設(shè)計(jì)，給出xp的點(diǎn)估計(jì)和置信上限、火工品作用裕度的區(qū)間估計(jì)、以及鑒定火工品發(fā)火可靠性的方法。

1.2.1 第1 段設(shè)計(jì)

第1 段設(shè)計(jì)的目的是: 合理確定試驗(yàn)水平的范圍。首先，猜測μ 的取值范圍(μmin，μmax)和σ 的值σg.以下將對(μmin，μmax)快速進(jìn)行調(diào)整，使試驗(yàn)水平較對稱地設(shè)置在μ 周圍。

在x1=3/4μmin+ 1/4μmax和x2=1/4μmin+ 3/4μmax兩水平進(jìn)行試驗(yàn)，獲得試驗(yàn)結(jié)果(y1，y2):

1)若(y1，y2)=(0，0)，在x3=μmax+1.5σg處進(jìn)行試驗(yàn)，若y3=1，跳入第2 段設(shè)計(jì); 若y3=0，在x4=μmax+3σg處進(jìn)行試驗(yàn)。如此繼續(xù)，直到某次試驗(yàn)結(jié)果為1 時(shí)，跳入第2 段設(shè)計(jì);否則在上一次試驗(yàn)水平的基礎(chǔ)上增加1.5 倍σg，繼續(xù)試驗(yàn)，直到出現(xiàn)發(fā)火結(jié)果。

2)若(y1，y2)=(1，1)，在x3=μmin-1.5σg處進(jìn)行試驗(yàn)，若y3=0，跳入第2 段設(shè)計(jì);若y3=1，則在x4=μmin-3σg處進(jìn)行試驗(yàn)。如此繼續(xù)，直到某次試驗(yàn)結(jié)果為0 時(shí)，跳入第2 段設(shè)計(jì);否則在上一次試驗(yàn)水平的基礎(chǔ)上減少1.5 倍σg，繼續(xù)試驗(yàn)，直到出現(xiàn)不發(fā)火結(jié)果。

3)若(y1，y2)=(0，1)，跳入第2 段設(shè)計(jì)。

4)若(y1，y2)=(1，0)，分別在x3=μmin-3σg和x4=μmax+3σg處各進(jìn)行1 試驗(yàn)。然后，跳入第2段設(shè)計(jì)。

1.2.2 第2 段設(shè)計(jì)

第2 段設(shè)計(jì)的目的是:合理布置試驗(yàn)水平，優(yōu)化模型信息。將該段設(shè)計(jì)分成3 部分，分別為尋找交錯(cuò)區(qū)間、加強(qiáng)交錯(cuò)區(qū)間和強(qiáng)化信息部分。

1)尋找交錯(cuò)區(qū)間。尋找交錯(cuò)區(qū)間分2 個(gè)步驟:

①計(jì)算出現(xiàn)發(fā)火結(jié)果的最小水平m1，以及出現(xiàn)不發(fā)火結(jié)果的最大水平M0.(a)若m1＜M0，停止試驗(yàn)，轉(zhuǎn)至步驟2); (b)若m1＞M0，在σ =σg的條件下，求μ 的最大似然估計(jì)，并在該點(diǎn)進(jìn)行試驗(yàn)，直至m1＜M0或m1-M0＜1.5σg.若m1＜M0，轉(zhuǎn)至步驟2);(c)在m1-M0＜1.5σg，而且m1＞M0情形，分別在m1+0.3σg和M0-0.3σg處進(jìn)行試驗(yàn)，如果出現(xiàn)不發(fā)火結(jié)果或發(fā)火結(jié)果，轉(zhuǎn)至步驟2);否則，轉(zhuǎn)至步驟②;

②令σg=2/3σg.更新①的m1和M0，并重復(fù)步驟①直至m1＜M0成立。

2)加強(qiáng)交錯(cuò)區(qū)間。如果M0- m1≥σg，則在(M0+m1)/2 處進(jìn)行試驗(yàn);如果0＜M0-m1＜σg，在(M0+m1)/2 +0.5σg和(M0+m1)/2 -0.5σg處各進(jìn)行1 次試驗(yàn)。

3)優(yōu)化模型信息。應(yīng)用D 最優(yōu)方法［7-9］分散試驗(yàn)水平，優(yōu)化模型信息。①利用已觀測到的試驗(yàn)數(shù)據(jù)(x1，y1)…(xk，yk)，計(jì)算(μ，σ)的最大似然估計(jì)(對該估計(jì)做如下限制性調(diào)整

式中: ak=min{ x1，…，xk}; bk=max{ x1，…，xk};②選擇新的試驗(yàn)水平xk+1，使得Fisher 信息矩陣的行列式在(k，actual，k，actual)處達(dá)到最大; ③步驟重復(fù)①～②，直到完成預(yù)定樣本量n1的試驗(yàn)。

1.2.3 第3 段設(shè)計(jì)

第3 段設(shè)計(jì)的目的是安排試驗(yàn)水平，使其更多、更快地聚集在xp附近。該段設(shè)計(jì)分為2 部分。

1)確定初始值。利用已觀測到的數(shù)據(jù)(x1，y1)，…，(xn1，yn1)，計(jì)算最小試驗(yàn)水平xmin和最大試驗(yàn)水平xmax，求出μ 和σ 最大似然估計(jì)(，)，對其進(jìn)行修正，即

2)優(yōu)化逼近xp.在第3 段設(shè)計(jì)中，首先在xn1+1處進(jìn)行試驗(yàn)?？紤]如下形式的試驗(yàn)設(shè)計(jì)

式中:xn1+i為第3 段設(shè)計(jì)的第i 次試驗(yàn)的水平; yn1+i為相應(yīng)的試驗(yàn)結(jié)果; 且

式中:Φ(·)為標(biāo)準(zhǔn)正態(tài)分布函數(shù); φ(·)為標(biāo)準(zhǔn)正態(tài)密度函數(shù)。

按(6)式和(7)式安排試驗(yàn)，直至完成容量為n2的試驗(yàn)。此時(shí)，xp的點(diǎn)估計(jì)為=xn+1，n=n1+n2.

1.3 火工品發(fā)火可靠性鑒定

假設(shè)火工品的工作水平為A，那么火工品的發(fā)火可靠度R 為該水平處的發(fā)火概率R=G((A-μ)/σ).按照文獻(xiàn)［11 -12］，針對火工品的發(fā)火特性，作用裕度γ =(A -θ)/θ.定理1 在置信水平1 -α下，給出xp的置信上限AU，定理2 給出鑒定火工品的發(fā)火可靠度不低于p 的原理，定理3 給出γ 的區(qū)間估計(jì)。

定理1 利用(6)式和(7)式，獲得試驗(yàn)數(shù)據(jù)(x1，y1)，…，(xn，yn).在置信水平1 -α 下，xp的置信上限為

式中:u1-α為標(biāo)準(zhǔn)正態(tài)分布的下側(cè)1 -α 分位數(shù)。

證明1 令θ =xp.根據(jù)文獻(xiàn)［4］，給定xn+1，θ的分布為由該分布，為了使

成立，AU需滿足

定理2 利用(6)式和(7)式，獲得試驗(yàn)數(shù)據(jù)(x1，y1)，…，(xn，yn).如果A≥AU，那么在置信水平1 -α 下，產(chǎn)品的可靠度R 不低于p.

證明2 根據(jù)定理1

1 -α=P(AU≥θ)，

于是，有

所以，

1 -α≤P(R≥p).

定理3 利用(6)式和(7)式，獲得試驗(yàn)數(shù)據(jù)(x1，y1)，…，(xn，yn).則在置信水平1 -α 下，γ 的區(qū)間估計(jì)為

取

式中: uα/2為標(biāo)準(zhǔn)正態(tài)分布的下側(cè)α/2 分位數(shù);u1-α/2為標(biāo)準(zhǔn)正態(tài)分布的下側(cè)1 -α/2 分位數(shù)。于是

2 模擬與應(yīng)用研究

該節(jié)主要模擬第1 節(jié)提出的優(yōu)化試驗(yàn)設(shè)計(jì)，給出火工品99.9%發(fā)火點(diǎn)x0.999的點(diǎn)估計(jì)、置信上限和該上限對x0.999真值的覆蓋率。

2.1 模擬

在模擬研究中，假設(shè)感度分布為F(x)=Φ((x-μ)/σ)，μ =10，σ =1.此時(shí)，99.9%發(fā)火點(diǎn)的真值為x0.999=13.0905.按照優(yōu)化試驗(yàn)設(shè)計(jì)規(guī)則，在試驗(yàn)前猜測μ 的取值范圍(μmin，μmax)和σ 的值σg.在實(shí)際應(yīng)用中，根據(jù)專家經(jīng)驗(yàn)，對μ 的猜測一般較準(zhǔn)確，對σ 的猜測波動(dòng)較大。所以，在模擬中μg取9，10，11; σg取0.5，1，1.5，2 和3; (μmin，μmax)=(μg-3σg，μg+3σg);樣本量n 取80 和100，其中第1 段和第2 段設(shè)計(jì)的總樣本量n1=30，第3 段設(shè)計(jì)的樣本量n2=50 或70.

在各種初始猜測下，模擬優(yōu)化試驗(yàn)設(shè)計(jì)，獲得容量為n 的數(shù)據(jù): 計(jì)算x0.999的點(diǎn)估計(jì)xn+1; 在置信水平1 -α=0.95 下，計(jì)算xp的置信上限AU;將上述過程重復(fù)1 000 次，計(jì)算xp估計(jì)的均值，AU的均值，計(jì)算{θ＜AU}的頻率，并與名義覆蓋率1 -α =0.95 比較。模擬結(jié)果如表1～表3.

表1 μguess =9，σguess =0.5、1、1.5、2 和3 猜測下的估計(jì)結(jié)果Tab.1 Estimation results in μguess =9，σguess =0.5，1，1.5，2 and 3

表2 μguess =10，σguess =0.5、1、1.5、2 和3 猜測下的估計(jì)結(jié)果Tab.2 Estimation results in μguess =10，σguess =0.5，1，1.5，2 and 3

表3 μguess =11，σguess =0.5、1、1.5、2 和3 猜測下的估計(jì)結(jié)果Tab.3 Estimation results in μguess =11，σguess =0.5，1，1.5，2 and 3

2.2 模擬研究結(jié)論

從上述表格可知，針對火工品的高發(fā)火可靠性指標(biāo)(置信水平0.95，可靠度下限0.999)，在樣本容量n=100 時(shí)，在各種可能的初始猜測下，應(yīng)用優(yōu)化試驗(yàn)設(shè)計(jì)收集數(shù)據(jù)，給出0.999 發(fā)火點(diǎn)的點(diǎn)估計(jì)和置信上限，點(diǎn)估計(jì)與真值偏差很小，置信上限均高于真值且接近真值，覆蓋率接近0.95 且高于0.95.這些都表明，應(yīng)用本文提出的試驗(yàn)方法，在n=100 時(shí)，x0.999的點(diǎn)估計(jì)、置信上限、以及置信上限的覆蓋率對初始猜測都是穩(wěn)健的，點(diǎn)估計(jì)較精確，置信上限可信，其覆蓋率大于0.95，表明置信上限略顯保守。

綜合而言，應(yīng)用本文提出的試驗(yàn)方法鑒定火工品的發(fā)火可靠性，樣本容量比現(xiàn)有方法大幅減小，結(jié)果可信且較精確。

References)

［1］GJB 376—87.火工品可靠性評估方法［S］.北京:國防科學(xué)技術(shù)工業(yè)委員會(huì)，1987.GJB376—87.Assessment method of reliability of initiating devices［S］.Beijing: National Standard Press，1987.(in Chinese)

［2］GJB/Z377A—94.感度試驗(yàn)用數(shù)據(jù)統(tǒng)計(jì)方法［S］.北京: 國防科學(xué)技術(shù)工業(yè)委員會(huì)，1995.GJB/Z377A—94.Statistical methods for sensitivity tests［S］.Beijing: National Standard Press，1995.(in Chinese)

［3］GJB6478—2008.火工品計(jì)數(shù)—計(jì)量可靠性評估方法［S］.北京:國防科學(xué)技術(shù)工業(yè)委員會(huì)，2008.GJB6478—2008.Variables-attributes synthetic assessment method of reliability of initiating devices［S］.Beijing: National Standard Press，2008.(in Chinese)

［4］Joseph V R.Efficient robbins-monro procedure for binary data［J］.Biometrika，2004，91: 461 -470.

［5］Joseph V R，Tian Yubin，Wu C F J.Adaptive design for stochastic rooting-finding［J］.Statistics Sinica，2007，17:1549 -1565.

［6］Silvapulle M J.On the existence of maximum likelihood estimators for the binomial response model［J］.Journal of the Royal Statistical Society，Ser B，1981，43: 310 -313.

［7］Neyer B T.A d-optimality-based sensitivity test［J］.Technometrics，1994，36: 61 -70.

［8］Jeff Wu，Tian Yubin.The optimal design for root finding with binary response data［C］.The First International Conference on the Interface between Statistics and Engineering，Beijing:Beijing University of Technology Press，2009: 9 -10.

［9］袁俊明，劉玉存.Neyer D-最優(yōu)化的新感度試驗(yàn)方法研究［J］.火工品，2005，2: 24 -27.YUAN Jun-ming，LIU Yu-cun.Study on new neyer doptimal sensitivity test［J］.Initiators ＆ Pyrotechnics，2005，2:24 -27.(in Chinese)

［10］Tian Yubin，Li Guoying，F(xiàn)ang Yongfei.The synthetically analytical method for data sets on pyrotechnics reliability test［J］.J Sys Sci ＆ Math Scis，2006，26: 147 -158.

［11］Laurence J.Bement.Pyrotechnic system failure: cause and prevention［R］.NASA TM100633，Washington D C: 1988.

［12］王鵬，杜志明.火工煙火裝置裕度研究與設(shè)計(jì)方法綜述［J］.火工品，2005，2: 34 -38.WANG Peng，DU Zhi-ming.Summarize of margin research and design method of pyrotechnic devices［J］.Initiators ＆ Pyrotechnics，2005，2:34 -38.(in Chinese)