1 引言

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

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

2 感度试验数据特点分析

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

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

$ \left\{ \begin{array}{l} {x_1}, {x_2}, \cdots, {x_k}\\ {n_1}, {n_2}, \cdots, {n_k}\\ {m_1}, {m_2}, \cdots, {m_k} \end{array} \right\} $

(1)

式中, k为刺激量个数, x_i(i=1, 2, …, k)为试验刺激量, m_i为在x_i试验的不响应数, n_i为在x_i试验的响应数。

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

3 基于广义线性模型的可靠性数据分析 3.1 感度分布函数

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

$ X = \mu + \sigma \varepsilon $

(2)

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

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

$ R\left( x \right) = G\left( {\frac{{x - \mu }}{\sigma }} \right) $

(3)

如果x服从正态分布, 则G(·)为标准正态分布函数, G(t)=Φ(t); 如果x服从Logistic分布, 则G(·)为标准Logistic分布函数, $G\left( t \right) = \frac{{{\rm{exp}}\left( t \right)}}{{1 + {\rm{exp}}\left( t \right)}}$。

3.2 分布参数极大似然估计

已知感度分布函数$G\left( {\frac{{x - \mu }}{\sigma }} \right)$, 结合式(1)的感度试验数据可得似然函数:

$ L = \prod\limits_{i = 1}^k {{C_i}} {\left( {G\left( {\frac{{{x_i} - \mu }}{\sigma }} \right)} \right)^{{n_i}}}{\left( {1 - G\left( {\frac{{{x_i} - \mu }}{\sigma }} \right)} \right)^{{m_i}}} $

(4)

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

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

令${\beta _1} = - \frac{\mu }{\sigma }, {\beta _2} = \frac{1}{\sigma }$, 则感度分布函数可变换为G(β₁+β₂x), 代入式(4)可得对数似然函数:

$ l = \sum\limits_{i = 1}^k {({n_i}{\rm{ln}}\left( {G\left( {{\beta _1} + {\beta _2}{x_i}} \right)} \right) + {m_i}{\rm{ln}}(1 - G({\beta _1} + {\beta _2}{x_i})))} $

(5)

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

$ l\left( \mathit{\boldsymbol{\beta }} \right) = \sum\limits_{i = 1}^k {{l_i}({\rho _i})} $

(6)

记${f_i}\left( \mathit{\boldsymbol{\beta }} \right) = \frac{{\partial l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {{\rm{ \mathsf{ β} }}_i}}}(i = 1, 2)$, 则参数β的极大似然估计${\mathit{\boldsymbol{\hat \beta }}} $满足:

$ {f_i}\left( {\mathit{\boldsymbol{\hat \beta }}} \right) = 0\;\;\;\;\;\;\;\;(i = 1, 2) $

(7)

记$f\left( \mathit{\boldsymbol{\beta }} \right) = ({f_1}\left( \mathit{\boldsymbol{\beta }} \right), {f_2}\left( \mathit{\boldsymbol{\beta }} \right))\prime, \mathit{\boldsymbol{\rho }} = \left( {{\rho _1}, {\rho _2}} \right)\prime, \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{\beta }} \right) = \left( {\frac{{\partial l\left( {\mathit{\boldsymbol{\beta }}} \right)}}{{\partial {\rho _1}}}, \frac{{\partial l\left( {\mathit{\boldsymbol{\beta }}} \right)}}{{\partial {\rho _2}}}} \right)\prime $, 则有$f\left( \mathit{\boldsymbol{\beta }} \right) = \frac{{\partial \rho \prime }}{{\partial \mathit{\boldsymbol{\beta }}}}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{\beta }} \right) $, 故式(7)可表示为矩阵形式:

$ \frac{{\partial \mathit{\boldsymbol{\rho }}\prime }}{{\partial \mathit{\boldsymbol{\beta }}}}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{\beta }} \right)\left| {_{\mathit{\boldsymbol{\beta }} = \mathit{\boldsymbol{\hat \beta }}}} \right. = 0 $

(8)

在β处对$\frac{{\partial l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial \mathit{\boldsymbol{\beta }}}}$进行Taylor展开可得:

