CN112949209B - 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法 - Google Patents

一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法 Download PDF

Info

Publication number
CN112949209B
CN112949209B CN202110330699.8A CN202110330699A CN112949209B CN 112949209 B CN112949209 B CN 112949209B CN 202110330699 A CN202110330699 A CN 202110330699A CN 112949209 B CN112949209 B CN 112949209B
Authority
CN
China
Prior art keywords
degradation
distribution
parameter
representing
stress
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110330699.8A
Other languages
English (en)
Other versions
CN112949209A (zh
Inventor
马小兵
闫冰心
王晗
周堃
吴建国
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
No 59 Research Institute of China Ordnance Industry
Original Assignee
Beihang University
No 59 Research Institute of China Ordnance Industry
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University, No 59 Research Institute of China Ordnance Industry filed Critical Beihang University
Priority to CN202110330699.8A priority Critical patent/CN112949209B/zh
Publication of CN112949209A publication Critical patent/CN112949209A/zh
Application granted granted Critical
Publication of CN112949209B publication Critical patent/CN112949209B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/04—Ageing analysis or optimisation against ageing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computing Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种退化速率‑波动联合更新的弹用密封橡胶贮存寿命评估方法，其步骤如下：步骤一：建立随机退化过程模型：步骤二：确定模型参数的先验分布；步骤三：确定模型参数的后验分布；步骤四：贮存寿命融合评估方法对比分析；通过以上步骤，本发明针对退化型高可靠、长寿命产品自然与加速退化试验数据，根据随机过程理论与加速模型建立产品加速退化模型，运用极大似然估计方法对模型参数进行估计，进而将加速数据模型参数等效到自然模型参数，进而确定模型参数先验分布，运用贝叶斯理论基于自然退化数据更新模型参数的后验分布，最终保证了贮存寿命估计的准确性。

Description

一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估 方法
技术领域
本发明涉及一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，它是一种基于随机过程理论和贝叶斯理论的退化型产品贮存寿命评估方法。它针对退化型产品自然与加速退化试验数据，通过建立维纳过程(“维纳过程”是指独立正态增量过程，也称作布朗运动过程)随机退化模型，对加速退化数据的分布参数进行极大似然估计，并结合阿伦尼斯加速模型(“阿伦尼斯加速模型”是指描述产品寿命和温度应力之间关系的物理模型)，将加速退化数据模型参数等效到自然退化模型参数，得到模型参数的先验信息，并运用贝叶斯理论(“贝叶斯理论”是指根据不确定性信息对事件作出主观推理和决策的理论)更新模型参数的后验分布，基于随机游走Metropolis方法获取参数的后验分布，最终评估产品的期望贮存寿命，最后将本发明方法与基于自然数据的橡胶寿命评估方法进行对比。适用于橡胶构件自然与加速退化试验中贮存寿命融合评估等领域，本发明属于可靠性工程技术领域。
背景技术
硅橡胶作为一种典型的密封材料，在机械系统中得到了广泛的应用，其性能直接影响系统的可靠性和安全性。对于此类退化型失效产品，通过自然退化试验测量产品随时间的性能变化，给定产品失效阈值，可以预测出产品的贮存寿命。然而，在产品使用过程中，积累得到的退化数据量较小，无法准确预测其寿命。加速试验作为高可靠、长寿命产品的重要试验之一，在获取产品寿命及退化信息中发挥着关键作用。加速退化试验可以有效弥补自然退化数据量少而导致的寿命预测不准确问题。基于贝叶斯理论，可以将自然退化数据与加速退化数据有效融合，传统贝叶斯信息融合方法(“贝叶斯信息融合方法”是指基于主观假设与融合观测数据更新主观假设的方法)一般采用共轭先验分布，虽然该假设可以带来较好的统计推断效果，但是难免出现假设偏差，导致寿命预测结果不准确。为提高先验分布与实际数据的拟合程度，应该尝试采用多种备选先验分布，而非共轭先验分布会带来后验分布不具备解析表达式的问题，需要采用仿真方法有效获取后验分布。因此，需要提出一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，准确评估产品寿命。
发明内容
(1)本发明的目的：针对高可靠、长寿命航天产品密封构件橡胶在加速寿命试验中具备较小失效样本，导致其贮存寿命评估无法准确进行的问题，提出一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，它是一种基于随机过程理论和贝叶斯理论的退化型橡胶寿命评估方法。通过建立维纳过程随机退化模型，对加速退化数据的分布参数进行极大似然估计，并结合阿伦尼斯加速模型，将加速退化数据模型参数等效到自然退化模型参数，得到模型参数的先验信息，并运用贝叶斯理论更新模型参数的后验分布，最终评估产品的贮存可靠寿命。
(2)技术方案：
本发明需要建立如下基本设置：
设置1：产品先验信息为加速退化数据，可以利用维纳过程描述其退化，且模型参数具有非共轭先验分布；
设置2：加速退化试验的应力为温度，加载方式为恒定应力，加速模型为阿伦尼斯模型(“阿伦尼斯加速模型”是指描述产品寿命和温度应力之间关系的物理模型)；
基于以上假设，本发明提出一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，通过如下步骤实现：
步骤一：建立随机退化过程模型：
产品在时刻t的性能退化可以表示为X(t)，若随机过程X(t)为维纳过程，则
X(t)＝μt+σB(t) (1)
其中，漂移参数μ表示退化速率，波动参数σ表示退化波动，B(t)表示标准布朗运动(“标准布朗运动”是指具备独立的零均值正态增量的随机过程)；
则产品退化性能首次达到失效阈值的时间为其寿命L，即
L＝inf{t|X(t)≥ω} (2)
其中，ω表示给定的失效阈值，L表示产品寿命，X(t)表示产品退化过程，inf{·}表示下确界；
产品寿命的概率密度函数f(l)如下：
Figure BDA0002994453780000021
其中，l表示产品寿命，ω表示给定的失效阈值，μ表示退化速率，σ表示退化波动；
产品寿命的累计分布函数F(l)如下：
Figure BDA0002994453780000022
其中，Φ(·)表示标准正态分布的分布函数，l表示产品寿命，ω表示给定的失效阈值，μ表示退化速率，σ表示退化波动；
产品寿命的期望如下：
Figure BDA0002994453780000031
其中，
Figure BDA0002994453780000032
表示期望算子，
Figure BDA0002994453780000033
表示产品寿命的期望，ω表示给定的失效阈值，μ表示退化速率；
步骤二：确定模型参数的先验分布；
Ⅰ加速退化试验数据参数估计；
针对恒定加速应力下退化试验数据，T0表示正常应力水平，Ti表示加速应力水平，其中i＝1,2,…,n，xi,j,k表示应力Ti下第j个样本在时刻k的退化观测值，其中j＝1,2,…,m，k＝1,2,…,l；根据维纳过程的性质，样本退化增量服从正态分布
Figure BDA0002994453780000034
可得第j个样本的似然函数为：
Figure BDA0002994453780000035
其中，ΔXi,j,k＝X(ti,j,k)-X(ti,j,k-1)表示产品退化增量，Δti,j,k＝ti,j,k-ti,j,k-1表示时间增量，μi,j表示第i个应力下第j个样本的退化速率，σi,j表示第i个应力下第j个样本的退化波动，
Figure BDA0002994453780000036
表示似然函数；
将似然函数对参数求倒数，并使得倒数等于零，可以得到参数的极大似然估计值为：
Figure BDA0002994453780000037
Figure BDA0002994453780000038
其中，
Figure BDA0002994453780000039
和
Figure BDA00029944537800000310
分别表示μi,j和σi,j的估计值，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μi,j表示第i个应力下第j个样本的退化速率，σi,j表示第i个应力下第j个样本的退化波动；
Ⅱ加速应力到正常应力的退化模型参数等效；
使用阿伦尼斯模型建立维纳退化模型参数μ,σ与加速应力T间的关系，即：
Figure BDA00029944537800000311
Figure BDA0002994453780000041
其中，AF表示加速因子，TS表示高应力下温度，TU表示低应力下温度，βμ和βσ表示加速模型参数；
通过最小二乘估计，可以得到加速模型参数的估计值为：
Figure BDA0002994453780000042
Figure BDA0002994453780000043
其中，Th和Tp表示第h个和第p个加速应力条件下的温度，
Figure BDA0002994453780000044
和
Figure BDA0002994453780000045
表示第h个和第p个加速应力条件下样本退化模型参数的均值，
即
Figure BDA0002994453780000046
同理
Figure BDA0002994453780000047
由式(8)和(9)可将加速应力下的参数估计值(μi,j,σi,j)转化为常规应力下的参数(μ0,j,σ0,j)；
Ⅲ确定模型参数先验分布；
根据常规应力下的模型参数等效值(μ0,j,σ0,j)，基于极大似然准则选择最适合参数(μ,σ)的分布模型，备选分布模型包括正态分布、对数正态分布、伽马分布和逆高斯分布；确定模型参数分布之后，可建立参数分布的似然函数，进而可以得到超参数的极大似然估计值；下面以参数μ0,j为例，给出分布的似然函数；
若参数μ0,j的最优拟合分布为正态分布，即
Figure BDA0002994453780000048
则似然函数为：
Figure BDA0002994453780000049
其中，参数θμ表示正态分布均值，参数
Figure BDA00029944537800000410
表示正态分布方差，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure BDA00029944537800000411
表示似然函数；
若参数μ0,j的最优拟合分布为对数正态分布，即
Figure BDA00029944537800000412
则似然函数为：
Figure BDA0002994453780000051
其中，参数
Figure BDA0002994453780000052
表示对数正态分布均值，参数
Figure BDA0002994453780000053
表示对数正态分布方差，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure BDA0002994453780000054
表示似然函数；
若参数μ0,j的最优拟合分布为伽马分布，即μ0,j～Ga(ρμ,ημ)，则似然函数为：
Figure BDA0002994453780000055
其中，ρμ表示对伽马分布形状参数，ημ表示伽马分布尺度参数，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure BDA00029944537800000515
表示似然函数；
若参数μ0,j的最优拟合分布为逆高斯分布，即
Figure BDA0002994453780000056
则似然函数为：
Figure BDA0002994453780000057
其中，
Figure BDA0002994453780000058
表示逆高斯分布形状参数，κμ表示逆高斯分布尺度参数，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure BDA0002994453780000059
表示似然函数；
为方便表达，将加速退化数据记为DA，将自然退化数据记为DN，两种数据统一记为DO；
步骤三：确定模型参数的后验分布；
I基于贝叶斯理论的后验分布推导
下面以模型参数先验分布
Figure BDA00029944537800000510
σ0,j～Ga(ρσ,ησ)为例，给出模型参数的后验分布；
自然退化数据DN的增量服从正态分布
Figure BDA00029944537800000511
其似然函数为：
Figure BDA00029944537800000512
其中，
Figure BDA00029944537800000513
表示未知参数，联合先验分布
Figure BDA00029944537800000514
表示根据加速退化数据DA得到的退化过程先验信息，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，参数θμ表示正态分布均值，参数
Figure BDA0002994453780000061
表示正态分布方差，
Figure BDA0002994453780000062
表示似然函数；
根据贝叶斯理论(“贝叶斯理论”是指根据不确定性信息对事件作出主观推理和决策的理论)，可以得到模型参数的后验分布为：
Figure BDA0002994453780000063
其中，
Figure BDA0002994453780000064
表示模型参数后验分布，π(θ)表示模型参数先验分布，
Figure BDA0002994453780000065
表示似然函数，参数θμ表示正态分布均值，参数
Figure BDA0002994453780000066
表示正态分布方差，ρσ和ησ表示伽马分布参数；
Ⅱ基于随机游走Metropolis方法的参数后验分布获取
基于上述推导，通过随机游走Metropolis方法(“Metropolis方法”是指以概率来接受新状态的重要性采样)，实现后验信息的获取，进一步基于后验信息给出产品寿命的期望估计；该方法的具体步骤为：
①设置初始值μ(0),σ(0)；
②对于p＝1,…,P，重复如下步骤：
[1]设置μ(p)＝μ(p-1),σ(p)＝σ(p-1)；
[2]从正态分布
Figure BDA0002994453780000067
和
Figure BDA0002994453780000068
中分别生成新的候选参数μ′和σ′，其中，
Figure BDA0002994453780000069
和
Figure BDA00029944537800000610
分别表示参数的方差
Figure BDA00029944537800000611
Figure BDA00029944537800000612
其中：
Figure BDA00029944537800000613
表示参数μ的方差，
Figure BDA00029944537800000614
表示参数lnσ的方差，μl和σl分别表示第l次采样得到的参数；
[3]计算接收概率log α＝min(0,A)，log β＝min(0,B)，其中，
Figure BDA00029944537800000615
Figure BDA00029944537800000616
其中，μ′表示第l次采样所得候选退化速率参数，μp-1表示第l-1次采样所得退化速率参数，σ′表示第l次采样所得候选退化波动参数，σp-1表示第l-1次采样所得退化波动参数；
[4]生成服从均匀分布(0，1)的随机数U，将该随机数与接收概率比较；若U＞α，则接收候选参数μ′，并更新仿真数据μ(p)＝μ′；若U＞β，则接收候选参数σ′，并更新仿真数据σ(p)＝σ′；否则拒绝候选参数，并保持仿真数据不更新；
根据随机游走Metropolis方法生成的参数μ和σ的后验分布记为
Figure BDA0002994453780000071
和
Figure BDA0002994453780000072
其中，μμpos表示参数μpos的后验分布均值，
Figure BDA0002994453780000073
表示参数μpos的后验分布方差，
Figure BDA0002994453780000074
表示参数σpos的后验分布均值，
Figure BDA0002994453780000075
表示参数σpos的后验分布方差；
步骤四：贮存寿命融合评估方法对比分析；
为说明本发明方法的先进性和有效性，将本发明方法与传统基于自然数据的橡胶寿命评估方法进行比较，此处对传统方法可以通过如下步骤实现：
Ⅰ建立随机退化过程模型；
同样的，产品在时刻t的性能退化可以表示为X(t)，如公式(1)所示；
Ⅱ模型参数估计；
针对自然退化数据，T0表示自然应力水平，x0,j,k表示应力T0下第j个样本在时刻k的退化观测值，其中j＝1,2,…,m，k＝1,2,…,l；根据维纳过程的性质，样本退化增量服从正态分布
Figure BDA0002994453780000076
可得第j个样本的极大似然函数为：
Figure BDA0002994453780000077
其中，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，
Figure BDA0002994453780000078
表示似然函数；
将似然函数对参数求倒数，并使得倒数等于零，可以得到参数的极大似然估计值为：
Figure BDA0002994453780000079
Figure BDA00029944537800000710
其中，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，
Figure BDA00029944537800000711
和
Figure BDA00029944537800000712
分别表示μ0,j和σ0,j的估计值；
产品寿命的期望如下：
Figure BDA0002994453780000081
其中，
Figure BDA0002994453780000082
表表示期望算子，L0,j表示正常工作应力下第j个样本的寿命，
Figure BDA0002994453780000083
表示L0,j的期望，ω表示失效阈值，μ0,j表示正常工作应力下第j个样本的退化速率；
Ⅲ寿命评估误差计算；
传统方法寿命评估结果的相对误差er1可表示为：
Figure BDA0002994453780000084
其中，L1表示传统方法寿命评估结果，LA表示产品自然条件下真实寿命，er1表示传统方法寿命评估结果的相对误差；
本发明方法寿命评估结果的相对误差er2可表示为：
Figure BDA0002994453780000085
其中，L2表示本发明方法寿命评估结果，LA表示产品自然条件下真实寿命；
寿命评估误差越小，表示寿命评估结果越准确；
通过以上步骤，本发明针对退化型高可靠、长寿命产品自然与加速退化试验数据，根据随机过程理论与加速模型建立产品加速退化模型，运用极大似然估计方法对模型参数进行估计，进而将加速数据模型参数等效到自然模型参数，进而确定模型参数先验分布，运用贝叶斯理论基于自然退化数据更新模型参数的后验分布，最终保证了贮存寿命估计的准确性。
(3)优点和功效：本发明是一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，其优点是：
①本方法针对退化型高可靠、长寿命产品自然与加速退化试验数据，根据随机过程理论与加速模型建立产品加速退化模型，运用极大似然估计方法对模型参数进行估计，进而将加速数据模型参数等效到自然模型参数，进而确定模型参数先验分布，运用贝叶斯理论基于自然退化数据更新模型参数的后验分布，最终保证了贮存寿命估计的准确性；
②本发明的方法无需给出参数初值，可操作性较强；
③本发明较共轭先验分布对退化模型参数分布形式的要求较低，评估结果更佳准确，体现了使用随机游走Metropolis算法进行后验分布获取的优越性；
④本评估方法科学，工艺性好，具有广阔推广应用价值。
附图说明
图1本发明所述方法流程图；
图2本发明3S-60橡胶加速老化试验数据示意图；
图3本发明3S-60橡胶自然贮存数据示意图；
图4本发明模型参数后验估计值的仿真迭代轨迹图示意图。
具体实施方式
下面结合实际案例对本发明做详细解释。
硅橡胶作为一种典型的密封材料，在机械系统中得到了广泛的应用，其性能直接影响系统的可靠性和安全性。高温下，热老化是3S-60硅橡胶一种关键的失效模式，其退化特征为压缩永久变形率。为研究3S-60硅橡胶的热老化过程，对其进行加速老化试验，并测量其压缩永久变形率。本数据集共有45个样本，分为5组，每组9个试样在不同温度水平下进行测试，并测量橡胶的压缩永久变形率。确定3S-60硅橡胶的实验室加速热老化试验的温度水平为150℃、140℃、130℃、120℃和110℃。3S-60硅橡胶在正常应力下工作的温度水平为56℃，失效阈值为0.3，即压缩永久变形率大于30％时橡胶失效；3S-60橡胶加速老化试验数据如图2所示，3S-60橡胶自然贮存数据如图3所示；
本发明一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，如图1所示，通过如下步骤实现：
步骤一：建立随机退化过程模型，如式子(1)所示；
步骤二：确定模型参数的先验分布；
Ⅰ根据步骤二，对于3S-60橡胶加速退化试验数据，得到参数的极大似然估计值如表2所示：
表2基于加速退化试验数据的模型参数极大似然估计结果
Figure BDA0002994453780000091
Figure BDA0002994453780000101
Ⅱ加速应力到正常应力的退化模型参数等效；
建立阿伦尼斯模型，分别描述退化模型参数μ，σ与加速应力T间的关系，通过最小二乘估计，可以得到加速模型参数的估计值为：
βμ＝-9831.6 (29)
βσ＝-6761.5 (30)
将加速应力下的参数估计值(μij,σij)转化为常规应力下的参数(μ0j,σ0j)，转化后的参数如表3所示：
表3加速应力参数等效为常规应力参数结果
Figure BDA0002994453780000102
Ⅲ确定模型参数先验分布；
根据加速应力参数等效后得到的常规应力参数(μ0j,σ0j),j＝1,2,…,30，基于拟合优度选择最适合参数(μ,σ)的分布模型，模型参数的拟合结果，即参数的对数似然函数如表4所示:
表4模型参数分布拟合结果
正态分布 对数正态分布 伽马分布 逆高斯分布
μ 531.581 528.113 529.516 528.046
σ 438.509 438.877 439.077 438.956
根据极大似然准则可知，参数μ的最优拟合分布为正态分布，即
Figure BDA0002994453780000103
根据极大似然估计，可以得到超参数的估计值为：
θμ＝0.999×10-5,λμ＝1.813×10-6 (31)
同理，参数σ的最优拟合分布为伽马分布，即
Figure BDA0002994453780000111
超参数的估计值为：
Figure BDA0002994453780000112
步骤三：确定模型参数的后验分布。
I基于贝叶斯理论的后验分布推导
根据加速退化数据DA得到的退化模型参数先验分布为
Figure BDA0002994453780000113
σ～Ga(17.8083,3.38×10-6)，则3S-60橡胶自然退化模型参数的联合先验分布记为
π(θ)＝π(0.999×10-5,(1.813×10-6)2,17.8083,3.38×10-6)。
根据贝叶斯原理，可以得到模型参数的后验分布为：
Figure BDA0002994453780000114
Ⅱ基于随机游走Metropolis算法的参数后验分布获取
基于上述推导，基于随机游走Metropolis算法实现后验信息的获取，进一步基于后验信息给出产品寿命的期望估计。仿真次数设置为1000，基于随机游走Metropolis方法可以获取得到参数μ后验估计结果如表5所示，得到参数σ的后验估计结果如表6所示：
表6模型参数后验估计结果
参数 均值 标准差 MC误差 2.5％分位点 中位数 97.5％分位点 仿真次数
μ<sub>pos</sub> 1.21×10<sup>-5</sup> 1.80×10<sup>-6</sup> 5.11×10<sup>-8</sup> 8.78×10<sup>-6</sup> 1.21×10<sup>-5</sup> 1.56×10<sup>-5</sup> 1000
σ<sub>pos</sub> 3.48×10<sup>-4</sup> 2.13×10<sup>-5</sup> 6.62×10<sup>-7</sup> 3.10×10<sup>-4</sup> 3.47×10<sup>-4</sup> 3.94×10<sup>-4</sup> 1000
图4给出了生成参数μ和参数σ后验估计值的仿真迭代轨迹图，结果说明参数后验估计值收敛性较好。
步骤四：贮存寿命融合评估方法对比分析。
为说明本发明方法的先进性和有效性，将本发明方法与传统基于自然数据的橡胶寿命评估方法进行比较，根据步骤四，表7给出了产品贮存寿命评估结果，其中本发明方法期望寿命相对误差小于传统方法的相对误差。
表7产品贮存寿命评估结果
评估方法 真实寿命/天 期望寿命/天 期望寿命相对误差
本发明方法 660 826 25.15％
传统方法 660 918 39.09％
结果表明，采用本发明方法能够合理融合自然与加速退化试验数据，并准确评估产品贮存寿命，达到预期目的。
综上所述，本发明涉及一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，它针退化型产品自然与加速退化试验数据融合问题，通过建立维纳过程随机退化模型，对加速退化数据的分布参数进行极大似然估计，并结合阿伦尼斯加速模型，将加速退化数据模型参数等效到自然退化模型参数，得到模型参数的先验信息，并运用贝叶斯理论更新模型参数的后验分布，最终评估产品的贮存可靠寿命。该方法的具体步骤是：首先建立随即退化过程模型；然后估计加速退化试验数据模型参数，建立加速应力到正常应力的退化模型参数等效，基于极大似然准则确定模型参数先验分布；其次基于贝叶斯理论推导参数后验分布，基于随机游走Metropolis算法获取参数的后验分布；最后给出产品贮存寿命期望，分析期望寿命与真实寿命间的相对误差。本发明适用于具备自然与加速试验数据的产品贮存寿命融合评估，具备较强的操作性。

Claims (1)

1.一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法，其特征在于：
设置1：产品先验信息为加速退化数据，利用维纳过程描述其退化，且模型参数具有非共轭先验分布；
设置2：加速退化试验的应力为温度，加载方式为恒定应力，加速模型为阿伦尼斯模型，阿伦尼斯加速模型是指描述产品寿命和温度应力之间关系的物理模型；
具体步骤如下：
步骤一：建立随机退化过程模型：
产品在时刻t的性能退化表示为X(t)，若随机过程X(t)为维纳过程，则
X(t)＝μt+σB(t) (1)
其中，漂移参数μ表示退化速率，波动参数σ表示退化波动，B(t)表示标准布朗运动；
则产品退化性能首次达到失效阈值的时间为其寿命L，即
L＝inf{t|X(t)≥ω} (2)
其中，ω表示给定的失效阈值，L表示产品寿命，X(t)表示产品退化过程，inf{·}表示下确界；
产品寿命的概率密度函数f(l)如下：
Figure FDA0002994453770000011
其中，l表示产品寿命，ω表示给定的失效阈值，μ表示退化速率，σ表示退化波动；产品寿命的累计分布函数F(l)如下：
Figure FDA0002994453770000012
其中，Φ(·)表示标准正态分布的分布函数，l表示产品寿命，ω表示给定的失效阈值，μ表示退化速率，σ表示退化波动；
产品寿命的期望如下：
Figure FDA0002994453770000013
其中，
Figure FDA0002994453770000021
表示期望算子，
Figure FDA0002994453770000022
表示产品寿命的期望，ω表示给定的失效阈值，μ表示退化速率；
步骤二：确定模型参数的先验分布；
2.1加速退化试验数据参数估计；
针对恒定加速应力下退化试验数据，T0表示正常应力水平，Ti表示加速应力水平，其中i＝1,2,…,n，xi,j,k表示应力Ti下第j个样本在时刻k的退化观测值，其中j＝1,2,…,m，k＝1,2,…,l；根据维纳过程的性质，样本退化增量服从正态分布
Figure FDA0002994453770000023
得到第j个样本的似然函数为：
Figure FDA0002994453770000024
其中，ΔXi,j,k＝X(ti,j,k)-X(ti,j,k-1)表示产品退化增量，Δti,j,k＝ti,j,k-ti,j,k-1表示时间增量，μi,j表示第i个应力下第j个样本的退化速率，σi,j表示第i个应力下第j个样本的退化波动，l(μ,σ)表示似然函数；
将似然函数对参数求倒数，并使得倒数等于零，得到参数的极大似然估计值为：
Figure FDA0002994453770000025
Figure FDA0002994453770000026
其中，
Figure FDA0002994453770000027
和
Figure FDA0002994453770000028
分别表示μi,j和σi,j的估计值，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μi,j表示第i个应力下第j个样本的退化速率，σi,j表示第i个应力下第j个样本的退化波动；
2.2加速应力到正常应力的退化模型参数等效；
使用阿伦尼斯模型建立维纳退化模型参数μ,σ与加速应力T间的关系，即：
Figure FDA0002994453770000029
Figure FDA0002994453770000031
其中，AF表示加速因子，TS表示高应力下温度，TU表示低应力下温度，βμ和βσ表示加速模型参数；
通过最小二乘估计，得到加速模型参数的估计值为：
Figure FDA0002994453770000032
Figure FDA0002994453770000033
其中，Th和Tp表示第h个和第p个加速应力条件下的温度，
Figure FDA0002994453770000034
和
Figure FDA0002994453770000035
表示第h个和第p个加速应力条件下样本退化模型参数的均值，
即
Figure FDA0002994453770000036
同理
Figure FDA0002994453770000037
由式(8)和(9)将加速应力下的参数估计值(μi,j,σi,j)转化为常规应力下的参数(μ0,j,σ0,j)；
2.3确定模型参数先验分布；
根据常规应力下的模型参数等效值(μ0,j,σ0,j)，基于极大似然准则选择最适合参数(μ,σ)的分布模型，备选分布模型包括正态分布、对数正态分布、伽马分布和逆高斯分布；确定模型参数分布之后，建立参数分布的似然函数，进而得到超参数的极大似然估计值；
若参数μ0,j的最优拟合分布为正态分布，即
Figure FDA0002994453770000038
则似然函数为：
Figure FDA0002994453770000039
其中，参数θμ表示正态分布均值，参数
Figure FDA00029944537700000310
表示正态分布方差，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure FDA00029944537700000311
表示似然函数；
若参数μ0,j的最优拟合分布为对数正态分布，即
Figure FDA00029944537700000312
则似然函数为：
Figure FDA0002994453770000041
其中，参数
Figure FDA0002994453770000042
表示对数正态分布均值，参数
Figure FDA0002994453770000043
表示对数正态分布方差，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure FDA0002994453770000044
表示似然函数；
若参数μ0,j的最优拟合分布为伽马分布，即μ0,j～Ga(ρμ,ημ)，则似然函数为：
Figure FDA0002994453770000045
其中，ρμ表示对伽马分布形状参数，ημ表示伽马分布尺度参数，μ0,j表示正常工作应力下第j个样本的退化速率，l(μ0,j|ρμ,ημ)表示似然函数；
若参数μ0,j的最优拟合分布为逆高斯分布，即
Figure FDA0002994453770000046
则似然函数为：
Figure FDA0002994453770000047
其中，
Figure FDA0002994453770000048
表示逆高斯分布形状参数，κμ表示逆高斯分布尺度参数，μ0,j表示正常工作应力下第j个样本的退化速率，
Figure FDA0002994453770000049
表示似然函数；
为方便表达，将加速退化数据记为DA，将自然退化数据记为DN，两种数据统一记为DO；
步骤三：确定模型参数的后验分布；
3.1基于贝叶斯理论的后验分布推导
模型参数先验分布为
Figure FDA00029944537700000410
σ0,j～Ga(ρσ,ησ)；
自然退化数据DN的增量服从正态分布
Figure FDA00029944537700000411
其似然函数为：
Figure FDA00029944537700000412
其中，
Figure FDA0002994453770000051
表示未知参数，联合先验分布
Figure FDA0002994453770000052
表示根据加速退化数据DA得到的退化过程先验信息，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，参数θμ表示正态分布均值，参数
Figure FDA0002994453770000053
表示正态分布方差，l(DN,μ,σ|θ)表示似然函数；
根据贝叶斯理论得到模型参数的后验分布为：
Figure FDA0002994453770000054
其中，
Figure FDA0002994453770000055
表示模型参数后验分布，π(θ)表示模型参数先验分布，l(DN,μ,σ|θ)表示似然函数，参数θμ表示正态分布均值，参数
Figure FDA0002994453770000056
表示正态分布方差，ρσ和ησ表示伽马分布参数；
3.2基于随机游走Metropolis方法的参数后验分布获取
基于上述推导，通过随机游走Metropolis方法，实现后验信息的获取，进一步基于后验信息给出产品寿命的期望估计；该方法的具体步骤为：
设置初始值μ(0),σ(0)；
对于p＝1,…,P，重复如下步骤：
3.21设置μ(p)＝μ(p-1),σ(p)＝σ(p-1)；
3.22从正态分布
Figure FDA0002994453770000057
和
Figure FDA0002994453770000058
中分别生成新的候选参数μ′和σ′，其中，
Figure FDA0002994453770000059
和
Figure FDA00029944537700000510
分别表示参数的方差
Figure FDA00029944537700000511
Figure FDA00029944537700000512
其中：
Figure FDA00029944537700000513
表示参数μ的方差，
Figure FDA00029944537700000514
表示参数lnσ的方差，μl和σl分别表示第l次采样得到的参数；
3.23计算接收概率logα＝min(0,A)，logβ＝min(0,B)，其中，
Figure FDA00029944537700000515
Figure FDA0002994453770000061
其中，μ′表示第l次采样所得候选退化速率参数，μp-1表示第l-1次采样所得退化速率参数，σ′表示第l次采样所得候选退化波动参数，σp-1表示第l-1次采样所得退化波动参数；
3.24生成服从均匀分布(0，1)的随机数U，将该随机数与接收概率比较；若U＞α，则接收候选参数μ′，并更新仿真数据μ(p)＝μ′；若U＞β，则接收候选参数σ′，并更新仿真数据σ(p)＝σ′；否则拒绝候选参数，并保持仿真数据不更新；
根据随机游走Metropolis方法生成的参数μ和σ的后验分布记为
Figure FDA0002994453770000062
和
Figure FDA0002994453770000063
其中，
Figure FDA0002994453770000064
表示参数μpos的后验分布均值，
Figure FDA0002994453770000065
表示参数μpos的后验分布方差，
Figure FDA0002994453770000066
表示参数σpos的后验分布均值，
Figure FDA0002994453770000067
表示参数σpos的后验分布方差；
步骤四：贮存寿命融合评估方法对比分析；具体步骤为：
4.1建立随机退化过程模型；
产品在时刻t的性能退化表示为X(t)，如公式(1)所示；
4.2模型参数估计；
针对自然退化数据，T0表示自然应力水平，x0,j,k表示应力T0下第j个样本在时刻k的退化观测值，其中j＝1,2,…,m，k＝1,2,…,l；根据维纳过程的性质，样本退化增量服从正态分布
Figure FDA0002994453770000068
得第j个样本的极大似然函数为：
Figure FDA0002994453770000069
其中，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，l(μ0,j,σ0,j)表示似然函数；
将似然函数对参数求倒数，并使得倒数等于零，得到参数的极大似然估计值为：
Figure FDA00029944537700000610
Figure FDA0002994453770000071
其中，ΔXi,j,k表示产品退化增量，Δti,j,k表示时间增量，μ0,j表示正常工作应力下第j个样本的退化速率，σ0,j表示正常工作应力下第j个样本的退化波动，
Figure FDA0002994453770000072
和
Figure FDA0002994453770000073
分别表示μ0,j和σ0,j的估计值；
产品寿命的期望如下：
Figure FDA0002994453770000074
其中，
Figure FDA0002994453770000075
表示期望算子，L0,j表示正常工作应力下第j个样本的寿命，
Figure FDA0002994453770000076
表示L0,j的期望，ω表示失效阈值，μ0,j表示正常工作应力下第j个样本的退化速率；
4.3寿命评估误差计算；
传统方法寿命评估结果的相对误差er1表示为：
Figure FDA0002994453770000077
其中，L1表示传统方法寿命评估结果，LA表示产品自然条件下真实寿命，er1表示传统方法寿命评估结果的相对误差；
寿命评估结果的相对误差er2表示为：
Figure FDA0002994453770000078
其中，L2表示寿命评估结果，LA表示产品自然条件下真实寿命；寿命评估误差越小，表示寿命评估结果越准确。
CN202110330699.8A 2021-03-26 2021-03-26 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法 Active CN112949209B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110330699.8A CN112949209B (zh) 2021-03-26 2021-03-26 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110330699.8A CN112949209B (zh) 2021-03-26 2021-03-26 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法

Publications (2)

Publication Number Publication Date
CN112949209A CN112949209A (zh) 2021-06-11
CN112949209B true CN112949209B (zh) 2022-05-17

Family

ID=76227090

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110330699.8A Active CN112949209B (zh) 2021-03-26 2021-03-26 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法

Country Status (1)

Country Link
CN (1) CN112949209B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113407896B (zh) * 2021-06-15 2022-10-14 电子科技大学 一种核电阀门加速寿命试验中的最佳应力和样本分配确定方法
CN113468757B (zh) * 2021-07-16 2022-08-30 西南石油大学 一种基于模糊随机理论评估腐蚀天然气管道可靠性的方法
CN113688513A (zh) * 2021-08-17 2021-11-23 中国电力科学研究院有限公司 Opgw光缆的寿命评估方法、系统、设备及存储介质
CN113705112B (zh) * 2021-09-23 2023-06-27 郑州航空工业管理学院 一种恒温器寿命doe重要因子的识别方法
CN114626248B (zh) * 2022-03-30 2023-04-18 北京航空航天大学 一种基于多应力加速退化数据的螺旋弹簧可靠性评估方法
CN114792053B (zh) * 2022-04-24 2022-10-11 哈尔滨工业大学 一种基于初值-速率相关退化模型的可靠性评估方法
CN114818348B (zh) * 2022-05-06 2022-10-11 哈尔滨工业大学 考虑多应力耦合作用对产品退化影响的可靠性评估方法
CN115828648B (zh) * 2023-02-15 2023-06-16 中国科学技术大学 一种火灾后电子封装互连结构寿命的预测方法及装置
CN115876681B (zh) * 2023-03-01 2023-05-23 中南大学 一种用于密封垫的安全度评估方法及测试装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102622473A (zh) * 2012-02-28 2012-08-01 北京航空航天大学 基于贝叶斯理论的步进应力加速退化试验优化设计方法
CN111832151A (zh) * 2020-05-27 2020-10-27 北京航空航天大学云南创新研究院 基于指数时间函数的Wiener加速退化模型构建方法及系统
CN111859658A (zh) * 2020-07-15 2020-10-30 北京强度环境研究所 一种产品贮存寿命与可靠性评估方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108959676B (zh) * 2017-12-22 2019-09-20 北京航空航天大学 一种考虑有效冲击的退化建模与寿命预测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102622473A (zh) * 2012-02-28 2012-08-01 北京航空航天大学 基于贝叶斯理论的步进应力加速退化试验优化设计方法
CN111832151A (zh) * 2020-05-27 2020-10-27 北京航空航天大学云南创新研究院 基于指数时间函数的Wiener加速退化模型构建方法及系统
CN111859658A (zh) * 2020-07-15 2020-10-30 北京强度环境研究所 一种产品贮存寿命与可靠性评估方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A Novel Parameter-Related Wiener Process Model With Its Estimation;Xiaobing Ma;《IEEE》;20201231;全文 *
Mechanism Equivalence in Designing Optimum;HAN WANG 等;《IEEE》;20180228;全文 *
Objective Bayesian analysis for the accelerated;Daojiang He 等;《STATISTICAL THEORY AND RELATED FIELDS》;20180515;全文 *
Remaining Useful Life Prediction Using;HONGYU WANG 等;《IEEE》;20181130;全文 *
基于Copula函数耦合性建模的二元加速退化数据统计分析方法;周源等;《兵器装备工程学报》;20180525(第05期);全文 *
基于多阶段-随机维纳退化过程的产品剩余寿命预测方法;魏高乐等;《科学技术与工程》;20150918(第26期);全文 *

Also Published As

Publication number Publication date
CN112949209A (zh) 2021-06-11

Similar Documents

Publication Publication Date Title
CN112949209B (zh) 一种退化速率-波动联合更新的弹用密封橡胶贮存寿命评估方法
CN110851980B (zh) 一种设备剩余寿命预测方法及系统
CN108959676B (zh) 一种考虑有效冲击的退化建模与寿命预测方法
CN108664700B (zh) 基于不确定数据包络分析的加速退化信息融合建模方法
CN102542155A (zh) 基于加速退化数据的粒子滤波剩余寿命预测方法
CN106570281A (zh) 基于相似产品信息的小子样产品贝叶斯可靠性评估方法
Ismail et al. RETRACTED ARTICLE: Optimal planning of failure-step stress partially accelerated life tests under type-II censoring
CN101710368A (zh) 基于多源退化数据的贝叶斯可靠性综合评估方法
CN107515965A (zh) 一种基于不确定过程的加速退化建模评估方法
CN105426647B (zh) 基于可靠度先验信息融合的冷备系统可靠度估计方法
CN115962797B (zh) 一种基于温度应力下的传感器可靠性测试方法及系统
CN114970157B (zh) 电子产品在电压应力作用下的小样本试验寿命预测方法
CN113486455A (zh) 一种自润滑关节轴承加速退化可靠性评估与寿命预测方法
Wang et al. A new class of mechanism-equivalence-based Wiener process models for reliability analysis
Bressel et al. Fuel cell remaining useful life prediction and uncertainty quantification under an automotive profile
CN114169128B (zh) 一种基于Bayes分析的可靠性强化试验定量评估方法
El-Adll Inference for the two-parameter exponential distribution with generalized order statistics
Anderson-Cook et al. Reliability modeling using both system test and quality assurance data
CN109359375B (zh) 一种数据产品贮存寿命预测方法
CN111400856B (zh) 基于多源数据融合的空间行波管可靠性评估方法
CN115344960A (zh) 一种基于贝叶斯信息融合的涡轮盘可靠性评估方法
CN111859296B (zh) 一种基于装备使用期间的测试性指标评估方法及系统
CN114996928B (zh) 电子产品的温压双应力小样本加速试验寿命预测方法
CN111460638B (zh) 一种考虑个体差异性和测量误差的产品剩余使用寿命预测方法
CN115576721B (zh) 基于不确定性量化的多源数据融合评估试验方法及设备

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant