CN102542155A - 基于加速退化数据的粒子滤波剩余寿命预测方法 - Google Patents

Abstract

本发明提出一种基于加速退化数据的粒子滤波剩余寿命预测方法，适用于机电产品、电子产品的寿命评估。该方法首先获取内场加速退化试验数据的相关信息，并建立性能退化方程，通过性能退化方程建立产品的状态方程和观测方程，将内场加速退化试验数据和少量外场数据作为先验信息，进行粒子滤波初始粒子抽样，并随时间对粒子权值进行更新和归一化粒子重采样，确定当前时刻产品的状态的最小均方估计，根据状态方程和观测方程预测下一刻时刻产品的状态，最后对产品的剩余寿命进行预测并确定预测结果的置信度。本发明方法避免单纯利用内场试验数据而忽略外场使用数据预测所造成的偏差，极大的减小了外场数据的收集时间，更具有工程使用价值。

Description

基于加速退化数据的粒子滤波剩余寿命预测方法

技术领域

本发明涉及一种产品的剩余寿命预测技术，适用于机电产品、电子产品的寿命评估，具体是一种基于加速退化数据的粒子滤波剩余寿命预测方法。

背景技术

近年来，故障预测与健康管理(PHM)无论在学术研究还是工业应用都变得越来越广泛。PHM方法允许在系统的实际的寿命周期条件下对其可靠性进行评估和进行剩余寿命预测，以预测故障将要发生的时间与地点，从而消除系统的风险。PHM过程中对在线数据的监测和学习是非常重要的一个环节，实现这一环节的方法有很多。粒子滤波由于其具有良好的逻辑推理性而成为了其中非常重要的一种方法。

粒子滤波是一种基于序贯蒙特卡罗方法和递推贝叶斯估计的统计滤波方法。它依据大数定理采用序贯蒙特卡罗方法来求解贝叶斯估计中的积分问题，其核心内容包含三方面：动态空间模型、贝叶斯估计理论和序贯蒙特卡罗抽样。

动态空间模型是指一个随时间相关随机系统，可以用两个模型来描述，一个为状态转移方程x_k＝f(x_k-1，v_k-1)，一个为观测方程z_k＝h(x_k，n_k)。其中x_k是代表系统在时间k的状态变量的理论值，x_k-1是代表系统在时间k-1的状态变量的理论值；z_k为系统在时间k的观测值，v_k-1为k-1时刻过程噪声；n_k为k时刻观测噪声。

贝叶斯估计理论是一种基于贝叶斯条件概率的估计理论。贝叶斯估计理论将未知量x看作是一个随机变量，并引入它的先验分布p(x)。在没有观测数据可以利用的情况下，只能根据以前的经验对x作出判断，进而估计x的分布。这样建立起来的分布即为关于x的先验分布p(x)；但是如果获得了观测数据z，z为后验信息，则可以根据以下公式对未知参数x的分布进行修订，p(x|z)为基于后验信息z的x的概率分布，简称为x的后验分布，p(z)为z的概率分布，p(z|x)为基于先验信息x的z的概率分布。

$p\left(x\right|z)=\frac{p\left(z\right|x)p\left(x\right)}{p\left(z\right)}=\frac{p\left(z\right|x)p\left(x\right)}{\∫p\left(z\right|x)p\left(x\right)\mathrm{dx}}$

序贯蒙特卡罗抽样将统计学中的序贯抽样与蒙特卡洛抽样方法相结合，从而实现后验概率密度的递推估计。序贯抽样方案是指在抽样时，不事先规定总的抽样个数，而是先抽少量样本，根据其结果，再决定停止抽样或继续抽样、抽多少，这样下去，直至决定停止抽样为止。蒙特卡罗方法是指当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。序贯蒙特卡罗抽样的基本思想是将实际求解的问题描述成某种随机变量，然后从已知的概率分布抽样，建立各种统计量，得到所求的解。

综上所述，粒子滤波的基本思想是：首先依据系统状态的经验条件分布在状态空间产生一组随机样本的集合，称这些样本为粒子，然后根据观测结果不断调节粒子的位置和权重，通过调整后的粒子信息修正最初的经验条件分布。其实质是用粒子及其权重组成的离散随机测度近似相关的概率分布，并且根据算法递推更新离散随机测度，当样本量足够大时，这种蒙特卡洛描述就近似于状态变量真实的后验概率密度函数。

粒子滤波在近年来发展迅速。Bhaskar，shane bulter，Marcos E.Orchard等人分别利用粒子滤波来对电池、半导体工艺中的减排设备涡轮发动机等进行剩余寿命预测，这些预测都取得了良好的结果。