$ \frac{{\partial \mathit{\boldsymbol{\rho }}\prime }}{{\partial \mathit{\boldsymbol{\beta }}}}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{\beta }} \right) + \left[{\sum\limits_{i = 1}^2 {\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}}}} \frac{{{\partial ^2}l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\rho _i}^2}}\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}\prime }}} \right]\left( {\mathit{\boldsymbol{\hat \beta }} - \mathit{\boldsymbol{\beta }}} \right) \cong 0 $

(9)

记$\mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{\beta }} \right) = - \sum\limits_{i = 1}^2 {\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}}}\frac{{{\partial ^2}l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\rho _i}^2}}\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}\prime }}} $, 假设J(β)可逆, 则有:

$ \mathit{\boldsymbol{\hat \beta }} \cong \mathit{\boldsymbol{\beta }} + \mathit{\boldsymbol{J}}{\left( \mathit{\boldsymbol{\beta }} \right)^{ - 1}}\frac{{\partial \mathit{\boldsymbol{\rho }}\prime }}{{\partial \mathit{\boldsymbol{\beta }}}}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{\beta }} \right) $

(10)

因此给定初值β, 利用式(10)进行迭代可得极大似然估计${\mathit{\boldsymbol{\hat \beta }}}$, 该迭代方法是Newton-Raphson算法的特殊应用^[11]。在实际应用中, 用J(β)的期望I(β)来代替J(β)将更为方便, I(β)为:

$ \mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{\beta }} \right) = \sum\limits_{i = 1}^2 {\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}}}E\left( { - \frac{{{\partial ^2}l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\rho _i}^2}}} \right)\frac{{\partial {\rho _i}}}{{\partial \mathit{\boldsymbol{\beta }}\prime }}} $

(11)

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

$ \mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{\beta }} \right) = \mathit{\boldsymbol{U}}\prime \mathit{\boldsymbol{W}}\left( \mathit{\boldsymbol{\beta }} \right)\mathit{\boldsymbol{U}} $

(12)

式中, $\mathit{\boldsymbol{U}} = \frac{{\partial \rho }}{{\partial \mathit{\boldsymbol{\beta }}\prime }}, \mathit{\boldsymbol{W}}\left( \mathit{\boldsymbol{\beta }} \right)$为对角矩阵, 其第i个元素为$E\left( { - \frac{{{\partial ^2}l\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\rho _i}^2}}} \right)$。用I(β)代替J(β), 则式(10)可转换为:

