CN110955963A - 一种航空电缆剩余寿命预测方法 - Google Patents

Abstract

一种航空电缆剩余寿命预测方法，从电缆失效概率所服从的威布尔分布出发，首先建立电缆实际寿命与表征寿命的特征参数间的威布尔模型，在利用极大似然法估计威布尔参数初值的基础上，采用贝叶斯估计法得到威布尔模型，从而减小在小样本失效数据情况下对参数估计的影响；其次建立表征寿命的特征参数与环境应力间的寿命预测模型，将失效物理模型与支持向量机相结合，通过计算权重分配模型比重得到最终的电缆寿命预测模型。本发明将贝叶斯估计中威布尔形状参数分为常数和随应力变化的两种情况，更加全面的考虑了电缆寿命在实际预测的情况；同时采用对数线性形式来描述失效物理模型，预测电缆在单应力或多应力耦合情况下的剩余寿命。

Description

一种航空电缆剩余寿命预测方法

技术领域

本发明属于电缆可靠性技术领域，具体涉及一种航空电缆剩余寿命预测方法。

背景技术

电缆将飞机内的电气系统和控制系统联系起来，是飞机内动力电源、控制信号以及数据信息传递的主要载体，健康的线缆是整个电气系统里最简单，却是最重要的部件，在飞机的飞行安全中有着相当重要的地位。近些年，一些在航空航天界发生的事件逐渐让人们意识到，线缆的绝缘老化是一个值得重视的问题，目前人们面对的最主要问题是缺乏一种方便可行的方法来评估航空线缆的剩余寿命。航空线缆老化绝缘状态在飞机服役期间，如果不能很好地被评估，有可能会造成经济上的损失甚至致命的灾难。例如有些电缆绝缘性能尚可满足正常工作的需要，在检修时就没有必要被更换掉；有的电缆在预计寿命内，但内部却发生了不可逆的损坏，有导致事故发生的危险。因此，对典型环境应力下的航空线缆的寿命进行准确的预计对保障飞行安全有着重要的意义，可以保证飞机的正常运行和飞行安全。

由于航空电缆可靠性很高，在正常工作状态获得寿命数据要耗费很大的时间、精力和成本，因此寿命数据的获取是建立在加速寿命试验的基础上。加速寿命试验(Accelerated life testing，ALT)，最早由美国罗姆航空发展中心于1967年给出统一定义：加速寿命试验是在合理的工程及统计假设的基础上，采用提高工作应力或环境应力的方法来使产品快速地暴露故障，再利用加速应力下的试验信息(失效数据、失效原因等)来外推产品在正常应力水平下的可靠性指标的一种寿命试验方法。简而言之，ALT是在保持失效机理不变的条件下，在高应力水平下对产品进行试验使其加速失效的一种寿命试验方法；其目的在于快速获取试验数据，迅速查明失效原因，用运加速加速模型对产品在正常应力水平下的各种可靠性指标进行统计推断。ALT有效地缩短了试验时间、降低了试验成本，使高可靠性长寿命产品在规定条件下的寿命与可靠性水平的快速评估成为可能。

不同于其他的故障类型，航空电缆绝缘老化是一个长期的、连续的甚至不能简单观测出来的复杂的状态，老化贯穿绝缘材料寿命的始终，很难找到精确表征电缆寿命的实际可检测特征参数。因此，本发明从概率分布的角度出发，在加速老化试验数据的基础上，通过参数估计法建立符合电缆失效分布规律的威布尔模型，得到表征电缆实际寿命的特征参数，从而建立典型应力影响下电缆单因素或多场耦合的多因素寿命模型，利用外推法预计电缆正常工作环境下的寿命。

威布尔分布是由最弱环节模型导出的，这个模型如同许多链环串联而成的一根链条，两端受拉力时，其中任意一个环断裂，则链条失效。对于遵循“最弱环”的产品，也就是由大量的、失效时间服从同一分布零件组成的产品，且该产品由各零件中最先发生失效的时间决定产品整体的失效时间，那么该产品的寿命就可以认为服从威布尔分布。威布尔分布对于产品寿命的“浴盆曲线”的三个失效期都有较强的适应力，能充分反映材料缺陷和应力集中源对材料疲劳寿命的影响，而且具有递增的失效率，所以将它作为航空电缆的寿命分布模型是十分合适的。在已知电缆寿命服从二参数威布尔分布后，需要根据总体分布和加速寿命数据求出具体参数。因此，研究威布尔分布模型的参数估计方法对航空电缆的可靠性分析是非常重要的。目前国内外主要基于两种理论研究二参数威布尔分布参数估计问题：经典统计理论和贝叶斯统计理论。