然而要利用粒子滤波进行准确地剩余寿命预测和故障诊断需要大量的产品外场使用信息，收集这些信息往往需要耗费相当长得时间和很高的成本。另一方面，在产品投入外场使用后，其在设计研制过程中进行的诸多内场试验信息往往就被忽视了，实际上即使在有了外场使用数据的条件下，内场试验数据尤其是加速退化试验数据由于包含大量的产品寿命信息，仍然可以用来对产品进行剩余寿命预测和故障诊断。

加速退化试验(ADT-Accelerated Degradation Testing)通过对受试产品施加比起正常工作条件更为严酷的应力，来获得产品在加速应力下的性能退化数据，从而进行寿命预测和可靠性评估。通常试验不会进行到产品失效。

在加速退化试验中，应力有两种施加方式。一种是恒定的，即应力水平不随时间变化；一种是步进的，即应力水平随时间呈阶梯状递增或递减变化。相较于加速寿命试验，没有必要观测到产品失效的发生，从这一点来讲，在加速退化试验过程中更容易保持产品在加速过程中的失效机理不变。更重要的一点是加速退化试验不仅显著地降低了试验时间和成本，同时增加了试验数据的合理性。对加速退化试验的研究最早是在1979年，由Shiomi和Yanagisawa在进行薄膜电阻的加速寿命试验(ALT-Accelerated Life Testing)中开始进行的。此后Suzuki、Carey、Tseng、Hamada和Chiao、Padgett等人又分别对薄膜电阻、MOS、逻辑芯片、发光二极管、碳素纤维等产品进行了ADT的试验设计和数据评估探讨，并获得了一系列的研究成果。这些都是ADT较为成功的尝试。

目前，在产品还没有投入到外场使用时，要对其进行寿命预测，数据来源只能是在设计和研制过程中进行的各种内场试验，尤其是加速退化试验。在这个阶段，单纯地利用ADT数据来进行寿命预测是无可厚非，然而当产品投入到外场使用后，可以获得外场使用数据的时候，再单纯地利用ADT数据进行寿命预测，显然是不科学的。然而完全利用外场数据进行预测，无疑又造成了这些内场数据的损失，这些数据的价值就被忽略了，实际上这些内场数据在产品拥有一定的外场数据的时候仍然具有很高的利用和借鉴价值，尤其是外场数据不是很充分的条件下。

发明内容

本发明针对目前在机电产品、电子产品等进行外场使用中，对机电产品或电子产品的外场使用的剩余寿命进行评估时，若单纯地利用内场ADT数据进行寿命预测不科学，不利用则造成内场数据的损失，而完全利用外场数据则需要采集大量外场使用数据，而实际工程中往往不能够获得足够的外场使用数据，提出了一种基于加速退化数据的粒子滤波剩余寿命预测方法。

本发明提出的一种基于加速退化数据的粒子滤波剩余寿命预测方法，该方法所针对的产品，在加速退化试验中满足如下两个条件：

(a)产品的性能退化过程是随机过程；

(b)产品的性能退化过程具有一阶马尔科夫性。

具体包括如下步骤：

步骤一、获取内场加速退化试验数据的相关信息，包括：应力水平数、总的样本数、每个应力水平下的产品的性能监测值的采集次数和数据采集间隔，并进行预处理，具体是：如果同一个应力水平下的样本量大于1，则对处于同一个应力水平下的样本的对应时刻的加速退化试验数据求均值。

步骤二、利用加速退化试验数据建立性能退化方程：Y_i(k)＝y₀+d(S_i)g(k)+ε；其中，Y_i(k)为第i个应力水平下产品在k时刻所有样本的性能参数值的均值；d(S_i)为加速模型；g(k)为关于时间k的函数，采用线性、幂指数、对数和指数四种模型拟合；ε为随机误差项，y₀为产品所有的退化初值的均值。

步骤三、通过性能退化方程建立产品的状态转移方程x_k＝f(x_k-1，v_k-1)和观测方程z_k＝h(x_k，n_k)：其中，x_k是代表系统在k时刻产品的状态值，x_k-1是代表系统在k-1时刻产品的状态值，即理论上的性能值，z_k为系统在k时刻的观测值，是通过仪器等设备采集到的性能值，v_k-1为k-1时刻的过程噪声，n_k为k时刻的观测噪声，f：x_k-1→x_k，h：x_k→y_k为有界非线性映射；具体所述的状态转移方程通过将k时刻产品的状态与k-1时刻产品的状态做差得到。