$ \mathit{\boldsymbol{\hat \beta }} = {\left( {\mathit{\boldsymbol{U'WU}}} \right)^{ - 1}}\mathit{\boldsymbol{U'W}}\left( {\mathit{\boldsymbol{U\beta + }}{\mathit{\boldsymbol{W}}^{{\rm{ - 1}}}}\mathit{\boldsymbol{u}}} \right) $

(13)

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

当火工品感度服从正态分布, 可取连接函数为G(β₁+β₂x)=Φ(β₁+β₂x), 由式(13)可得参数极大似然估计${\mathit{\boldsymbol{\hat \beta }}}$; 当感度服从Logistic分布时, 取连接函数为 $G({\beta _1} + {\beta _2}x) = \frac{{{\rm{exp}}({\beta _1} + {\beta _2}x)}}{{1 + {\rm{exp}}({\beta _1} + {\beta _2}x)}}$, 由式(13)可得参数极大似然估计${\mathit{\boldsymbol{\hat \beta }}}$。已知参数估计${\mathit{\boldsymbol{\hat \beta }}}$, 可得感度分布参数估计$(\hat \sigma = {{\hat \beta }_2}, \hat \mu = - \frac{{{{\hat \beta }_1}}}{{{{\hat \beta }_2}}})$。

3.3 可靠性分析

给定工作刺激量x, 把感度分布参数估计$\left( {\hat u, \hat \sigma } \right)$代入式(3)可得火工品的可靠度估计:

$ \hat R\left( x \right) = G\left( {\frac{{x - \hat \mu }}{{\hat \sigma }}} \right) $

(14)

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

由Wald统计量^[12]可知, 在一定条件下, 极大似然估计${\mathit{\boldsymbol{\hat \beta }}} $具有渐近正态性

$ \frac{{\mathit{\boldsymbol{\hat \beta }} - \mathit{\boldsymbol{\beta }}}}{{\mathit{\boldsymbol{I}}{{\left( {\mathit{\boldsymbol{\hat \beta }}} \right)}^{ - 1}}}} \approx N(0, 1) $

(15)

其中$\mathit{\boldsymbol{I}}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)$为Fisher信息矩阵的估计。记h(β)为β的函数, 由式(15)可得:

$ \frac{{h\left( {\mathit{\boldsymbol{\hat \beta }}} \right) - h\left( \mathit{\boldsymbol{\beta }} \right)}}{{\sqrt {\mathit{\boldsymbol{J}}{\prime _h}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)\mathit{\boldsymbol{I}}{^{ - 1}}{\mathit{\boldsymbol{J}}_h}(\mathit{\boldsymbol{\hat \beta }})} }} \approx N(0, 1) $

(16)

式中, ${\mathit{\boldsymbol{J}}_{\rm{h}}}\left( \mathit{\boldsymbol{\beta }} \right) = \left( {\frac{{\partial h\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\mathit{\boldsymbol{\beta }}_1}}}, \frac{{\partial h\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\mathit{\boldsymbol{\beta }}_2}}}} \right)\prime $, 利用式(16)可得h(β)的近似区间估计^[13]。由于火工品可靠度R(x)为关于参数β的连续可导函数, 故可取h(β)=R(x), 由式(16)来求可靠度R(x)的置信区间。由于可靠度函数的取值区间为(0, 1], 而正态分布随机变量的取值为(-∞, ∞), 为了使参数估计更好地满足正态分布的性质, 需要对其进行logit变换, 取$h\left( \mathit{\boldsymbol{\beta }} \right) = \ln \left( {\frac{{R\left( x \right)}}{{1 - R\left( x \right)}}} \right)$, 给定置信水平1-α, 可得火工品在工作刺激量x的可靠性置信下限:

$ {R_{\rm{L}}}\left( x \right) = {\left( {1 + \frac{{1 - R\left( x \right)}}{{R\left( x \right)}}{\rm{exp}}\left( {{z_{1 - \alpha }}{\sigma _h}} \right)} \right)^{ - 1}} $

(17)

式中, ${\sigma _h} = \sqrt {{\rm{ }}\mathit{\boldsymbol{J}}{\prime _h}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{}}^{ - 1}}{\mathit{\boldsymbol{J}}_h}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)} $, z_1-α为标准正态分布的1-α分位点。

4 实验验证

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

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

根据工程经验可知该火工品的感度服从对数正态分布。通过把对数正态分布变换为正态分布, 综合利用表 1中的三组升降法试验数据进行可靠性分析。取连接函数为Φ(β₁+β₂x), 利用广义线性模型可得参数β的极大似然估计$({{\hat \beta }_1} = - 3.33, {{\hat \beta }_2} = 4.54)$, 进而可得分布参数的估计$(\hat \mu = 0.734, \hat \sigma = 0.22) $。给定置信水平为γ=0.95和工作刺激量水平x=6cm, 结合参数估计的渐近正态性质, 由式(17)可得可靠性下限R_1L=0.9996, 故该火工品的可靠性达到了指标要求。

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

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

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

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

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

根据工程经验可知该火工品的感度服从对数Logistic分布。通过把对数Logistic分布变换为Logistic分布, 综合利用表 3中的三组升降法试验数据进行可靠性分析。取连接函数为$\frac{{{\rm{exp}}\left( {{\beta _1} + {\beta _2}x} \right)}}{{1 + {\rm{exp}}\left( {{\beta _1} + {\beta _2}x} \right)}}$, 利用广义线性模型可得参数β的极大似然估计$({{\hat \beta }_1} = - 130.19, {{\hat \beta }_2} = 22.27)$, 进而可得分布参数的估计$(\hat \mu = 5.846, \hat \sigma = 0.045) $。给定置信水平为γ=0.95和工作刺激量水平x=700 mA, 结合参数估计的渐近正态性质, 由式(17)可得该火工品的可靠性下限R_3L=0.99998, 故该火工品的可靠性达到了指标要求。

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

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

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

5 结论

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