经典统计理论的观点是把获得的数据或样本视作来自特定概率分布的总体，该理论下的参数估计只利用了样本信息和总体信息。总体信息是指已确定的总体分布形式等信息；样本信息是指在总体中抽样所获取信息。但是来自于以前经验和历史资料的信息，作为重要信息来源却没有加以利用。在实际工作中，在抽样之前我们常常拥有有关统计问题的信息，也就是先验信息。贝叶斯统计学综合利用了以上三种信息来进行统计推断，对于小样本数据的参数估计具有很大的优势。在采用贝叶斯估计求解威布尔参数的过程中，游达章等把代表失效机理的形状参数视为常数，然而实际情况中，形状参数可能随应力变化而变化。因此，我们根据实际试验数据将情况分为两种，若每组代表失效机理的形状参数差别很小时，将形状参数视为常数；若每组试验形状参数差别较大时，将形状参数视为随应力变化的变量。分别在这两种情况下采用先验矩方法和bootstrap自助法求解先验分布，然后根据共轭先验分布的特性和MCMC方法得到后验分布并求出寿命特征参数。

在得到电缆寿命的威布尔分布后，需要建立表征电缆寿命的特征参数和环境应力之间的关系模型。目前有基于失效物理模型和数据驱动两种方法。基于失效物理模型的寿命预测是通过数学或物理模型来描述产品的退化过程，由于某些关键零件或材料发生了蒸发、扩散、氧化、腐蚀、移位、结晶等失效机理，对零件或材料中的基本物质产生分子、原子等微观级别的化学与物理变化，从而使产品的特性逐渐恶化，最终超出界限而失效。数据驱动方法使用事件数据和状态监视(CM)数据，并结合机器学习(Machine Learning，ML)技术，例如人工神经网络(Artificial Neural Network，ANNs)，支持向量机(Support VectorMachine，SVM)，关联向量机(Relevance Vector Machine，RVM)，来训练一个预测模型，然后利用训练后的模型来估计电缆的剩余寿命。由于这两种方法在建立寿命预测模型中各有利弊，越来越多的学者尝试将这两种方法结合在一起，封杨等采用将威布尔分布模型与SFAM神经网络融合的方法来预测轴承的剩余寿命；Linxia Liao等把粒子滤波法分别和两种数据驱动方法结合，第一种数据融合方法用于建立寿命模型以获得产品内在状态，第二种融合方法用于对产品的粒子更新及长阶段的寿命预测。这种数据融合的方法相较于单独采用基于失效模型的方法或数据驱动的方法具有更好的预测效果和优势，因此本发明提出将失效物理模型与数据驱动相结合的数据融合的方法，分别利用对数线性化的失效物理模型和支持向量机得到有关电缆寿命的预测模型，并通过最小化误差的权重法得到更为准确的寿命预测模型。

发明内容

本发明针对现有对航空电缆寿命预测方法的局限，提供小样本失效数据的贝叶斯参数估计和剩余寿命预测模型建立方法，使得电缆在加速寿命试验数据基础上，预测单一应力或多种应力耦合作用下的电缆剩余寿命。考虑了用贝叶斯估计威布尔参数中形状参数不变和变化的情况，分别用不同的方法去求解后验分布和寿命特征参数。此外，还采用了将失效物理模型和支持向量机的数据融合方法获得寿命预测模型，通过最小化寿命特征参数和模型预测结果的误差，得到分配权重，建立寿命预测模型，提高模型的适用性和准确性。

为实现上述目的，本发明采用以下技术方案：

一种航空电缆剩余寿命预测方法，其特征在于，包括如下步骤：

步骤一、利用极大似然估计法得到电缆寿命分布的似然函数和威布尔参数的初始值，威布尔参数包括代表失效机理的形状参数和代表实际剩余寿命的特征参数；

步骤二、通过小样本贝叶斯估计法求解威布尔特征参数，根据加速寿命试验数据，考虑形状参数随环境应力不变的情况，将威布尔分布转换为指数形式，采用先验矩方法求解先验参数，通过共轭先验分布得到特征参数的点估计结果；

步骤三、根据加速寿命试验数据，考虑形状参数随环境应力强度变化而变化的情况，利用bootstrap自助法求解含有两个先验参数的先验分布，通过贝叶斯估计原理得到后验分布的基础上，采用马尔科夫链蒙特卡罗中的Metropolis-Hastings采样算法求解后验分布，得到特征参数的值；

步骤四、基于得到的特征参数，通过电缆失效的机理，建立环境应力下电缆剩余寿命的失效物理模型；

步骤五、基于得到的特征参数，利用环境应力与特征参数关系，建立表征特征参数与环境应力间关系的支持向量机模型；

步骤六、根据失效物理模型和支持向量机模型，建立电缆特征寿命预测模型。

为优化上述技术方案，采取的具体措施还包括：

进一步地，所述步骤一的具体过程为：

通过NASA航空工程与安全中心(NESC)对电缆剩余寿命数据进行立项研究，统计得到电缆的剩余寿命t的分布数据符合二参数的威布尔分布，其概率分布函数和概率密度分布函数分别为：