步骤四、根据以上三步骤获取的结果作为先验信息，进行粒子滤波初始粒子抽样，具体是：产品在粒子滤波的0时刻，产品的状态服从高斯分布N(y₀，σ²Δk)，从高斯分布N(y₀，σ²Δk)中利用序贯蒙特卡罗抽样产生粒子群

为0时刻第i个粒子的状态值，N_s为抽样产生的粒子的总数，各粒子权值均为1/N_s。

步骤五、粒子权值的更新和归一化，具体是：首先，更新k时刻的粒子权值： $w_k^i=w_{k-1}^{i}p\left(z_k\right|\stackrel{\‾}{x_k^i})=w_{k-1}^{i}p(z_k-h\left(\stackrel{\‾}{x_k^i}\right)),i=\mathrm{1,2},...N_s,$ 其中，为k时刻第i个粒子的权值，

为k-1时刻第i个粒子的权值，

为k时刻第i个粒子的状态的估计值，

为基于

的z_k的概率分布，h(.)为状态值到观测值的映射，为k时刻第i个粒子对应的观测值，由观测方程可知

为噪声项。然后对更新后的粒子权值进行归一化处理：

步骤六、粒子重采样，得到新的粒子集合

为k+1时刻第i个粒子。

步骤七、确定当前时刻产品的状态的最小均方估计。

步骤八、根据状态转移方程和观测方程预测下一刻时刻产品的状态。

步骤九、对产品的剩余寿命进行预测。

步骤十、确定产品状态能够接受的置信度。

本发明的优点和积极效果在于：

[1]本发明方法中使用的数据既包括内场加速退化数据也包括外场使用数据，避免单纯利用内场试验数据而忽略外场使用数据预测所造成的偏差。

[2]本发明方法也适用于那些难以建立性能退化模型的产品，虽然有些产品的性能退化过程不适合建立性能退化模型来进行预测，但是本发明方法中建立性能退化模型更主要的是为了使在粒子滤波初期通过抽样产生的粒子具有足够高的精度，接下来的更新和预测部分则主要利用的是粒子滤波算法，从而避免了由于不能建立模型而被迫使用非参数方法所带来的误差，工程应用范围更广。

[3]本发明方法中利用加速退化试验数据辅助建立状态转移方程和观测方程，使得两个方程的参数在一开始就具有较高的准确性，极大的减小了外场数据的收集时间，更具有工程使用价值。

[4]本发明方法通过状态转移方程可以获得产品在某时刻预测结果的可信度，在实际工程中，根据具体产品确定可接受的概率，在采用本发明方法得到的某时刻的产品的状态可接受的可信度，从而为工程人员提供可信度参考。

附图说明

图1为本发明的粒子滤波剩余寿命预测方法的整体步骤流程图；

图2为本发明实施例中对加速退化试验数据进行预处理后得到的性能数据分布图；

图3为本发明实施例采集的外场数据的分布图；

图4为本发明实施例中进行在线粒子滤波所得的结果示意图；

图5为本发明实施例进行离线剩余寿命预测的结果示意图；

图6为本发明实施例在线粒子滤波阶段的均方误差的示意图；

图7为本发明实施例离线阶段粒子的均方误差的示意图；

图8为本发明实施例选择的置信区间的示意图。

具体实施方式

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

结合背景技术中的技术，构想是否可以利用ADT数据把加速退化试验寿命评估方法与粒子滤波剩余寿命预测方法有机的结合起来以解决当前剩余寿命技术中存在的问题呢？如果可以实现上述构想，就可以解决由于外场使用数据少而不足以进行剩余寿命预测的问题，解决需要花费很长的时间去采集足够多的外场使用数据的问题，还解决了加速退化试验寿命评估方法而未考虑产品在正常工作应力下的性能表现数据的问题以及内场ADT数据得不到充分利用而造成信息量无形的损失的问题。

鉴于以上的分析和思考，本发明提出了把加速退化试验数据以及少量的外场使用数据作为先验信息，辅助建立产品的状态转移方程和观测方程，并把加速退化试验寿命评估方法与粒子滤波剩余寿命预测方法有机的结合起来，进而实现产品的剩余寿命预测。类似的研究目前还未见报道。而由于本发明的剩余寿命预测方法同时考虑了内场ADT数据和外场使用数据，通过增加对产品进行评估预测的数据量，提高了预测结果的准确度，所以其预测结果的精确度要好于单纯利用内场或者外场数据进行预测得到的结果。

