本发明提供一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法。首先，它针对定时截尾加速寿命试验中的小样本数据，以威布尔分布为基础，在保证失效机理一致(形状参数不变)的前提下，构建尺度参数(特征寿命)与温度应力的对数线性加速模型。其次，使用极大似然估计和理论失效数，将定时截尾试验数据转换为定数截尾试验数据。最后，基于信仰推断方法，给出寿命分布参数和加速模型参数的点估计和区间估计，外推给出正常应力水平下可靠性指标的点估计和区间估计，完成固体推进剂可靠性评估。所提出的方法可以有效指导基于定时截尾加速寿命试验数据的固体推进剂可靠性评估。本发明适用于装备剩余寿命预测与可靠性评估等相关。
背景技术：
1、固体推进剂是一种具有特定性能的含能复合材料，它是战略导弹、战术导弹、火箭、空间飞行器等各类武器固体发动机的动力源。固体推进剂是保证发动机正常工作、可靠履行使命的关键基础。在贮存过程中，随着贮存时间的增长以及由于温度变化而产生的热应力，固体推进剂与衬层之间的界面强度会逐渐退化，从而直接影响发动机贮存寿命与可靠性，导致武器无法正常使用。为更加全面了解固体推进剂的性能和寿命特性，需要对固体推进剂在不同温度水平下进行加速寿命试验，对试验数据进行分析，进而评估固体推进剂的寿命和可靠性，为推进剂的设计、应用和改进提供有价值的信息和依据。
2、固体推进剂具有高可靠、长寿命特点，使得正常贮存环境下开展寿命试验所需时间过长、成本过高，有效样本量过少。这时，定时截尾加速寿命试验被广泛使用以加速产品失效过程。实际中，过小的样本量会严重影响参数估计的准确性。因此，亟需解决小样本下、加速寿命试验数据的可靠性评估准确性过低的问题。
3、基于此，本发明提出一种基于定时截尾加速寿命试验数据，综合运用理论失效数、极大似然估计与信仰推断方法的固体推进剂可靠性评估方法。
技术实现思路
1、(1)本发明的目的：
2、本发明针对固体推进剂定时截尾加速寿命数据分析过程中，由固体推进剂高可靠、长寿命导致样本量过低而产生的一系列估计问题，提供一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，具体包含加速寿命试验数据收集、对数线性加速模型建立、定数截尾数据转换及模型参数估计等。考虑高应力对退化过程的加速效应，以威布尔分布为基础，在保证失效机理一致(形状参数不变)的前提下，构建尺度参数(特征寿命)与温度应力的对数线性加速模型，并运用理论失效数、极大似然估计，将定时截尾试验数据转换为定数截尾试验数据，进而基于信仰推断方法分别给出寿命分布参数和加速模型参数的点估计和区间估计，并外推出正常应力水平下可靠性指标的点估计和区间估计，实现固体推进剂可靠性评估。
3、(2)技术方案：基于上述理论和思路，本发明提供一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，具体的实施步骤如下：
4、步骤一：收集加速寿命试验数据
5、基于产品寿命服从威布尔分布的定时截尾加速寿命试验安排如下：设定加速应力为温度，正常应力水平为室内温度t0。设置k个不同应力水平的定时截尾加速寿命试验，每组试验分别对应一个应力水平ti，且满足如下关系：
6、t0＜t1＜t2＜…＜tk, (1)
7、每组各测试ni个产品，i＝1，2，...，k。设定各组截尾时间为ti，s，记录相应的失效时间ti，j，且满足如下关系：
8、
9、其中，ri为应力水平ti下的产品失效数，0≤ri≤ni；而ti，j为第i个试验水平下第j个试验产品的失效时间，j＝1,2,…,ri。
10、步骤二：构建尺度参数与温度应力的对数线性加速模型
11、考虑到固体推进剂的主要敏感应力为温度，而工程中与温度相关的主要加速模型为阿伦尼斯模型，其形式如下：
12、
13、式中，l表示产品的寿命特征，如平均寿命、中位寿命等，k代表波耳兹曼常数、其值为8.617×10-5ev/℃，a代表模型参数，e代表激活能。因此，对此模型两边作对数变换，将阿伦尼斯加速模型转换为对数线性形式，
14、
15、其中，a＝ln a，为待定参数。由式(4)可以看出，寿命特征的对数形式与开尔文温度的倒数呈线性关系。
16、因此，在每个应力水平ti下，假设产品寿命xi服从weibull分布wb(θi，β)，其尺度参数为θi(θi＞0)，i＝0，1，...，k，形状参数为β(β>0)。其概率密度函数和分布函数为：
17、
18、和
19、
20、根据工程经验，在保证失效机理一致(形状参数β不变)的前提下，加速应力仅会改变产品寿命特征(即尺度参数θi)，因此，将加速应力对尺度参数θi的影响表示为：
21、ln(θi)＝α0+α1xi, (7)
22、其中，α0，α1为待估未知参数，而为第i试验应力(温度)水平倒数的归一化取值。因此，在应力水平ti下，产品寿命的概率密度函数和分布函数分别为：
23、
24、和
25、
26、步骤三：基于理论失效数的定数截尾数据转换
27、对于定时截尾加速寿命试验数据，似然函数可表示为：
28、
29、其中，θi＝exp(α0+α1xi)。接着，基于上述似然函数，采用极大似然估计，对产品寿命分布fi(t)进行初步估计。
30、对于某组定时截尾寿命试验，截尾时间为ti，s，定义分布函数fi(t)与经验分布函数差异的统计量di为：
31、
32、在寿命试验中，抽取样品数为n，当假设h0：f(t)＝fi(t)为真时，从理论上，到截尾时间ti，s时，应该有nfi(ti，s)个产品失效，记，
33、ri＝nfi(t0), (12)
34、称为理论截尾数或失效数。
35、本步骤通过极大似然估计获得分布参数初步估计值，进而估计理论分布函数。基于公式(12)计算定时截尾时间的等效失效数，进而将定时截尾试验数据转换为定数截尾试验数据。
36、步骤四：基于信仰推断方法的可靠性评估
37、由于产品寿命服从威布尔分布，根据威布尔分布与指数分布之间的关系，将威布尔分布失效时间转化为指数分布失效时间的形式，即
38、(ti，j/θi)β～exp(1), (13)
39、其中，ti，j为第i个试验水平下第j个试验产品的失效时间。进一步，对于指数分布，由于两倍总试验时间和平均寿命之比服从分布χ2(2ri)，因此，
40、
41、即
42、
43、令
44、
45、则
46、
47、根据对数伽马分布的性质，si的对数变换满足：
48、e[log(si)-βlog(θi)]＝ψ(ri), (18)
49、var[log(si)-βlog(θi)]＝ψ′(ri), (19)
50、其中，μ(x)＝dlog(γ(x))/dx，ψ′(x)＝d2log(γ(x))/dx2。
51、因此，考虑以下回归模型：
52、e(ui/β)＝log(θi)＝α0+α1xi, (20)
53、
54、其中，ui＝log(si)-ψ(ri)。
55、根据高斯-马尔可夫定理，α0，α1的无偏估计分别由下式给出：
56、
57、
58、其中，
59、对于每一个i＝1，2，…，k，j＝1，2，…，ri，ri为第i个应力水平下的理论失效个数，令
60、
61、对于i＝1，2，…，k，有，
62、
63、且wi(β)是一个关于β严格单调的函数。
64、另一方面，注意到w1(β)，w2(β)，…，wk(β)是独立的，由此，定义：
65、
66、可以证明
67、进一步，根据逆变换方法，可得，
68、
69、同时，由于上述等式左边是关于β的严格单调递增函数，因此，求解上式便可以得到参数β的点估计值用估计值代入式(22)和(23)，可得参数α0，α1的点估计值
70、此时，对于任意的0＜γ＜1有，
71、
72、为形状参数β的置信水平为γ的区间估计。其中，是自由度为v的χ2分布的γ分位数；而w-1(w)是等式w(β)＝w，w＞0的解。
73、接着，给出其他参数以及在正常应力t0下的一些重要指标(平均寿命、寿命分位数以及可靠度)的广义置信区间估计。
74、令
75、
76、
77、其中，xi＝1/ti。进一步，可以推导得到，
78、
79、
80、其中，从式(31)和(32)可以明显看出：v1和v2的分布不依赖于任何未知参数。因此，可以将v1和v2作为枢轴量。
81、此外，w(β)是β的严格单调递增函数。令将w(β)＝w的唯一解g(w，t)代替β。此外，从式(29)和(30)中，可以得到，
82、
83、
84、在α0和α1表达式中，将β替换为g(w，t)和代入式(33)和(34)中，得到参数α0，α1的广义枢轴量：
85、
86、
87、其中，
88、令y0，γ和y1，γ分别为y0和y1的γ分位数，那么可以得到α0和α1的置信度为γ的置信区间分别为[y0，(1-γ)/2，y0，1-(1-γ)/2]和[y1，(1-γ)/2，y1，1-(1-γ)/2]。由式(35)和(36)获得的y0，y1的γ分位数，可以通过下述蒙特卡洛仿真算法得到：
89、①对于一组给定的样本量与观测到的失效时间，随机且独立地生成t1～χ2(2r1)，…，tk～χ2(2rk)，由此可以根据w(β)＝w计算g(w，t)＝w-1(w)的值；根据式(31)～(32)计算v1，v2的值。
90、②根据式(35)和(36)，可以求解得到y0，y1的值。
91、③重复步骤1和步骤2，重复次数b≥10000。此时，可以得到b组y0，y1的值。
92、④分别将所有的y0，y1的所有实现值从小到大排序：y0，1＜y0，2＜…＜y0，b，y1，1＜y1，2…＜y1，b，从而得到y0，y1的γ分位数。
93、通过以上方法，可得加速模型中的未知参数α0，α1的广义区间估计。进一步，由于在威布尔分布下，正常应力条件t0下的平均寿命μ、p分位点tp以及可靠度r(t0)分别由以下公式给出，
94、μ＝θ0γ(1+1/β), (37)
95、
96、r(t0)＝exp[-(t0/θ0)β], (39)
97、其中，θ0＝exp(α0+α1x0)，x0＝1/t0。根据上述公式，将对应参数的点估计代入可得平均寿命μ、p分位点tp以及可靠度r(t0)的点估计。另外，与推导α0，α1的广义枢轴量相似，同样可以得到μ，tp和r(t0)的广义枢轴量y2，y3，y4：
98、
99、
100、
101、其中，y2，y3，y4的γ分位数均可通过蒙特卡洛算法获得，进而可得平均寿命μ、p分位点tp以及可靠度r(t0)的区间估计结果。
102、(3)优点和功效：
103、①本发明针对固体推进剂定时截尾加速寿命试验数据分析过程中，固体推进剂高可靠、长寿命、有效样本量过低而产生的一系列估计问题，运用极大似然估计、理论失效数方法以及信仰推断方法进行寿命分布的参数估计，分别给出尺度参数和加速模型参数的点估计和区间估计，并外推出正常应力水平下可靠性指标的置信区间，从而有效支持固体推进剂的贮存管理决策。
104、②本发明提出的方法计算简便，容易实现，且更加符合工程实际，方便工程技术人员掌握使用，方法科学合理，便于应用推广。
1.一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：
2.根据权利要求1所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤二中，在每个应力水平ti下，设产品寿命xi服从weibull分布wb(θi,β)，其尺度参数为θi(θi>0)，i＝0,1,…,k，形状参数为β，β>0；概率密度函数和分布函数为：
3.根据权利要求1或2所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤二中，将加速应力对尺度参数θi的影响表示为：
4.根据权利要求1所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤三中，对于某组定时截尾寿命试验，截尾时间为ti,s，定义分布函数fi(t)与经验分布函数差异的统计量di为：
5.根据权利要求1所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，由于产品寿命服从威布尔分布，根据威布尔分布与指数分布之间的关系，将威布尔分布失效时间转化为指数分布失效时间的形式，即
6.根据权利要求5所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，根据对数伽马分布的性质，si的对数变换满足：
7.根据权利要求6所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，根据高斯-马尔可夫定理，α0,α1的无偏估计分别由下式给出：
8.根据权利要求7所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，接着，给出参数以及在正常应力t0下的重要指标的广义置信区间估计；
9.根据权利要求8所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，令y0,γ和y1,γ分别为y0和y1的γ分位数，得到α0和α1的置信度为γ的置信区间分别为[y0,(1-γ)/2,y0,1-(1-γ)/2]和[y1,(1-γ)/2,y1,1-(1-γ)/2]；由式(35)和(36)获得的y0,y1的γ分位数，通过下述蒙特卡洛仿真算法得到：
10.根据权利要求9所述的一种基于定时截尾加速寿命试验数据的固体推进剂可靠性评估方法，其特征在于：在步骤四中，正常应力条件t0下的平均寿命μ、p分位点tp以及可靠度r(t0)分别由以下公式给出，