其中，m为形状参数，表示加速寿命试验中产品的失效机理；η为特征参数，一般用于表征产品的实际剩余寿命；由于电缆可靠性很高，在短时间内施加正常工作环境应力很难获得电缆寿命数据，而加速老化试验是目前被广泛用于研究电缆老化失效的试验方法；设样品受到p(p≥1)种不同类型的典型环境应力s＝(s₁，s₂，…，s_p)，并将试验条件分为k(k≥2)组环境应力，在第i组定时截尾试验环境下，有n_i个电缆样品参与试验，t_iτ为试验截止的失效时间，r_i为试验样品的失效数，得到有关失效时间的试验数据为

其中i＝1，2，…，k；因此，完全寿命下的r_i个子样失效样品的似然函数为：

其中，下标i表示第i组定时截尾试验环境下；

n_i-r_i个未发生失效的子样似然函数为：

将完全寿命数据集和未发生失效数据集结合，定时截尾试验中关于m_i，η_i的似然函数为：

将似然函数取对数，得到：

通过极大似然函数，对m_i，η_i求偏导：

利用matlab求解上述超越方程，通过二分法对根所在的区间进行逐次二分，不断接近真实解的方法，得到关于威布尔参数m_i，η_i的初始值。

进一步地，所述步骤二的具体过程为：

贝叶斯理论是利用其验前分布和样本信息来计算其后验分布，从而估算该变量的点估计和置信区间，并进一步推导其他相关可靠性特征量的估计值；已知待估计的参数θ的先验分布和样本似然函数分别为π(θ)和p(x|θ)，其后验分布表达式为：

其中，样本似然函数p(x|θ)为L(m_i，η_i)；

在利用贝叶斯方法估计小样本的威布尔参数时，需要知道参数的先验分布，相关的先验信息包括相似产品加速应力试验信息、专家经验信息等；在权利要求2中提到，威布尔参数m，η分别代表产品的失效机理和实际剩余寿命；根据实际加速试验结果，情况分为两种，一种情况是根据极大似然法得到的各组的形状参数m_i基本不变，将形状参数视为常数，其值取各组结果加权平均值；另一种是将形状参数m_i看作随环境应力强度变化的变量，在此基础上求解相应的特征参数η_i；步骤二中主要考虑第一种情况；

在形状参数不变的条件下，取各组结果加权平均值为公共形状参数：

下述内容均在第i组试验中，为方便叙述和理解，表示试验组数的下标i已省略；

将威布尔分布转换为指数形式：

样品的似然函数为：

其中，f(x)是F(x)对应的概率密度函数，n为每组试验的样本数，r为每组截尾时间内失效样品数，

根据大量统计信息得到，θ的先验分布符合倒伽马形式，其具体分布表达式为：

其中，a，b为先验分布中的参数，Γ表示伽马函数，采用先验矩方法求取倒伽马分布中的超参数；通过先验信息计算出θ的均值E(θ)和方差V(θ)，根据先验矩方法E(θ)和V(θ)又同时满足：

利用matlab求解上述二元方程组，得到θ的先验分布π(θ)～IGa(a，b)；在贝叶斯统计中，由于伽马分布是指数分布的共轭先验分布，则后验分布与先验分布属于同类，均为伽马分布形式，即θ的后验分布形式同样为倒伽马分布IGa(a+r，b+Xs)：

其中选择平方损失函数，θ的均值即为其点估计：

由此得到形状参数为常数时每组试验表征寿命的特征参数。

进一步地，所述步骤三的具体过程为：

在形状参数随应力强度改变的情况下，形状参数和特征参数的估计都需要通过贝叶斯估计求解；根据电缆失效方面的大量统计信息，形状参数和特征参数的先验分布均为倒伽马形式，同样，表示试验组数的下标i已省略，具体如下：

利用参数化的bootstrap自助法求得先验分布，具体如下：首先，利用先验矩方法得到初始倒伽马分布π₀(η)～IGa(a，b)，π₀(m)～IGa(c，d)；通过参数化bootstrap法抽取N组再生样本X(η)＝(X₁，X₂，…，X_N)，Y(m)＝(Y₁，Y₂，…，Y_N)；N组再生样本超参数对的求解方法为最小二乘法，从而得到N组超参数对的值(a_l，b_l)，(c_l，d_l)，l＝1，2，…，N；最后，对N组超参数对求取平均值

得到最终的先验分布

其具体步骤如图2；

在得到关于威布尔参数m，η的先验分布后，关于电缆剩余寿命的后验分布为：

π(m，η|x)∝p(x|m，η)π(m)π(η)

由于此贝叶斯估计中含有两个参数，在计算后验分布时涉及相关参数的双重积分，很难用数值方法计算出来，因此采用马尔可夫链蒙特卡罗(MCMC)方法实现，其基本原理是利用一个马尔可夫链来对目标分布进行采样，由于采样的目标分布无法直接进行采样，为马尔可夫链设计一个状态转移算子：θ⁽¹⁾→θ⁽²⁾→…θ^(t)→，然后利用一个建议分布生成一个新的候选状态，根据一定的概率选择接受这个新值或者拒绝这个新值，最终使马尔可夫链的稳态分布与目标分布吻合；采用马尔可夫链蒙特卡罗方法中的Metropolis-Hastings算法对后验分布函数进行抽样，对于含有m，η两个参数的后验分布的求解如下：