在产品的寿命预测过程中，通常采用随机过程来描述产品的性能退化过程。目前，布朗运动及其推广已广泛的出现在许多纯科学领域中，如物理、经济、通信理论、生物、管理科学与数理统计等。同时，由于布朗运动与微分方程(如热传导方程等)有密切的联系，它已经成为概率与分析联系的重要渠道。为此，本发明有以下假设：

a产品的性能退化过程是随机过程；

b产品的性能退化过程具有一阶马尔科夫性。

本发明集合产品的内场步进加速退化试验数据和外场在线收集到的使用数据为例，介绍如何利用加速退化试验数据作为先验信息辅助建立状态转移方程和观测方程以及如何结合离线剩余寿命预测方法进行产品的基于加速退化试验数据的粒子滤波剩余寿命的预测方法。如图1所述，具体步骤流程如下：

步骤一、获取关于加速退化试验数据相关信息，假设应力水平数为m，总共有n个样本，数据采集间隔为Δt，应力水平从低到高依次为S₁＜S₂＜…S_m，S_i(i＝1，2，…，m)下的产品的性能监测值，即观测值z_k。在第i个应力下的采集次数为M_i，则各个应力下采集次数的总和为

如果同一个应力水平下样本量大于1，则首先对同一应力水平下的不同样本对应时刻的加速退化试验数据求均值以减小产品差异带来的影响。

步骤二、对加速退化试验数据进行回归分析处理，以此作为先验信息，建立性能退化方程并估计其参数，所述的性能退化方程为：Y_i(k)＝y₀+d(S_i)g(k)+ε；其中，Y_i(k)为第i个应力水平S_i下产品在k时刻的性能参数值，d(S_i)为加速模型；g(k)为关于时间k的函数，通常采用线性模型、幂指数模型、对数模型和指数模型四种模型拟合；y₀为退化初值，ε为随机误差项，假设ε为布朗运动，服从高斯分布N(0，σ²Δk)，σ为退化过程中的扩散系数。

利用线性模型、幂指数模型、对数模型和指数模型四种模型对性能退化数据进行拟合，并采用回归分析中的相关系数来决定时间函数g(k)的具体形式，相关系数越接近1则表明所选回归模型对性能退化数据的拟合程度越好。

步骤三、由于退化方程是基于时间的随机过程，而且前面又假设该产品的性能退化过程具有一阶马尔科夫性，所以通过把k时刻产品的状态值与k-1时刻产品的状态值做差就可以得到状态转移方程x_k＝f(x_k-1，v_k-1)。本发明方法中产品的状态值为理论上的性能值，观测值为通过仪器观测得到的性能值，所以本发明中认为产品的观测值是产品的状态值加上观测仪器的系统误差，从而可以确定观测方程z_k＝h(x_k，n_k)，其中，x_k是代表产品在时间k的状态值，x_k-1是代表产品在时间k-1的状态值；z_k为产品在时间k的观测值，v_k-1为k-1时刻的过程噪声，n_k为k时刻的观测噪声，f：x_k-1→x_k，h：x_k→y_k为有界非线性映射。

步骤四、再以步骤一、步骤二和步骤三骤获取的结果作为先验信息，进行粒子滤波初始化。本发明把内场加速退化试验数据以及少量的外场使用数据作为先验数据，步骤一至三中所得到各方程及其参数作为先验信息。

由于在粒子滤波的最初时刻，产品的观测值即为退化初值y₀，又在步骤二估计得到了ε服从高斯分布N(0，σ²Δk)，所以认为产品在最初时刻即粒子滤波的0时刻，产品的状态服从高斯分布N(y₀，σ²Δk)。从高斯分布N(y₀，σ²Δk)中利用序贯蒙特卡罗抽样产生粒子群

设定所有粒子权值均为1/N_s，N_s为抽样产生的粒子的总数，

为0时刻第i个粒子。

步骤五、粒子权值的更新和归一化。

在任意k时刻，更新粒子权值 $w_k^i=w_{k-1}^{i}p\left(z_k\right|x_k^i)=w_{k-1}^{i}p(z_k-h\left(x_k^i\right)),i=\mathrm{1,2},...N_s,$ 为k时刻第i个粒子的权值，

为k-1时刻第i个粒子的权值，为k时刻第i个粒子的状态估计值，

为基于

的z_k的概率分布，h(.)为状态值到观测值的映射，

为k时刻第i个粒子对应的观测值，由观测方程可知

为噪声项。因此在知道噪声项的分布类型的时候就可以利用该分布来进行粒子权值的更新。更新粒子权值后进行归一化

通过归一化处理，既保持了粒子的特性，又保持了抽样原理的准确。

步骤六、粒子重采样。

得到新的粒子集合

