CN104615866A - 一种基于物理统计模型的寿命预测方法 - Google Patents

Abstract

一种基于物理统计模型的寿命预测方法，它包括以下具体步骤：步骤一：构建产品退化的物理模型；步骤二：将产品退化的物理模型转换成线性模型的形式；步骤三：对物理模型的线性形式进行参数估计；步骤四：根据产品性能参数退化轨迹的特点，选择合适的随机过程，构造物理统计模型；步骤五：选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新；步骤六：对含卡尔曼滤波更新的物理统计模型中的参数进行参数估计；步骤七：求解产品的寿命分布和可靠度函数。本发明在可靠性技术领域里有较好的应用前景。

Description

一种基于物理统计模型的寿命预测方法

技术领域

本发明提出一种基于物理统计模型的寿命预测方法，它是一种基于物理统计结合卡尔曼滤波更新的加速退化试验数据建模方法，属于可靠性技术领域中的加速退化试验数据处理。

背景技术

随着产品设计、材料、制造工艺的发展，产品的可靠性和寿命已经可以达到一个很高的水平。随之而来亟需解决的难题是如何评估这些高可靠性产品的寿命。对于寿命少则几年多则数十年的设备或系统，加速试验在试验技术上为解决这一问题提供了可行的理论支撑。然而，如何对加速试验数据进行合理的处理，得出产品寿命和可靠性较为精确的估计依旧是一个问题。

目前加速试验数据处理主要有两种思路：

1.在加速试验中产品出现失效，利用失效数据对产品寿命符合的概率分布形式进行估计，从而得到产品的寿命分布和可靠度分布函数；

2.在加速试验中产品未出现失效，但是能从产品的性能参数监测到具有退化特征，利用退化数据对产品的退化轨迹进行建模，外推得到产品的寿命。

在第二种思路中，对于退化轨迹进行建模又主要有两种思路：

1.依据产品发生退化的机理建立物理模型，之后利用退化数据对模型中的参数进行估计，得到产品的退化轨迹模型；

2.依据产品退化轨迹的特点，选择合适的统计模型，并利用退化数据对模型中的参数进行估计，最终得到产品的退化轨迹模型。

利用物理模型对退化轨迹进行建模的方法中，对产品寿命和可靠度的预测的结果与建立的物理模型是否准确直接相关。然而在实际工程中，针对产品建立精确的物理模型是较为复杂的。利用统计模型对退化轨迹进行建模的方法中，对产品寿命和可靠度的预测的结果与退化轨迹的特点是否突出有密切联系，不仅需要数据量较大，而且需要数据能够体现出统计规律。一种新的思路是将上述提到的物理方法与统计方法相结合，这种混合建模方法使得能够利用对产品的结构、机理认识建立的物理模型，同时结合统计手段，在建模是能够将样品间个体差异作为随机性的因素包含在模型中。

发明内容

针对加速退化试验的退化数据处理，发明了一种基于物理统计结合卡尔曼滤波更新的加速退化数据建模方法。本发明的为了精确的对加速退化试验得到的退化轨迹进行外推，全面考虑了产品的物理退化机理、产品个体特性、产品退化各阶段的退化轨迹差异性，对于高可靠长寿命的产品，能够在指定失效阈值的前提下进行寿命分布和可靠度函数的预测。

本发明提出的退化过程有依据以下四点假设前提：

1.所研究的产品是具备性能参数退化的产品，而不是单一功能的成败型产品，即产品的退化状态是一系列性能参数采样点，而不是功能正常与功能不正常两点；

2.产品的性能参数退化趋势总体是具有单调性的，即退化轨迹应是逐渐趋近失效阈值；

3.产品的性能参数具有失效阈值，即产品的某项性能参数或输出参数达到性能指定的失效阈值后，产品本身或装备了该产品的系统将不能正常稳定工作；

4.导致产品发生性能退化的应力及机理在退化过程是确定不变的，即产品自退化开始到达失效阈值为止，促使其退化的应力及机理在退化过程中始终促使产品发生退化，并且不存在其它应力及机理仅在退化过程中某一阶段起作用。

基于上述假设，本发明提供的一种基于物理统计模型的寿命预测方法，主要包括以下具体几个步骤：

步骤一：构建产品退化的物理模型。