设置马尔科夫链初始次数t_s＝1，初始参数值μ＝(m，η)，状态转移算子

循环下列步骤直至t_s＝T，T表示迭代次数：

1)从建议分布

中生成θ^*；

2)计算接受率

3)从(0，1)均匀分布中分别生成u；

4)若u≤α，则接受建议分布，并令

否则令

由此得到形状参数随环境应力改变时每组试验的特征参数。

进一步地，所述步骤四的具体过程为：

基于失效物理模型的寿命预测是通过数学或物理模型来描述产品的退化过程，由于某些关键零件或材料发生了蒸发、扩散、氧化、腐蚀、移位、结晶等失效机理，对零件或材料中的基本物质产生分子、原子等微观级别的化学与物理变化，从而使物品的特性逐渐恶化，最终超出界限而失效。对于如热应力、电应力、振动应力等典型单环境应力下，电缆的寿命特征参数服从Arrhenins方程、逆幂律模型、指数分布模型；在热-电、热-振动、热-电-振动等多种应力耦合作用的情况下服从广义Eyring模型、FALLOU模型、RAMU模型等；为了统一描述在单应力或多应力环境下的电缆失效模型，在p种典型环境应力作用下，用对数线性化形式描述其失效物理模型：

其中，η为威布尔参数中表征电缆剩余寿命的特征参数，X_j为与p种典型环境应力s＝(s₁，s₂，…，s_p)相关的表达式，h_j为系数，b_η为常数项；在已知加速寿命试验的每组应力大小和特征参数(s_i，η_i)情况下，利用最小二乘法拟合对数线性化模型，求出每组的参数h_ij，b_iη，从而得到电缆剩余寿命的失效物理模型为：

进一步地，所述步骤五的具体过程为：

用失效物理模型来表征环境应力与寿命特征参数的关系可能会存在一定的误差，因此从数据驱动角度出发，利用SVM方法从结构风险最小化的思想来求解对应关系，其基本原理如下：

假定n个维数为d的样本，则可以表示为(x_i，y_i)，x∈R^d，i＝1，2，…，n，存在一个超平面H：w·x+b＝0满足以下条件：

式中：w为权向量，ξ_i≥0为误差变量，C为优化惩罚参数，b为偏置，

为核函数，实现了从低维空间到高维空间的映射，从而将非线性问题转换为线性分类；根据上述条件建立拉格朗日函数：

其中，α_i为拉格朗日乘子；为了求拉格朗日函数的最小值，需要对w，b，ξ_i求偏导，并令等式为零，得出：

根据Mercer条件，内积

用核函数κ(x_i，x_j)来表示，因此，通过对偶性质，超平面约束条件可变为：

经过求解，可以得出最优回归函数为：

由此，通过SVM方法，可以建立表征寿命的特征参数与环境应力间的关系，从而得到电缆剩余寿命的支持向量机模型为：

其中，s代表p种典型环境应力s＝(s₁，s₂，…，s_p)，α_i为拉格朗日乘子，η_i为由贝叶斯估计得到的每组的特征寿命，

为核函数，b为偏置。

进一步地，所述步骤六中，在已经得到关于特征参数和环境应力的失效物理模型

和支持向量机模型

的基础上，引入权重的概念，用于分配这两个模型在最终电缆特征寿命预测模型中的比重，引入权重w后特征寿命的具体表达式为：

通过最小化特征参数和模型预测结果

间的误差

计算得到权重，从而得到电缆特征寿命预测模型

其具体流程如图3。

本发明的有益效果是：

1、不同于需要大量试验数据的传统参数估计方法，将更适用于航空电缆小样本寿命数据的贝叶斯估计用于威布尔分布模型的确定中，通过贝叶斯参数估计法得到电缆实际寿命与表征寿命特征参数的关系；

2、不同于将代表失效机理的形状参数视为常数的大部分文献，本发明在考虑加速寿命实际试验数据下，将形状参数分为不变和随应力变化的两种情况，在形状参数不变的情况下，将威布尔分布转换为指数分布，利用先验矩方法求出先验分布，再通过共轭先验分布得到后验分布；在形状参数随应力变化的情况下，利用bootstrap自助法求出先验分布，再通过MCMC方法中Metropolis-Hastings算法得到后验分布；

3、将单应力或者多应力环境应力下的失效物理模型统一为对数线性化形式，并引入权重分配失效物理模型和SVM模型在预测模型中的比重，通过数据融合的方法，使得模型的适用性更加广泛，并且提高了寿命预测的精确度。

附图说明

图1是航空电缆剩余寿命预测方法流程图。

图2是步骤三中求得先验分布的方法流程图。

图3是步骤六中求得电缆特征寿命预测模型的方法流程图。