表示k+1时刻的第i个粒子。

步骤七、k时刻未知参数x的最小均方估计为

计算公式为

所采用的最小均方估计为众多的估计方法之一。

步骤八、预测。利用状态转移方程f预测未知参数

通过以上步骤，实现了对加速退化试验数据的预处理，建立了性能退化方程。根据退化方程和前面的假设1和2可以推导出状态转移方程x_k＝f(x_k-1，v_k-1)和观测方程z_k＝h(x_k，n_k)，还可以根据状态转移方程和观测方程预测下一时刻的产品的状态并且通过持续收集在线数据，不断地更新和修正状态转移方程和观测方程中的参数。

步骤九、剩余寿命预测。重复上面的步骤五至步骤八直至通过抽样产生的粒子具有足够的精确性进行剩余寿命预测。由于认为粒子已经具有了足够的精确性即通过抽样产生的粒子与对应时刻的状态值非常接近并且保持稳定，所以自进行剩余寿命预测的一刻起粒子的权值保持不变，如公式 $x_{t+k+1}=E\left[f(x_{t+k}^{\left(i\right)},v_{t+k})\right]=E\left[x_{t+k+1}^{\left(i\right)}\right],k=\mathrm{0,1,2}...,E\left[f(x_{t+k}^{\left(i\right)},v_{t+k})\right]$ 和都表示对t+k时刻所有粒子预测得到t+k+1时刻的产品的状态值的期望，为t+k时刻第i个粒子的产品状态的理论值；v_t+k为t+k时刻的过程噪声；x_t+k+1为t+k+1时刻产品的状态的预测结果。依据此公式进行递推就可以得到k时刻之后任意时刻的产品状态的预测结果。在得到这些结果后，能够绘制出如图5所示的图，其中纵坐标表示的是产品退化百分比，横坐标为产品退化到对应退化百分比的时间。从图中观察就可以得到产品性能退化到某一个百分比的时间。如果确定了一个产品退化到某一个百分比的时候产品失效，那么对应的时间就是产品的寿命。

步骤十、确定产品状态能够接受的置信度。

k+1时刻产品性能参数x_k+1服从正态分布N(u，R)，u＝x_k+d(T₀)Δt，R＝2σ²Δt，在工程实际中，小概率事件发生时，认为产品出现了异常，通常这个概率P₀可以根据具体产品或者系统而设定。设k+1时刻以P₀的概率接受当前产品的状态，即Pr(u-q≤x_k+1≤u+q)＝P₀，计算得到q，则[u-q，u+q]为k+1时刻产品状态的可以接受的区间。

本发明方法的基本思想就是：在进行剩余寿命预测的一刻起，即保持所有粒子的权值不变，任何一时刻的预测结果为当前时刻所有粒子的加权和。

实施例：本发明方法采用超辐射发光二级管的内场加速退化试验数据和外场使用数据来验证其正确性。

超辐射发光二极管步进应力加速退化试验应力为温度，样本量为3，试验监测参数为光功率，对所有样本的初始性能监测值求均值得到y₀，数据采集间隔为Δt₁＝0.5分钟，M1，M2，M3，M4分别为对应温度应力下的数据采集个数，3个样本都参与了4个温度应力水平的试验，对所有样本在同一时刻的性能监测值求均值为Y(k)。具体试验信息如下表1所示：

表1试验参数信息

T0	T1	T2	T3	T4
					25℃	60℃	80℃	100℃	110℃

y0	M1	M2	M3	M4
					121.6	4519	4293	1805	2052

1.建立退化方程和加速模型，采用漂移布朗运动来描述其退化过程。漂移布朗运动是在布朗运动的基础上添加了一个趋势项，这样既有布朗运动的随机性又具有趋势项的趋势性。根据内场产品的退化失效机理分析，考虑其加速试验应力为温度，选用Arrennius模型作为超辐射发光二极管的加速模型，即漂移系数抑或退化率d(S_i)与温度的关系为d(T_i)＝exp(A-Ea/(κ*T_i))，i＝1，2，…，m，m为应力水平数。Arrhenius模型是一个基于实验结果的经验公式，能够定量地给出化学反应速率与温度的关系。其中，将应力S写作绝对温度T；A为待估计参数，Ea为激活能，κ为波尔兹曼常数。根据步骤二可知，ε～N(0，σ²Δt)，则该产品的性能退化漂移布朗运动可表示为：

Y_i(k)＝ε+exp(A-Ea/(κ*T_i))·k+y₀，i＝1，2，…，m