根据产品的使用信息，确定产品可能经受的环境应力和工作应力，依此对产品退化的物理模型进行调研，并结合对产品结构、工况等方面的分析，确定可能的退化模型。对于单应力退化，通常其模型的形式是y＝f(t,S)，其中y是性能参数退化量，t是时间，S是加速应力。

步骤二：将产品退化的物理模型转换成线性模型的形式。

根据产品退化的物理模型，转换成线性模型形式的具体方法主要有以下三种：

1.若在物理模型中，含有应力水平S的项与含有时间t的项是乘积关系，则可以将两项分别做整体代换构建成线性模型。即若产品的加速退化物理模型形式为y＝Ag(S)h(t)，则可以分别作代换B_S＝g(S)，C_t＝h(t)，将物理模型代换成线性模型y＝AB_SC_t，其中B_S为与应力水平有关的系数，C_t是经过形式转换的时间量。在此模型中，性能参数的退化量与经过形式转换的时间量具有线性关系。

2.若在物理模型中，含有应力水平S的项与含有时间t的项是指数关系，则可以对模型等式两边先取对数，再分别作整体代换构建线性模型。即若产品的加速退化物理模型形式为y＝Ag(S)^h(t)或y＝Ah(t)^g(S)，先对等式两边取对数得ln y＝h(t)ln g(S)+ln A或ln y＝g(S)ln h(t)+ln A。对第一种情况作代换B_S＝ln g(S)，C_t＝h(t)；对第二种情况作代换B_S＝g(S)，C_t＝ln h(t)，都可以得到线性模型形式ln y＝B_SC_t+ln A。在此模型中，性能参数退化量的对数与经过形式转换的时间量具有线性关系。

3.若在物理模型中，含有应力水平S的项与含有时间t的项是其他关系，则可以根据具体情况用包括取对数、整体代换等方法变换为线性模型的形式。

步骤三：对物理模型的线性形式进行参数估计。

本发明中，采用了最小二乘法，利用历史样本数据对变换成线性形式的物理模型进行参数估计。以线性模型y＝AB_SC_t为例说明：

1.整理产品性能参数的退化数据，将试验的真实时间量度按照线性形式中时间转换形式进行转换，从而将退化轨迹中退化量与真实时间的对应关系转换为退化量与经时间转换形式转换后的计算时间的对应关系；

2.计算物理模型在每一时刻求得的退化量与加速退化试验测得的退化量偏差。偏差其中表示在第j个样本在对应应力水平S_j下在第i个性能退化监测点时刻t_ji时测得的性能参数退化量。

3.计算物理模型在每一时刻求得的退化量与加速退化试验测得的退化量偏差的平方和，即对各应力水平下所有样本的偏差进行平方加和。偏差的平方和 $\mathrm{SS}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{\left[e\left(C_{t_{\mathrm{ji}}}\right)\right]}^2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{[{\mathrm{AB}}_{S_j}{C}_{t_{\mathrm{ji}}}-y_{t_{\mathrm{ji}}}]}^2.$

4.利用奈尔德-米德单纯性法求解偏差的平方和的最小值，当偏差的平方和取得最小值时的对应的参数值即为物理模型线性形式的参数估计值。

步骤四：根据产品性能参数退化轨迹的特点，选择合适的随机过程，构造物理统计模型。

在线性模型的基础上，依据产品性能退化的特点，选择合适的随机过程，构造得到的随机过程模型即为最终的物理统计模型。在选择适用于产品退化过程的随机过程时，首先需要对产品退化过程的特点进行分析：

1.产品性能参数的退化量是严格单调递增或递减的，适合选择伽马过程进行建模；

2.产品性能参数的退化量增量是相互独立，且符合正态分布，适合选择布朗运动进行建模；

3.产品性能参数的退化量增量是相互独立，但不符合正态分布，可以尝试选择逆高斯过程。

例如产品的性能参数退化的线性模型形式是y＝AB_SC_t，且退化量增量是相互独立且符合正态分布，则可以构造布朗运动随机过程模型y＝μC_t+σB(C_t)，其中μ为漂移布朗运动的漂移系数，μ＝AB_S，σ为漂移布朗运动的扩散系数，B(C_t)是以经过形式转换的时间量C_t为时间量的标准布朗运动，其服从N(0,C_t)正态分布。并且根据代换关系可以利用线性模型形式中已知的参数估计值计算得到物理统计模型中的参数估计。在上例中，即可由μ＝AB_S计算得到μ的取值。