图4a、4b分别是参数m，η的分布密度图。

图5a、5b分别是参数m，η迭代过程中的马尔科夫状态链。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

如图1所示的一种航空电缆剩余寿命预测方法，从电缆失效概率所服从的威布尔分布出发，首先建立电缆实际寿命与表征寿命的特征参数间的威布尔模型，在利用极大似然法估计威布尔参数初值的基础上，采用贝叶斯估计法得到威布尔模型，从而减小在小样本失效数据情况下对参数估计的影响，并且在求解威布尔分布中，将形状参数分为常数和随应力变化的两种情况，更加全面的考虑了电缆寿命在实际预测的情况；其次建立表征寿命的特征参数与环境应力间的寿命预测模型，采用对数线性形式来统一描述在单应力或多应力耦合作用下的失效物理模型，并将失效物理模型与支持向量机相结合，通过计算权重分配模型比重得到最终的电缆寿命预测模型。通过加速寿命试验数据建立表征电缆实际寿命的特征参数的预测模型，从而利用模型预测电缆在正常工作环境下的寿命特征参数。

一、极大似然法求似然函数和参数初值

本次加速寿命试验数据是对航空电缆进行定时截尾恒定应力试验获取的，在温度和振动应力两种典型环境应力(T，S)作用下，将应力强度分为五组，每组各有8个电缆样品，相关试验数据如表1所示：

表1

在测试周期(t_i-1，t_i)之间估计失效样品的失效时间时，采用插值法，即在(t_i-1，t_i)均匀插入r_i个点，则r_i个样品的失效时间为：

(其中j＝1，2，…，r_i)

在第i组定时截尾试验环境下，有n_i个电缆样品参与试验，t_iτ为试验截止的失效时间，r_i为试验样品的失效数，得到有关失效时间的试验数据为

其中i＝1，2，…，k。因此，完全寿命下的r_i个子样失效样品的似然函数为：

n_i-r_i个未发生失效的子样似然函数为：

将完全寿命数据集和未发生失效数据集结合，定时截尾试验中关于m_i，η_i的似然函数为：

将似然函数取对数，得到：

通过极大似然函数，对m_i，η_i求偏导：

利用matlab求解上述超越方程，通过二分法对根所在的区间进行逐次二分，不断接近真实解的方法，得到关于威布尔参数的m_i，η_i的初值。各组威布尔参数值如表2所示：

表2

组号	1	2	3	4	5
						形状参数m<sub>o</sub>	3.01	3.358	2.67	3.51	3.127
特征参数η<sub>o</sub>	2.29	24.84	1.54	15.88	1.406

以得到的威布尔参数初值为基础，步骤二、三分别将形状参数为常数和随应力变化的两种情况展开讨论，但如表2所示，每组的形状参数估计值差别较大，所以步骤四、五、六均以形状参数随应力变化的威布尔参数值为实际标准值。

二、形状参数为常数情况下的贝叶斯估计

在假定形状参数为常数的情况下，公共形状参数为：

在第i组试验中，将威布尔分布转换为指数形式：

样品的似然函数为：

其中，r为每组截尾时间内失效样品数，

根据大量统计信息得到，θ的先验分布符合倒伽马形式，其具体分布表达式为：

其中，a，b为先验分布中的参数，这里采用先验矩方法求取倒伽马分布中的超参数。通过先验信息计算出θ的均值E(θ)和方差V(θ)，根据先验矩方法E(θ)和V(θ)又同时满足：

在查询航空电缆相似产品加速应力试验信息后，利用matlab求解上述二元方程组，得到θ的先验分布π(θ)～IGa(a，b)。

在贝叶斯统计中，由于伽马分布是指数分布的共轭先验分布，则后验分布与先验分布属于同类，均为伽马分布形式，即θ的后验分布形式同样为倒伽马分布IGa(a+r，b+Xs)：

其中选择平方损失函数，θ的均值即为其点估计：

在对每组试验数据进行形状参数为常数的贝叶斯估计后，得到的寿命特征参数的值如下表3所示：

表3

组号	1	2	3	4	5
						形状参数m	3.135	3.135	3.135	3.135	3.135
特征参数η	2.520	25.994	1.477	16.781	1.515

三、形状参数随应力变化情况下的贝叶斯估计

在形状参数随应力强度改变的情况下，形状参数和寿命特征参数的估计都需要通过贝叶斯估计求解。根据电缆失效方面的大量统计信息，形状参数和特征参数的先验分布均为倒伽马形式，同样，表示试验组数的下标i已省略：

由于加速老化试验数据较少，为了获得更为准确的先验分布，在利用先验矩方法得到先验参数a，b，c，d后，采用参数化的bootstrap自助法，该方法是根据已获得的先验数据求得先验分布，并从中作有放回的再抽样，将小样本的参数估计转化为大样本的参数估计问题，利用多次模拟抽样的结果逼近实际值，更适用于小样本的数据统计推断。其具体步骤如图2。