2.采用退化百分比的方式来衡量产品的退化过程，得到三个样本的性能数据分布图如图2所示。退化百分比ζ为试验时某时刻的绝对性能退化量除以试验前性能初值，其计算公式为ζ＝[y₀-Y_i(k)]/y₀。由观察可以发现，产品的退化过程近似线性，所以对各个应力水平S_i下的试验数据以最小二乘法进行线性回归分析，得到的直线的斜率即为退化率d(T_i)。将Arrennius模型两边取对数，便得到了lg(d(T_i))＝A-Ea/(κ*T_i)，这是一个因变量为lg(d(T_i))关于自变量为1/(κ*T_i)的线性方程，由于lg(d(T_i))中d(T_i)已经求得，另外1/(κ*T_i)已知，所以再一次利用最小二乘法进行线性回归便可以得到A和Ea的估计值。A为回归直线的截距，Ea为斜率。将T₀＝25℃带入到Arrennius模型，可以求出T₀＝25℃下的退化率d(T₀)。然后利用步进应力加速退化试验关于求σ²极大似然估计的公式：

${\widehat{\mathrm{\σ}}}^2=\frac{1}{n\·(M-m)\mathrm{\Δt}}\underset{l=1}{\overset{n}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_i-1}{\mathrm{\Σ}}}{[q_{\mathrm{lij}}-\exp (\hat{A}-\widehat{\mathrm{Ea}}/(\mathrm{\κ}*T_i))]}^2$

其中，M表示各个应力下采集次数的总和；l＝1，2，…，n，表示第l个样本；i＝1，2，…，m，表示第i个应力；j＝1，2，…，m_i-1表示第j个性能监测值，M_i表示第i个应力下总的采集次数。q_lij＝y_lij-y_li(j-1)，其中y_lij为第l个样本在第i个应力水平下的第j个性能监测值，y_li(j-1)为第l个样本在第i个应力水平下的第j-1个性能监测值。由上述公式可以求得σ的估计值。详见表2。

表2参数估计结果

参数	σ	A	Ea	d(T0)
					估计值	5.5471e-08	9.0682	0.6267	2.1909e-08

3.利用试验数据估计得到的结果建立粒子滤波状态转移方程和观测方程。由于数据为计算机自动通过传感器采集记录，而传感器的测量精度非常高，在本实施例中甚至可以忽略，所以本实施例中认为观测方程不存在观测误差，即不需要建立观测方程。得到的状态转移方程为：

$x_{k+1}~N(x_k+d\left(T_0\right)\mathrm{\Δ}{t}_2,\mathrm{\Δ}{t}_1\mathrm{\σ}\sqrt{2})$

其中，x_k为超辐射发光二极管在k时刻的状态值，x_k+1为超辐射发光二极管在k+1时刻的状态值，d(T₀)为超辐射发光二极管在正常温度应力T₀＝25℃时的加速模型值，Δt₁为对超辐射发光二极管在内场加速退化试验时的数据采集间隔，Δt₂为对超辐射发光二极管在外场加速退化试验时的数据采集间隔，σ为超辐射发光二极管在加速退化试验过程中的扩散系数。

采集到产品在外场使用293320分钟的数据，数据采集间隔Δt₂＝10分钟，这些外场数据的分布如图3所示。选取该外场数据前1000分钟的数据，模拟在线更新，并进行在线粒子滤波，在线滤波结果如图4所示。

4.重复本发明的步骤五至步骤八直至到第1000分钟，通过抽样产生的粒子具有足够的精确性来进行剩余寿命预测，自第1000分钟起进行剩余寿命预测预测(Remaining UsefulLife，简称RUL)，得到结果如图5所示，图中：粗虚线为自第10000分钟剩余寿命预测阶段所得到的结果，细虚线为外场数据。

本发明方法通过以加速退化试验数据作为先验信息可以极大地缩短在线数据的收集时间，原因是在滤波的开始阶段，抽样产生的粒子已经具有足够高的精度，而如果不以加速退化试验数据为先验信息的话，则需要在线采集大量的数据进行更新修正，使抽样产生的粒子达到一定的精度。以下选取均方误差作为来衡量粒子滤波阶段(即前10000分钟)，抽样产生的粒子是否足够稳定，也就是抽样产生的粒子具有足够高的精确性来进行剩余寿命预测。均方误差的计算公式为

RMSE表示均方误差的大小，num是粒子的个数，