步骤五：选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新。

选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新，使得物理统计模型能更好的拟合退化轨迹从而更精确的预测。

卡尔曼滤波算法能够利用目标的历史信息，滤去噪声的影响，得到一个较为准确的估计值。以步骤四中的服从漂移布朗运动的物理统计模型y＝μC_t+σB(C_t)为例，在此模型中与样本间个体差异密切相关且与退化轨迹形态密切相关的是漂移系数μ，因此用卡尔曼滤波算法对漂移系数μ进行随时间的更新计算。

假设初始时刻的服从正态分布N(μ₀,P₀)，假设对于漂移系数从到存在的更新过程，其中η～N(0,Q)。则卡尔曼滤波从到的更新过程按下述式进行：

${\widehat{\mathrm{\μ}}}_{t_{i+1}}={\widehat{\mathrm{\μ}}}_{t_i}+P_{t_{i+1}|t_i}(t_{i+1}-t_i){F_{t_{i+1}}}^{-1}({\mathrm{\Δy}}_{t_{i+1}}-{\mathrm{\Δy}}_{t_i}-{\widehat{\mathrm{\μ}}}_{t_i}(t_{i+1}-t_i))$

$P_{t_{i+1}|t_i}=P_{t_i|t_i}+Q$

$F_{t_{i+1}}={(t_{i+1}-t_i)}^2{P}_{t_{i+1}|t_i}+{\mathrm{\σ}}^2(t_{i+1}-t_i)$

$P_{t_{i+1}|t_{i+1}}=P_{t_{i+1}|t_i}-P_{t_{i+1}|t_i}{(t_{i+1}-t_i)}^2{F_{t_{i+1}}}^{-1}{P}_{t_{i+1}|t_i}$

当漂移系统更新到时，对于t_i后的任意时刻t_k，预测退化轨迹的物理统计模型为： $y_{t_k}=y_{t_i}+{\widehat{\mathrm{\μ}}}_{t_i}(t_k-t_i)+\mathrm{\σB}(t_k,t_i).$

步骤六：对含卡尔曼滤波更新的物理统计模型中的参数进行参数估计。

利用产品性能参数的退化数据对含卡尔曼滤波更新物理统计模型中的参数进行参数估计。本发明中在此使用的参数估计方法是极大似然估计法。参数估计方法的主要流程是：

1.根据物理统计模型，确定一服从某一分布形式的统计量，例如漂移布朗运动模型中，退化量的增量即为服从正态分布的统计量；

2.根据物理统计模型，求出该统计量的表达式及其所服从的分布，例如步骤五的物理统计模型中，统计量的表达式即为两时刻点退化量的之差服从正态分布 $N({\mathrm{\μ}}_{t_i}(t_{i+1}-t_i),P_{t_{i+1}|t_i}{(t_{i+1}-t_i)}^2+{\mathrm{\σ}}^2(t_{i+1}-t_i));$

3.列出极大似然函数，续前例，服从均值为方差为的正态分布，则由连续分布的极大似然函数法可得对数极大似然函数为：