在得到关于威布尔参数m，η的先验分布后，关于电缆剩余寿命的后验分布为：

π(m，η|x)∝p(x|m，η)π(m)π(η)

由于此贝叶斯估计中含有两个参数，在计算后验分布时涉及相关参数的双重积分，很难用数值方法计算出来，因此采用马尔可夫链蒙特卡罗(MCMC)方法实现，其基本原理是利用一个马尔可夫链来对目标分布进行采样，由于采样的目标分布无法直接进行采样，为马尔可夫链设计一个状态转移算子：θ⁽¹⁾→θ⁽²⁾→…θ^(t)→，然后利用一个建议分布生成一个新的候选状态，根据一定的概率选择接受这个新值或者拒绝这个新值，最终使马尔可夫链的稳态分布与目标分布吻合。此文我们采用MCMC方法中的Metropolis-Hastings算法对目标函数(后验分布函数)进行抽样。对于含有m，η两个参数的后验分布的求解如下：

设置马尔科夫链初始次数t_s＝1，初始参数值μ＝(m，η)，状态转移算子

循环下列步骤直至t_s＝T，T表示迭代次数：

1)从建议分布

中生成θ^*；

2)计算接受率

3)从(0，1)均匀分布中分别生成u；

4)若u≤α，则接受建议分布，并令

否则令

以第一组试验为例，由于先验信息数据较少，利用matlab中的wblrnd(2.29，3.01)生成先验信息，通过先验矩方法得到关于m，η的先验分布参数值(a，b)，(c，d)，在得到先验分布的基础上，利用bootstrap自助法重复N＝5000次，通过matlab自带gamfit函数求解N组参数值，最后取平均值作为最终的先验分布参数。接下来需要根据贝叶斯公式求解后验分布和特征参数的值，通过Metropolis-Hastings算法编写matlab程序，本例中选取迭代次数T＝10000，得到后验分布的概率分布图、威布尔参数的马尔科夫状态链和估计值，为了更好的了解迭代过程中的参数变化，采用OpenBUGS软件辅助得出威布尔参数的分布密度图和马尔科夫状态链。图4a、4b分别给出了m，η参数的分布密度图，图5a、5b分别是m，η迭代过程中的马尔科夫状态链。其它四组参数估计过程与第一组相同。最后给出在形状参数随应力变化的情况下五组试验的威布尔参数值如表4所示：

表4

组号	1	2	3	4	5
						形状参数m	2.94	3.287	2.432	3.259	2.869
特征参数η	2.489	25.068	1.649	16.037	1.35

四、建立失效物理模型

基于失效物理模型的寿命预测是通过数学或物理模型来描述产品的退化过程，由于某些关键零件或材料发生了蒸发、扩散、氧化、腐蚀、移位、结晶等失效机理，对零件或材料中的基本物质产生分子、原子等微观级别的化学与物理变化，从而使物品的特性逐渐恶化，最终超出界限而失效。对于如热应力、电应力、振动应力等典型单环境应力下，电缆的寿命特征参数服从Arrhenins方程、逆幂律模型、指数分布模型；在热-电、热-振动、热-电-振动等多种应力耦合作用的情况下服从广义Eyring模型、FALLOU模型、RAMU模型等；为了统一描述在单应力或多应力环境下的电缆失效模型，在p种典型环境应力作用下，用对数线性化形式描述其失效物理模型：

其中，η为威布尔参数中表征电缆剩余寿命的特征参数，X_j为与某种应力相关的表达式，h_j为系数，b_η为常数项。在已知加速寿命试验的每组应力大小和寿命特征参数(s_i，η_i)情况下，利用最小二乘法拟合对数线性化模型，求出每组的参数h_ij，b_iη，从而得到电缆剩余寿命的失效物理模型。

本例中，我们选取的是在热应力和振动应力耦合作用下的情况，一般服从广义Eyring方程，其对数线性化形式可写为：

利用最小二乘法可以得到相关的系数：h₁＝1751.9，h₂＝-1.3，b_η＝-3.96

所以热-振动应力下电缆的失效物理模型表达式为：

五、建立SVM模型

用失效物理模型来表征环境应力与寿命特征参数的关系可能会存在一定的误差，因此我们从数据驱动角度出发，利用SVM方法从结构风险最小化的思想来求解对应关系。其基本原理如下：

假定n个维数为d的样本，则可以表示为(x_i，y_i)，x∈R^d，i＝1，2，…，n，存在一个超平面H：w·x+b＝0满足以下条件：

式中：w为权向量，ξ_i≥0为误差变量，C为优化惩罚参数，b为偏置，

为核函数，实现了从低维空间到高维空间的映射，从而将非线性问题转换为线性分类。根据上述条件建立拉格朗日函数：

其中，α_i为拉格朗日乘子。为了求拉格朗日函数的最小值，需要对w，b，ξ_i，求偏导，并令等式为零，得出：

根据Mercer条件，内积

用核函数κ(x_i，x_j)来表示，因此，通过对偶性质，超平面约束条件可变为：