表示k时刻的第i个粒子的状态值，x_k表示在线滤波阶段是k时刻的在线采集到的产品的状态值，在剩余寿命阶段是由所有粒子加权得到的结果。计算得到在线粒子滤波阶段的均方误差如图6所示，离线阶段粒子的均方误差如图7所示。从图6和图7中可以发现，仅仅需要很少的时间就可以使抽样产生的粒子的精度达到很高的一个程度并且保持稳定。本发明方法确实可以做到通过将加速退化试验数据与外场使用数据相结合进行寿命预测，减少了对预测对外场数据的需求，缩短了进行寿命预测所需要的时间。

5.由状态转移方程可以得到k+1时刻产品性能参数x_k+1服从N(u，R)，u＝x_k+d(T₀)Δt，R＝2σ²Δt，在工程实际中，小概率事件生时，认为产品出现了异常，通常这个概率由使用方依据实际要求而设定。假设k+1时刻以P₀的概率接受当前状态，即P(u-q≤x_k+1≤u+q)＝P₀，根据下述公式计算得到q：

p(|x_k+1-u|≤q)＝p₀

$p\left(\right|\frac{x_{k+1}-u}{\sqrt{R}}|\≤\frac{q}{\sqrt{R}})=p_0$

显然

服从标准高斯分布N(0，1)，通过查询标准高斯分布表，可以得到进而计算得到

则[u-q，u+q]为该时刻产品状态可以接受的置信区间。

以滤波阶段第5000分钟时刻数据为例，取置信度为P₀＝0.8，x₅₀₀₀＝8.0814e(-4)，则x₅₀₁₀～N(x₅₀₁₀+d(T₀)Δt，R)。求解得到u＝8.0814e(-5)+2.1909e(-7)＝8.1033e(-5)，R＝1.2270e(-4)，q＝5.4873e(-4)。x₅₀₁₀置信区间[-4.6792e(-4)，6.2954e(-4)]，如图8所示。第5010分钟采集到的实际产品状态值为x₅₀₁₀＝7.1219e(-5)，落入了可接受区间，也就是说在这一步骤中，小概率事件并没有发生。产品第5010分钟的状态可接受，并未出现异常。

Claims (5)

1.一种基于加速退化数据的粒子滤波剩余寿命预测方法，其特征在于，假设产品的加速退化试验中，满足下面条件：

(a)产品的性能退化过程是随机过程；

(b)产品的性能退化过程具有一阶马尔科夫性；

该剩余寿命预测方法具体包括如下步骤：

步骤一、获取加速退化试验数据的相关信息，包括：应力水平数、总的样本数、每个应力水平下的产品的性能监测值的采集次数和数据采集间隔，并进行预处理，具体是：如果同一个应力水平下的样本量大于1，则对处于同一个应力水平下的样本的对应时刻的加速退化试验数据求均值；

步骤二、利用加速退化试验数据建立产品的性能退化方程：

Y_i(k)＝y₀+d(S_i)g(k)+ε                (1)

其中，Y_i(k)为第i个应力水平下产品在k时刻的所有样本的性能参数值的均值；d(S_i)为加速模型；g(k)为关于时间k的函数，采用线性、幂指数、对数和指数四种模型拟合；y₀为产品所有的退化初值的均值，ε为随机误差项，假设ε为布朗运动，服从分布N(0，σ²Δk)，σ为退化过程中的扩散系数，Δk为数据采集间隔；

步骤三、根据性能退化方程建立状态转移方程x_k＝f(x_k-1，v_k-1)和观测方程z_k＝h(x_k，n_k)；其中，x_k代表产品在k时刻的状态值；x_k-1代表产品在k-1时刻的状态值；z_k为产品在k时刻的观测值，是通过仪器观测得到的性能值；v_k-1为k-1时刻的过程噪声；n_k为k时刻的观测噪声；f：x_k-1→x_k，h：x_k→y_k为有界非线性映射；所述的状态转移方程通过将k时刻产品的状态值与k-1时刻产品的状态值做差得到；

步骤四、根据以上三步骤获取的结果作为先验信息，进行粒子滤波初始粒子抽样，具体是：产品在粒子滤波的0时刻，产品的状态服从高斯分布N(y₀，σ²Δk)，从高斯分布N(y₀，σ²Δk)中利用序贯蒙特卡罗抽样产生粒子群

为0时刻第i个粒子的状态值，N_s为抽样产生的粒子的总数，各粒子权值均为1/N_s；

步骤五、粒子权值的更新和归一，具体是：

首先，根据式(2)更新k时刻的粒子权值：

$w_k^i=w_{k-1}^{i}p\left(z_k\right|\stackrel{\‾}{x_k^i})=w_{k-1}^{i}p(z_k-h\left(\stackrel{\‾}{x_k^i}\right)),i=\mathrm{1,2},...N_s,---\left(2\right)$