$\begin{array}{c}\log \left(L\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}(-\frac{1}{2}\ln \left(2\mathrm{\π}(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))\right))\\ -\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}\left(\frac{{(y_{t_{j(i+1)}}-y_{t_{\mathrm{ji}}}-{\widehat{\mathrm{\μ}}}_{t_{\mathrm{ji}}}(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}^2}{2(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}\right)\end{array}$

其中j表示第j个样本，i表示共m_j个样品的性能参数退化轨迹点中的第i个性能参数退化轨迹点。

4.利用搜索算法对于对数极大似然函数求极值，并求出取得极值时对应的参数值。在本发明中，此处使用粒子群寻优算法。粒子群寻优算法是一种启发式智能搜索算法，其是通过沿搜索得到的解中最优解方向继续迭代搜索，能够快速有效的收敛于极值点。

步骤七：求解产品的寿命分布和可靠度函数。

依据物理统计模型可以求解产品的寿命分布和可靠度函数。运用不同的统计模型构造的物理统计模型其寿命分布和可靠度函数的形式是不同的。以利用漂移布朗运动构造的物理统计模型为例，产品的寿命分布，也即首穿时分布服从逆高斯分布，并可由寿命分布求得可靠度分布。对与寿命分布为L＝f(t)的产品，其可靠度函数为

本发明的优点在于：

1.本发明以产品的物理模型为基础，这是建立退化轨迹模型的理论基础。在后续的模型构建和推导都是以物理模型为基础，从而在每一步构建得到的模型中参数都是物理模型的一个或几个参数的演变，是有与其对应的明确物理意义，保证了模型是确实可信的服从物理理论的，是模型能够得到准确预测结果的根源。

2.本发明在物理模型的基础上，引入了统计学中的随机过程模型。随机过程反应了一系列随机现象中的统计特性。在物理模型的基础上引入随机过程模型，使得物理统计模型既具备了预测退化轨迹准确的退化趋势又能够符合客观实际的体现性能参数在退化过程中的随机波动。此处将物理与统计相结合，是对有限的性能退化数据的充分挖掘利用。

3.本发明将卡尔曼滤波引入物理统计模型，使得物理统计模型中参数能够随退化的进程进行更新，从而使得模型能够更好的拟合各个时刻的退化轨迹，得到更准确的预测结果。

4.本发明提出的物理统计模型能够方便的求解出产品的寿命分布和可靠度函数。借助随机过程的统计性质，物理统计模型能够方便的求解得到寿命分布，亦称作首穿时分布，首达时分布，其是指退化轨迹首次达到指定的失效阈值的时刻。依据寿命分布，可以方便的求出可靠度函数。

1.本发明创新的将产品退化中的退化机理、性能参数退化的随机性、退化各阶段可能具备不同特点这三点同时进行考虑，克服了传统建模方法考虑不全面只针对单一的特性进行建模的缺点。

附图说明

图1是本发明方法流程图；

图2加速度计标度因素置信区间动态趋势图；

图3四号加速度计标度参数的退化轨迹图；

图4四号加速度计的寿命分布图；

图5四号加速度计的可靠度函数图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明一种基于物理统计模型的寿命预测方法，它是一种基于物理统计结合卡尔曼滤波更新的加速退化数据建模方法，首先对本发明方法的意义进行详细说明。

对于高可靠性产品，产品的退化轨迹数据是宝贵的，使用物理模型对退化轨迹进行拟合只能得到确定性模型，不能够体现退化过程中的随机性也无法得到寿命分布；使用统计学的随机过程模型能够体现退化过程中的随机性，但是其模型参数的物理意义无法确定因而不利于判定对寿命影响较大的物理参数，从而无法对改进产品提供明确的方向；在传统的退化轨迹建模中，通常认定退化量或退化量增量是服从某一固定分布，而实际中产品在各个阶段退化的速率是可能有所差异的，其服从的分布随退化的进城是可能发生改变的。本发明结合了物理、统计方法，利用里卡尔曼滤波进行参数的更新，有效的解决了上述方法的缺陷和问题。

以下实施例是按照附图1所示的流程进行实施的，本实施例选取了加速度计标度因数退化试验数据处理作为范例。标度因素是加速度计的关键性能参数，当标度因数退化至失效阈值，加速度计测得的加速度会超出误差允许范围，从而判定加速度计发生了失效。本方法为高可靠长寿命的加速度计标度因素退化轨迹建模提供了行之有效的建模方法。

首先检查加速度计的标度因数退化轨迹，确认其是否满足本发明的假设前提。由加速退化试验数据可知，加速度计的标度因数是具备性能参数退化的，其退化趋势逐渐上升，根据设计指标有明确的失效阈值，在试验条件下，温度是促使其发生退化的加速应力，以上四点事实表明加速度计的退化符合本发明的四点假设前提。

其具体实施步骤如下：

步骤一：构建产品退化的物理模型。

根据对加速度计的分析和调研，可以确定加速度计标度因数的退化模型是其中，a与b是与加速度计中某材料相关的两个模型参数，T是温度应力，t是时间。

步骤二：将产品退化的物理模型转换成线性模型的形式。

根据加速度计的物理模型可知，其含有温度T的项与含有时间t的项是乘积关系，因而可以分别将两项做代换得物理模型的线性模型形式。即做代换C_t＝ln(t+1)，可得y＝aB_SC_t。

步骤三：对物理模型的线性形式进行参数估计。

本发明中，采用了最小二乘法对变换成线性形式的物理模型进行参数估计。对于第j个样本的第i个性能退化点，模型与实验数据的偏差可得所有样本的偏差的平方和：

${\mathrm{SS}}_{\mathrm{Ej}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{\left[e\left(t_{\mathrm{ji}}\right)\right]}^2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{[\mathrm{a\; exp}\left(\frac{b}{T_j}\right)t_{\mathrm{ji}}-K_{1e}^{t_{\mathrm{ji}}}]}^2$

将图2所示3个加速度计的历史数据代入式中，利用奈尔德-米德单纯性法求解偏差的平方和的最小值。可以求得当a为39.09，b为-18.83时偏差的平方和取得最小值。

步骤四：根据产品性能参数退化轨迹的特点，选择合适的随机过程，构造物理统计模型。

加速度计的标度因数退化量的增量是相互独立，考虑使用漂移布朗运动过程构造统计模型。在此做代换μ＝aB_S，并为物理模型的线性形式附加上布朗运动的扩散项，可以加速度计的物理统计模型为y＝μC_t+σB(C_t)，其中μ＝aB_S，C_t＝ln(t+1)。当确定温度应力T取为90℃时，μ取值为39.09。

步骤五：选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新。

显然，在加速度计标度因数的物理统计模型中，漂移系数μ决定的退化轨迹的趋势和退化速率的快慢，并且是与加速度计中某材料相关的参数。因此选择用卡尔曼滤波算法对漂移系数μ进行更新，以能够更准确的估计每一时刻的漂移系数μ_t。更新漂移参数的物理统计模型可以表述为：

$y_{t_{i+1}}=y_{t_i}+{\widehat{\mathrm{\μ}}}_{t_i}(C_{t_{i+1}}-C_{t_i})+\mathrm{\σB}(C_{t_{i+1}},C_{t_i})$

其中，是t_i时刻漂移参数的估计值，服从的正态分布。

步骤六：对含卡尔曼滤波更新的物理统计模型中的参数进行参数估计。

由加速度计的带更新漂移参数过程的物理统计模型可知，退化量增量服从正态分布，即有 $y_{t_{i+1}}-y_{t_i}~N({\widehat{\mathrm{\μ}}}_{t_i}(C_{t_{i+1}}-C_{t_i}),P_{t_{i+1}|t_i}{(C_{t_{i+1}}-C_{t_i})}^2+{\mathrm{\σ}}^2(C_{t_{i+1}}-C_{t_i})).$

因而可以求得极大似然函数为：

$L=\underset{j=1}{\overset{n}{\mathrm{\Π}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Π}}}\left(\frac{1}{\sqrt{2\mathrm{\π}(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}}\exp \right\{-\frac{{(y_{t_{j(i+1)}}-y_{t_{\mathrm{ji}}}-{\widehat{\mathrm{\μ}}}_{t_{\mathrm{ji}}}(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}^2}{2(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}\left\}\right)$

对极大似然函数两边求对数可得对数极大似然函数：

$\begin{array}{c}\log \left(L\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}(-\frac{1}{2}\ln \left(2\mathrm{\π}(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))\right))\\ -\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}\left(\frac{{(y_{t_{j(i+1)}}-y_{t_{\mathrm{ji}}}-{\widehat{\mathrm{\μ}}}_{t_{\mathrm{ji}}}(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}^2}{2(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}\right)\end{array}$

其中和分别有卡尔曼滤波更新过程可以计算得到：

${\widehat{\mathrm{\μ}}}_{t_{i+1}}={\widehat{\mathrm{\μ}}}_{t_i}+P_{t_{i+1}|t_i}(t_{i+1}-t_i){F_{t_{i+1}}}^{-1}({\mathrm{\Δy}}_{t_{i+1}}-{\mathrm{\Δy}}_{t_i}-{\widehat{\mathrm{\μ}}}_{t_i}(t_{i+1}-t_i))$

$P_{t_{i+1}|t_i}=P_{t_i|t_i}+Q$

$F_{t_{i+1}}={(t_{i+1}-t_i)}^2{P}_{t_{i+1}|t_i}+{\mathrm{\σ}}^2(t_{i+1}-t_i)$

$P_{t_{i+1}|t_{i+1}}=P_{t_{i+1}|t_i}-P_{t_{i+1}|t_i}{(t_{i+1}-t_i)}^2{F_{t_{i+1}}}^{-1}{P}_{t_{i+1}|t_i}$

需要进行寿命预测的四号样本温度应力T为90℃，由步骤四中的代换求得μ取值为39.09。以此为漂移参数更新的起点，将图3所示的试验数据代入对数极大似然函数，利用粒子群优化算法对多参数的对数极大似然函数求取最值，可得在参数P₀为64.75，Q为366，σ为15.61时，可以取得极大值。因而由此确定了模型的参数取值为39.09，P₀为64.75，Q为366，σ为15.61。

步骤七：求解产品的寿命分布和可靠度函数。

设加速退化试验截止于t_m，则可由卡尔曼滤波更新过程求得此时的漂移参数估计值并设加速度计标度因数的失效阈值y_L为280。以此为依据，可以计算加速度计标度因数寿命所服从的概率密度函数为：

$f\left(t\right)=\frac{(y_L-y_{t_m})}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}^2{(C_t-C_{t_m})}^3}}\exp \{-\frac{{((y_L-y_{t_m})-{\widehat{\mathrm{\μ}}}_{t_m}(C_t-C_{t_m}))}^2}{2{\mathrm{\σ}}^2(C_t-C_{t_m})}\}\frac{{\mathrm{dC}}_t}{\mathrm{dt}}$

从而可得图4四号加速度计寿命分布图。

利用可靠度函数与寿命分布概率密度函数的关系式可以求得：

$R\left(t\right)=\mathrm{\Φ}\left(\frac{(y_L-y_{t_m})-{\mathrm{\μ}}_{t_m}(C_t-C_{t_m})}{\sqrt{{\mathrm{\σ}}^2(C_t-C_{t_m})}}\right)-\exp \left(\frac{2(y_L-y_{t_m}){\mathrm{\μ}}_{t_m}}{{\mathrm{\σ}}^2}\right)\mathrm{\Φ}\left(\frac{-(y_L-y_{t_m})-{\widehat{\mathrm{\μ}}}_{t_m}(C_t-C_{t_m})}{\sqrt{{\mathrm{\σ}}^2(C_t-C_{t_m})}}\right)$

从而可得图5四号加速度计可靠度函数图。

本发明的研究分析和实例应用表面，对于符合本发明运用条件的产品性能退化数据，使用本发明中的物理统计模型能够很好的结合物理机理和试验数据的随机性，同时结合了卡尔曼滤波使得模型中的参数能够进行更新，更好的拟合了退化各阶段的特征。本发明切实有效的克服了传统的物理模型构建确定性退化模型的缺陷，传统的统计模型构建的退化模型无法明确参数的物理意义的缺陷，传统的退化模型忽视了退化各阶段退化特征可能不一致的缺陷，并能够便捷的得到产品的寿命分布和可靠度函数。

实际上，对于许多工程实际产品，建立精确的物理模型是有难度的，但是在物理统计模型中，由于随机过程模型的被融入其中，使得即使物理模型有所偏差的条件下，也可以通过随机项的弥补，更好的拟合退化轨迹，并且本发明的物理统计模型的参数具备了更新功能，对退化轨迹的拟合能力得到了进一步提升。

Claims (1)

1.一种基于物理统计模型的寿命预测方法，其特征在于：它包括以下具体步骤：

步骤一：构建产品退化的物理模型；

根据产品的使用信息，确定产品可能经受的环境应力和工作应力，依此对产品退化的物理模型进行调研，并结合对产品结构、工况方面的分析，确定可能的退化模型；对于单应力退化，通常其模型的形式是y＝f(t,S)，其中y是性能参数退化量，t是时间，S是加速应力；

步骤二：将产品退化的物理模型转换成线性模型的形式；

根据产品退化的物理模型，转换成线性模型形式的具体方法有以下三种：

1.若在物理模型中，含有应力水平S的项与含有时间t的项是乘积关系，则将两项分别做整体代换构建成线性模型，即若产品的加速退化物理模型形式为y＝Ag(S)h(t)，则分别作代换B_S＝g(S)，C_t＝h(t)，将物理模型代换成线性模型y＝AB_SC_t，其中B_S为与应力水平有关的系数，C_t是经过形式转换的时间量，在此模型中，性能参数的退化量与经过形式转换的时间量具有线性关系；

2.若在物理模型中，含有应力水平S的项与含有时间t的项是指数关系，则对模型等式两边先取对数，再分别作整体代换构建线性模型；即若产品的加速退化物理模型形式为y＝Ag(S)^h(t)或y＝Ah(t)^g(S)，先对等式两边取对数得ln y＝h(t)ln g(S)+lnA或ln y＝g(S)ln h(t)+ln A；对第一种情况作代换B_S＝lng(S)，C_t＝h(t)；对第二种情况作代换B_S＝g(S)，C_t＝ln h(t)，都能得到线性模型形式ln y＝B_SC_t+ln A；在此模型中，性能参数退化量的对数与经过形式转换的时间量具有线性关系；

3.若在物理模型中，含有应力水平S的项与含有时间t的项是其他关系，则根据具体情况用包括取对数、整体代换方法变换为线性模型的形式；

步骤三：对物理模型的线性形式进行参数估计；

采用最小二乘法，利用历史样本数据对变换成线性形式的物理模型进行参数估计，以线性模型y＝AB_SC_t说明如下：

1.整理产品性能参数的退化数据，将试验的真实时间量度按照线性形式中时间转换形式进行转换，从而将退化轨迹中退化量与真实时间的对应关系转换为退化量与经时间转换形式转换后的计算时间的对应关系；

2.计算物理模型在每一时刻求得的退化量与加速退化试验测得的退化量偏差，偏差其中表示在第j个样本在对应应力水平S_j下在第i个性能退化监测点时刻t_ji时测得的性能参数退化量；

3.计算物理模型在每一时刻求得的退化量与加速退化试验测得的退化量偏差的平方和，即对各应力水平下所有样本的偏差进行平方加和；偏差的平方和 $\mathrm{SS}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{[{e\left(C_{t_{\mathrm{ji}}}\right)]}^2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j}{\mathrm{\Σ}}}{[\mathrm{AB}}_{S_j}{C}_{t_{\mathrm{ji}}}-y_{t_{\mathrm{ji}}}]}^2;$

4.利用奈尔德-米德单纯性法求解偏差的平方和的最小值，当偏差的平方和取得最小值时的对应的参数值即为物理模型线性形式的参数估计值；

步骤四：根据产品性能参数退化轨迹的特点，选择合适的随机过程，构造物理统计模型；

在线性模型的基础上，依据产品性能退化的特点，选择合适的随机过程，构造得到的随机过程模型即为最终的物理统计模型；在选择适用于产品退化过程的随机过程时，首先需要对产品退化过程的特点进行分析：

1.产品性能参数的退化量是严格单调递增或递减的，选择伽马过程进行建模；

2.产品性能参数的退化量增量是相互独立，且符合正态分布，选择布朗运动进行建模；

3.产品性能参数的退化量增量是相互独立，但不符合正态分布，尝试选择逆高斯过程；

产品的性能参数退化的线性模型形式是y＝AB_SC_t，且退化量增量是相互独立且符合正态分布，则构造布朗运动随机过程模型y＝μC_t+σB(C_t)，其中μ为漂移布朗运动的漂移系数，μ＝AB_S，σ为漂移布朗运动的扩散系数，B(C_t)是以经过形式转换的时间量C_t为时间量的标准布朗运动，其服从N(0,C_t)正态分布；并且根据代换关系利用线性模型形式中已知的参数估计值计算得到物理统计模型中的参数估计；即可由μ＝AB_S计算得到μ的取值；

步骤五：选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新；

选择物理统计模型中与体现样本个体差异且与退化轨迹形态密切相关的参数作为重点参数，用卡尔曼滤波算法对其进行随退化时间的更新，使得物理统计模型能更好的拟合退化轨迹从而更精确的预测；

卡尔曼滤波算法能够利用目标的历史信息，滤去噪声的影响，得到一个较为准确的估计值，步骤四中的服从漂移布朗运动的物理统计模型y＝μC_t+σB(C_t)与样本间个体差异密切相关且与退化轨迹形态密切相关的是漂移系数μ，因此用卡尔曼滤波算法对漂移系数μ进行随时间的更新计算；

设初始时刻的服从正态分布N(μ₀,P₀)，设对于漂移系数从到存在的更新过程，其中η～N(0,Q)，则卡尔曼滤波从到的更新过程按下述式进行：

${\widehat{\mathrm{\μ}}}_{t_{i+1}}={\widehat{\mathrm{\μ}}}_{t_i}+P_{t_{i+1}|t_i}(t_{i+1}-t_i){F_{t_{i+1}}}^{-1}({\mathrm{\Δy}}_{t_{i+1}}-{\mathrm{\Δy}}_{t_i}-{\widehat{\mathrm{\μ}}}_{t_i}(t_{i+1}-t_i))$

$P_{t_{i+1}|t_i}=P_{t_i|t_i}+Q$

$F_{t_{i+1}}={(t_{i+1}-t_i)}^2{P}_{t_{i+1}|t_i}+{\mathrm{\σ}}^2(t_{i+1}-t_i)$

$P_{t_{i+1}|t_{i+1}}=P_{t_{i+1}|t_i}-P_{t_{i+1}|t_i}{(t_{i+1}-t_i)}^2{F_{t_{i+1}}}^{-1}{P}_{t_{i+1}|t_i}$

当漂移系统更新到时，对于t_i后的任意时刻t_k，预测退化轨迹的物理统计模型为： $y_{t_k}=y_{t_i}+{\widehat{\mathrm{\μ}}}_{t_i}(t_k-t_i)+\mathrm{\σB}(t_k,t_i);$

步骤六：对含卡尔曼滤波更新的物理统计模型中的参数进行参数估计；

利用产品性能参数的退化数据对含卡尔曼滤波更新物理统计模型中的参数进行参数估计，在此使用的参数估计方法是极大似然估计法，参数估计方法的流程是：

1.根据物理统计模型，确定一服从某一分布形式的统计量，漂移布朗运动模型中，退化量的增量即为服从正态分布的统计量；

2.根据物理统计模型，求出该统计量的表达式及其所服从的分布，步骤五的物理统计模型中，统计量的表达式即为两时刻点退化量的之差服从正态分布 $N({\mathrm{\μ}}_{t_i}(t_{i+1}-t_i),P_{t_{i+1}|t_i}{(t_{i-1}-t_i)}^2+{\mathrm{\σ}}^2(t_{i+1}-t_i));$

3.列出极大似然函数，服从均值为方差为的正态分布，则由连续分布的极大似然函数法得对数极大似然函数为：

$\begin{array}{c}\log \left(L\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}(-\frac{1}{2}\ln (2\mathrm{\π}(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))))\\ -\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m_j-1}{\mathrm{\Σ}}}\left(\frac{{(y_{t_{j(i+1)}}-y_{t_{\mathrm{ji}}}-{\widehat{\mathrm{\μ}}}_{t_{\mathrm{ji}}}(C_{t_{j(j+1)}}-C_{t_{\mathrm{ji}}}))}^2}{2(P_{t_{i+1}|t_i}{(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}})}^2+{\mathrm{\σ}}^2(C_{t_{j(i+1)}}-C_{t_{\mathrm{ji}}}))}\right)\end{array}$

其中j表示第j个样本，i表示共m_j个样品的性能参数退化轨迹点中的第i个性能参数退化轨迹点；

4.利用搜索算法对于对数极大似然函数求极值，并求出取得极值时对应的参数值；此处使用粒子群寻优算法，粒子群寻优算法是一种启发式智能搜索算法，其是通过沿搜索得到的解中最优解方向继续迭代搜索，能够快速有效的收敛于极值点；

步骤七：求解产品的寿命分布和可靠度函数；

依据物理统计模型能求解产品的寿命分布和可靠度函数，运用不同的统计模型构造的物理统计模型其寿命分布和可靠度函数的形式是不同的；利用漂移布朗运动构造的物理统计模型，产品的寿命分布，也即首穿时分布服从逆高斯分布，并由寿命分布求得可靠度分布，对与寿命分布为L＝f(t)的产品，其可靠度函数为 $R\left(t\right)=1-{\∫}_0^{t}f\left(s\right)\mathrm{ds}.$