经过求解，可以得出最优回归函数为：

详细的建立SVM预测的步骤如下：

1)对训练数据归一化处理到特定区间[y_min，y_max]，归一化公式为：

2)选择并应用核函数，本例选用高斯核函数：

κ(x，x_i)＝exp(-γ||x-x_i||²)

其中，γ为核函数参数；

3)通过交叉验证和网格搜索法获得惩罚因子和核函数参数；

4)根据获得的参数训练数据并预测结果。

以每组热-振动应力强度为自变量，寿命特征参数为因变量，利用上述步骤，得到电缆寿命预测方程：

主要的SVM参数选取如下表5所示：

表5

相关参数	选值
		核函数	高斯核函数
参数选择算法	交叉验证和网格搜索法
		惩罚因子	0.5
核函数参数	16
		核函数	1.291

六、权重法建立寿命预测模型

在步骤四和步骤五的基础上，已经得到关于寿命特征参数和环境应力的失效寿命模型和SVM模型，引入权重的概念，用于分配这两个模型在最终预测模型中的比重，通过最小化寿命特征参数和模型预测结果间的误差得到权重，其具体流程如图3。

通过权重法得到权重的值w＝0.836，计算每一种方法各组试验的误差平方与实际值比值的总和

得到最终的预测结果与失效物理模型、SVM模型结果如表6所示：

表6

	失效物理模型	SVM	权重法	实际特征参数
					1	2.344	2.73	2.4076	2.489
2	24.668	19.536	23.827	25.068
					3	1.540	1.87	1.5942	1.649
4	14.520	15.796	14.730	16.037
					5	2.153	1.592	1.961	1.35
计算总误差	0.643	1.321	0.547

从误差结果来看，利用权重法将失效物理模型与SVM方法结合在一起，对于特征参数的预测准确度更好，并且对于各种加速寿命试验下的电缆寿命预测模型具有更好的通用性和适用性。

需要注意的是，发明中所引用的如“上”、“下”、“左”、“右”、“前”、“后”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (7)

1.一种航空电缆剩余寿命预测方法，其特征在于，包括如下步骤：

步骤一、利用极大似然估计法得到电缆寿命分布的似然函数和威布尔参数的初始值，威布尔参数包括代表失效机理的形状参数和代表实际剩余寿命的特征参数；

步骤二、通过小样本贝叶斯估计法求解威布尔特征参数，根据加速寿命试验数据，考虑形状参数随环境应力不变的情况，将威布尔分布转换为指数形式，采用先验矩方法求解先验参数，通过共轭先验分布得到特征参数的点估计结果；

步骤三、根据加速寿命试验数据，考虑形状参数随环境应力强度变化而变化的情况，利用bootstrap自助法求解含有两个先验参数的先验分布，通过贝叶斯估计原理得到后验分布的基础上，采用马尔科夫链蒙特卡罗中的Metropolis-Hastings采样算法求解后验分布，得到特征参数的值；

步骤四、基于得到的特征参数，通过电缆失效的机理，建立环境应力下电缆剩余寿命的失效物理模型；

步骤五、基于得到的特征参数，利用环境应力与特征参数关系，建立表征特征参数与环境应力间关系的支持向量机模型；

步骤六、根据失效物理模型和支持向量机模型，建立电缆特征寿命预测模型。

2.如权利要求1所述的一种航空电缆剩余寿命预测方法，其特征在于：所述步骤一的具体过程为：

电缆的剩余寿命t的分布数据符合二参数的威布尔分布，其概率分布函数和概率密度分布函数分别为：

其中，m为形状参数，表示加速寿命试验中产品的失效机理；η为特征参数，用于表征产品的实际剩余寿命；

采用加速寿命试验作为用于研究电缆老化失效的试验方法，设样品受到p(p≥1)种不同类型的典型环境应力s＝(s₁，s₂，…，s_p)，并将试验条件分为k(k≥2)组环境应力，在第i组定时截尾试验环境下，有n_i个电缆样品参与试验，t_iτ为试验截止的失效时间，r_i为试验样品的失效数，得到有关失效时间的试验数据为

其中i＝1，2，…，k；因此，完全寿命下的r_i个子样失效样品的似然函数为：

其中，下标i表示第i组定时截尾试验环境下；

n_i-r_i个未发生失效的子样似然函数为：

将完全寿命数据集和未发生失效数据集结合，定时截尾试验中关于m_i，η_i的似然函数为：

将似然函数取对数，得到：

通过极大似然函数，对m_i，η_i求偏导：

利用matlab求解上述超越方程，通过二分法对根所在的区间进行逐次二分，不断接近真实解的方法，得到关于威布尔参数m_i，η_i的初始值。

3.如权利要求2所述的一种航空电缆剩余寿命预测方法，其特征在于，所述步骤二的具体过程为：

贝叶斯理论是利用其验前分布和样本信息来计算其后验分布，从而估算该变量的点估计和置信区间，并进一步推导其他相关可靠性特征量的估计值；已知待估计的参数θ的先验分布和样本似然函数分别为π(θ)和p(x|θ)，其后验分布表达式为：