其中，

为k时刻第i个粒子的权值，

为k-1时刻第i个粒子的权值，

为k时刻第i个粒子的状态的估计值，

为基于

的z_k的概率分布，为k时刻第i个粒子对应的观测值，

为噪声项；

然后，对更新后的粒子权值进行归一化处理：

步骤六、粒子重采样，得到新的粒子集合

为k+1时刻第i个粒子的状态值；

步骤七、采用式(3)确定k时刻产品的状态的最小均方估计

$\widehat{x_k}\≈\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}{w}_k^i{x}_k^i---\left(3\right)$

步骤八、根据状态转移方程和观测方程预测下一时刻产品的状态，并且通过持续收集在线数据，不断地更新和修正状态转移方程和观测方程中的参数；

步骤九、重复上面的步骤五至步骤八直至通过抽样产生的粒子具有足够的精确性，然后开始进行剩余寿命预测；所述的抽样产生的粒子具有足够的精确性，是指通过抽样产生的粒子与对应时刻的产品的状态值非常接近并且保持稳定；

步骤十、确定产品状态能够接受的置信度，具体是：k+1时刻产品的状态变量x_k+1服从正态分布N(u，R)，u＝x_k+d(S₀)Δk，R＝2σ²Δk，设k+1时刻以P₀的概率接受当前产品的状态，则有：

Pr(u-q≤x_k+1≤u+q)＝P₀                   (4)

根据式(4)得到q，则[u-q，u+q]为k+1时刻产品的状态能够接受的置信区间。

2.根据权利要求1所述的一种基于加速退化数据的粒子滤波剩余寿命预测方法，其特征在于，所述的产品为超辐射发光二极管，加速退化试验的应力为温度，选用Arrennius模型作为超辐射发光二极管的加速模型：d(T_i)＝exp(A-Ea/(κ*T_i))，i＝1，2，…，m，T_i为第i个应力水平下的温度，A为待估计参数，Ea为激活能，κ为波尔兹曼常数，则超辐射发光二极管的性能退化方程为：

Y_i(k)＝ε+exp(A-Ea/(κ*T_i))·k+y₀，i＝1，2，…，m                    (5)

其中，ε～N(0，σ²Δt)，Δt为数据采集间隔，y₀为温度为T₀时的退化初值，y₀通过对所有样本的初始性能监测值求均值得到。

3.根据权利要求2所述的一种基于加速退化数据的粒子滤波剩余寿命预测方法，其特征在于，所述的超辐射发光二极管的状态转移方程为：

$x_{k+1}~N(x_k+d\left(T_0\right)\mathrm{\Δ}{t}_2,\mathrm{\Δ}{t}_1\mathrm{\σ}\sqrt{2})---\left(6\right)$

其中，x_k为超辐射发光二极管在k时刻的状态值，x_k+1为超辐射发光二极管在k+1时刻的状态值，d(T₀)为超辐射发光二极管在正常温度应力T₀＝25℃时的加速模型值，Δt₁为对超辐射发光二极管在内场加速退化试验时的数据采集间隔，Δt₂为对超辐射发光二极管在外场加速退化试验时的数据采集间隔，σ为超辐射发光二极管在加速退化试验过程中的扩散系数。

4.根据权利要求2所述的一种基于加速退化数据的粒子滤波剩余寿命预测方法，其特征在于，所述的超辐射发光二极管在步骤九中采用均方误差RMSE来判断抽样产生的粒子是否具有足够的精确性，具体均方误差RMSE为：

$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{num}-1}\underset{i=1}{\overset{\mathrm{num}}{\mathrm{\Σ}}}{\left(\frac{x_k^i-x_k}{x_k}\right)}^2}---\left(7\right)$

其中，num是抽样粒子的总个数，

表示k时刻的第i个粒子的状态值，x_k表示在线滤波阶段是k时刻时在线采集到的超辐射发光二极管的状态值。

5.根据权利要求1所述的一种基于加速退化数据的粒子滤波剩余寿命预测方法，其特征在于，步骤九中所述的进行剩余寿命预测具体是：预测t+k+1时刻产品的状态

根据所得到的产品的状态预测值x_t+k+1判断该产品在t+k+1时刻是否出现故障，若出现故障，则对应的时间值就是该产品的寿命，根据所得到的寿命值能够得到产品当前的剩余寿命；其中，

和

都表示对t+k时刻所有粒子进行预测得到t+k+1时刻的产品的状态值的期望，

为t+k时刻第i个粒子的采样值，v_t+k为t+k时刻的过程噪声。

Images