其中，样本似然函数p(x|θ)为L(m_i，η_i)；

根据实际加速试验结果，情况分为两种，一种情况是根据极大似然法得到的各组的形状参数m_i基本不变，将形状参数视为常数，其值取各组结果加权平均值；另一种是将形状参数m_i看作随环境应力强度变化的变量，在此基础上求解相应的特征参数η_i；

在形状参数不变的条件下，取各组结果加权平均值为公共形状参数：

第i组试验中，将威布尔分布转换为指数形式：

样品的似然函数为：

其中，f(x)是F(x)对应的概率密度函数，n为每组试验的样本数，r为每组截尾时间内失效样品数，

θ的先验分布符合倒伽马形式，其具体分布表达式为：

其中，a，b为先验分布中的参数，Γ表示伽马函数，采用先验矩方法求取倒伽马分布中的超参数；通过先验信息计算出θ的均值E(θ)和方差V(θ)，根据先验矩方法E(θ)和V(θ)又同时满足：

利用matlab求解上述二元方程组，得到θ的先验分布π(θ)～IGa(a，b)；在贝叶斯统计中，由于伽马分布是指数分布的共轭先验分布，则后验分布与先验分布属于同类，均为伽马分布形式，即θ的后验分布形式同样为倒伽马分布IGa(a+r，b+Xs)：

其中选择平方损失函数，θ的均值即为其点估计：

由此得到形状参数为常数时每组试验表征寿命的特征参数。

4.如权利要求3所述的一种航空电缆剩余寿命预测方法，其特征在于，所述步骤三的具体过程为：

在形状参数随应力强度改变的情况下，形状参数和特征参数的估计都需要通过贝叶斯估计求解；形状参数和特征参数的先验分布均为倒伽马形式，具体如下：

利用参数化的bootstrap自助法求得先验分布，具体如下：首先，利用先验矩方法得到初始倒伽马分布π₀(η)～IGa(a，b)，π₀(m)～IGa(c，d)；通过参数化bootstrap法抽取N组再生样本X(η)＝(X₁，X₂，…，X_N)，Y(m)＝(Y₁，Y₂，…，Y_N)；N组再生样本超参数对的求解方法为最小二乘法，从而得到N组超参数对的值(a_l，b_l)，(c_l，d_l)，l＝1，2，…，N；最后，对N组超参数对求取平均值

得到最终的先验分布

在得到关于威布尔参数m，η的先验分布后，关于电缆剩余寿命的后验分布为：

π(m，η|x)∝p(x|m，η)π(m)π(η)

采用马尔可夫链蒙特卡罗方法中的Metropolis-Hastings算法对后验分布函数进行抽样，对于含有m，η两个参数的后验分布的求解如下：

设置马尔科夫链初始次数t_s＝1，初始参数值μ＝(m，η)，状态转移算子

循环下列步骤直至t_s＝T，T表示迭代次数：

1)从建议分布

中生成θ^*；

2)计算接受率

3)从(0，1)均匀分布中分别生成u；

4)若u≤α，则接受建议分布，并令

否则令

由此得到形状参数随环境应力改变时每组试验的特征参数。

5.如权利要求2所述的一种航空电缆剩余寿命预测方法，其特征在于，所述步骤四的具体过程为：

为了统一描述在单应力或多应力环境下的电缆失效模型，在p(p≥1)种典型环境应力作用下，用对数线性化形式描述其失效物理模型：

其中，η为威布尔参数中表征电缆剩余寿命的特征参数，X_j为与p种典型环境应力s＝(s₁，s₂，…，s_p)相关的表达式，h_j为系数，b_η为常数项；在已知加速寿命试验的每组应力大小和特征参数(s_i，η_i)情况下，利用最小二乘法拟合对数线性化模型，求出每组的参数h_ij，b_iη，从而得到电缆剩余寿命的失效物理模型为：

6.如权利要求2所述的一种航空电缆剩余寿命预测方法，其特征在于，所述步骤五中，从数据驱动角度出发，利用SVM方法从结构风险最小化的思想来求解表征寿命的特征参数与环境应力间的关系，从而得到电缆剩余寿命的支持向量机模型为：

其中，s代表p种典型环境应力s＝(s₁，s₂，…，s_p)，α_i为拉格朗日乘子，η_i为由贝叶斯估计得到的每组的特征寿命，

为核函数，b为偏置。

7.如权利要求1所述的一种航空电缆剩余寿命预测方法，其特征在于，所述步骤六中，在已经得到关于特征参数和环境应力的失效物理模型

和支持向量机模型

的基础上，引入权重的概念，用于分配这两个模型在最终电缆特征寿命预测模型中的比重，引入权重w后特征寿命的具体表达式为：

通过最小化特征参数和模型预测结果

间的误差

计算得到权重，从而得到电缆特征寿命预测模型

Images