CN114091790B - 一种融合现场数据和两阶段加速退化数据的寿命预测方法 - Google Patents

一种融合现场数据和两阶段加速退化数据的寿命预测方法 Download PDF

Info

Publication number
CN114091790B
CN114091790B CN202210068042.3A CN202210068042A CN114091790B CN 114091790 B CN114091790 B CN 114091790B CN 202210068042 A CN202210068042 A CN 202210068042A CN 114091790 B CN114091790 B CN 114091790B
Authority
CN
China
Prior art keywords
stage
degradation
accelerated
temperature stress
parameter
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
CN202210068042.3A
Other languages
English (en)
Other versions
CN114091790A (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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202210068042.3A priority Critical patent/CN114091790B/zh
Publication of CN114091790A publication Critical patent/CN114091790A/zh
Application granted granted Critical
Publication of CN114091790B publication Critical patent/CN114091790B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • 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
    • 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/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)
  • Business, Economics & Management (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Economics (AREA)
  • Strategic Management (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Development Economics (AREA)
  • Quality & Reliability (AREA)
  • Databases & Information Systems (AREA)
  • Game Theory and Decision Science (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Software Systems (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Control Of Electric Motors In General (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了一种融合现场数据和两阶段加速退化数据的寿命预测方法，根据不同热应力下电机绝缘材料试验样本的加速退化数据，利用SIC和残差平方和最小准则确定退化过程变点位置；构建热应力作用下基于两阶段随机过程的电机绝缘性能退化模型，根据Arrhenius方程得到模型参数和变点分布参数与温度间的关系，外推产品在正常温度下退化过程变点位置的概率密度函数及模型参数；利用待测电机绝缘材料的现场退化数据，用CUSUM控制图辨识变点，结合变点的概率密度函数，确定变点出现的概率；将加速退化数据作为模型参数的先验信息，再基于Bayes方法根据现场退化数据得到模型参数的后验估计值，得到首达时间意义下产品寿命的累积分布函数，实现电机绝缘的实时寿命预测。

Description

一种融合现场数据和两阶段加速退化数据的寿命预测方法
技术领域
本发明涉及设备寿命预测技术领域，尤其涉及一种融合现场数据和两阶段加速退化数据的寿命预测方法。
背景技术
为了应对能源短缺及大气污染问题，我国高度重视新能源发电和电动汽车的发展。风力发电机作为新能源发电的重要设备，因其长期运行于复杂、恶劣的环境工况下，发电机绝缘破坏是造成风力发电机故障的重要原因之一。电动汽车驱动电机因其频繁启动、变速、路况多变等运行特点，驱动电机故障主要是由电机绕组损坏引起的，而造成绕组损坏的重要原因之一是绝缘破坏。因此对电机绝缘进行寿命预测对提高产品的性能具有重大意义。
电机温度老化特性是影响绝缘寿命的关键因素，为保证电机绝缘的可靠性，往往在设计中提高安全裕度，例如增加铜导体截面积，以保证运行温度最低，或者增加主绝缘的厚度以避免电气击穿等。然而，提高安全裕度以牺牲材料的性能为代价，其最直接后果就是牺牲了电机的功率密度指标。通过降低不必要的安全裕度，可以降低制造成本，减小电机重量和体积，显著提升电机功率密度。降低绝缘材料安全裕度使电机的工作应力尽量接近其极限水平，对其可靠性提出了更高要求。如果能在电机运行过程中精准预测电机绝缘剩余寿命，及时做出合理、精确的检修、更换方案，可以大大提高电机的可靠性。
长寿命、高可靠性的产品在正常应力下工作时，其失效和故障时间的数据很少，并且大都是截尾数据，因此采用统计分析方法对电机进行可靠性评估时，常因历史样本数据太小而造成评估结果不准确。对电机进行加速退化试验，可以在较短的时间内获取更多有用的退化数据，有效提高试验效率及预测准确性。长寿命、高可靠性的产品性能由于损伤累加性，在退化不同阶段的变化具有不同的统计特征，其退化过程往往呈多阶段性，如平稳退化、快速退化、故障三个阶段。因此，研究退化过程中的变点估计问题具有很大的工程价值。目前，已存在一些基于加速退化试验进行寿命预测的研究，但未考虑实际产品退化过程表现出的多阶段性。
不同产品个体由于异质性，其退化过程不可能完全一致。对产品个体进行剩余寿命预测时，现有的研究方法主要是根据历史退化数据获取参数的先验估计值，再融合个体的现场退化数据，获取参数的后验估计值。在多阶段退化过程中，如何根据变点的先验概率分布函数对个体在线检测退化过程中估计变点出现的概率缺乏研究。对长寿命、高可靠性产品进行寿命预测，关键在于对退化过程变点的估计以及模型的建立，而根据历史数据对多阶段退化过程建模并进行变点估计，因其统计特征差，已不再适用于高可靠性产品领域。
发明内容
鉴于此，本发明提供一种融合现场数据和两阶段加速退化数据的寿命预测方法，解决了历史失效故障数据少、收集困难以及利用加速退化数据建模中未考虑含变点的多阶段退化过程对可靠性评估的不准确性问题。
为实现上述目的，本发明提供一种融合现场数据和两阶段加速退化数据的寿命预测方法，所述方法包括步骤：
S1、选取N个电机绝缘材料作为试验样本进行恒应力加速退化试验，设定n个加速温度应力，将N个电机绝缘材料分为n组，记录第k个加速温度应力下对应n k 个试验样本，在每一个加速温度应力下对相应的试验样本进行恒应力加速退化试验，得到N组加速退化数据，每一组加速退化数据对应为该电机绝缘材料在第k个加速温度应力下并基于各个测量时刻下所得到的加速退化数据序列，k取值为1~n；
S2、根据每一组加速退化数据建立两阶段Wiener过程退化模型，通过SIC确定每一个电机绝缘材料退化过程中变点的区间，基于残差平方和最小准则得到变点的估计值，进而得到了第k个加速温度应力下n k 个电机绝缘材料的n k 个变点位置，假设每一个加速温度应力下的变点位置服从两参数威布尔分布，对第k个加速温度应力下的n k 个变点位置利用极大似然估计法得到该加速温度应力下的变点分布参数，进而得到n组变点分布参数；
S3、根据变点分布参数和加速温度应力之间的变化规律构建变点分布参数的加速模型，建立变点分布参数与加速温度应力之间的矩阵方程，通过求解矩阵方程得到最小二乘估计值，并根据所述最小二乘估计值得到在正常温度应力下两阶段退化过程中的变点分布参数值，并得到正常温度应力下变点的概率密度函数；
S4、对N组加速退化数据构建两阶段Wiener过程退化模型, 在Wiener过程第一阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数，求解得到第一阶段的极大似然估计值，在Wiener过程第二阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数，求解得到第二阶段的极大似然估计值；
S5、根据第一阶段的极大似然估计值和第二阶段的极大似然估计值，得到正常温度应力下两阶段Wiener过程退化模型中两阶段的漂移参数和扩散参数，结合正常温度应力下变点的概率密度函数，并根据预设的性能参数失效阈值，得到正常温度应力下产品寿命的累积分布函数；
S6、获取待测电机绝缘材料在正常温度应力下测量得到的现场退化数据，根据正常温度应力下变点的概率密度函数计算现场退化数据出现变点的置信度，以及根据CUSUM控制图辨识退化过程变点，判断现场退化数据中是否出现变点；
S7、若现场退化数据中出现变点，根据两阶段维纳过程退化模型中漂移参数和扩散参数服从共轭先验分布，将加速退化数据作为先验信息，利用变点前的现场退化数据计算第一阶段的漂移参数和扩散参数的后验期望值，利用变点后的现场退化数据计算第二阶段的漂移参数和扩散参数的后验期望值；
S8、根据第一阶段的漂移参数和扩散参数的后验期望值、以及第二阶段的漂移参数和扩散参数的后验期望值更新正常温度应力下两阶段维纳过程退化模型中两阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品寿命的累积分布函数。
优选的，所述步骤S7包括：
若现场退化数据中变点未出现，利用现场退化数据计算第一阶段的漂移参数和扩散参数的后验期望值，更新正常温度应力下两阶段维纳过程退化模型第一阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品的累积分布函数。
优选的，所述步骤S2包括：
根据加速退化数据建立两阶段Wiener过程退化模型，并建立似然函数，利用SIC确定电机绝缘材料退化过程中变点的区间，SIC定义为：
Figure 345599DEST_PATH_IMAGE001
；
其中，
Figure 571568DEST_PATH_IMAGE002
是退化模型的极大似然函数，
Figure 135405DEST_PATH_IMAGE003
是
Figure 709474DEST_PATH_IMAGE004
的极大似然估计，p是退化模型 中的参数个数，m是样本的大小，设在第k个加速温度应力下，第j个电机绝缘材料实验样本 的变点
Figure 751380DEST_PATH_IMAGE005
位于时间间隔
Figure 891243DEST_PATH_IMAGE006
内，建立退化增量的似然函数；
根据SIC假设检验，H 0 表示无变点，H 1 表示有变点，则，
Figure 625981DEST_PATH_IMAGE007
Figure 156188DEST_PATH_IMAGE009
；
式中，
Figure 64101DEST_PATH_IMAGE010
；
计算原假设H 0 下的SIC(n jk )值以及备择假设H 1下的SIC(a jk )值；如果
Figure 746886DEST_PATH_IMAGE011
，则拒绝原假设H 0，即接受H 1，认为存在变点，变点的估计 区间位于
Figure 904723DEST_PATH_IMAGE012
，
Figure 735276DEST_PATH_IMAGE013
；
第k个加速温度应力下，第j个产品的变点
Figure 102672DEST_PATH_IMAGE014
位于时间间隔
Figure 967860DEST_PATH_IMAGE015
内，对 两阶段Wiener过程退化模型利用残差平方和最小准则获取变点
Figure 778821DEST_PATH_IMAGE016
的估计值，建立关于
Figure 549200DEST_PATH_IMAGE017
的一元二次方程RSS j ，并且
Figure 736598DEST_PATH_IMAGE018
，根据方程性质求出当RSS j 取得最小值时， 得到变点的取值为
Figure 908823DEST_PATH_IMAGE019
。
优选的，所述步骤S3包括：
设T k 为第个k加速温度应力，求出T k 下n k 个电机绝缘材料的变点估计值
Figure 156264DEST_PATH_IMAGE020
，假设变点服从两参数威布尔分布，则变点的概率密度函数 为：
Figure 428588DEST_PATH_IMAGE021
；
式中，m为形状参数，η为尺度参数；
利用极大似然法进行参数估计，构建参数的似然函数，求解得到形状参数m，尺度 参数η的估计值
Figure 419677DEST_PATH_IMAGE022
，得到n个加速温度应力下变点服从两参数威布尔的参数估计值
Figure 243146DEST_PATH_IMAGE023
；
变点分布参数与加速温度应力之间的变化规律符合Arrhenius模型，第k个加速温度应力T k 下的形状参数m表示为：
Figure 395910DEST_PATH_IMAGE024
；
第k个加速温度应力
Figure 875301DEST_PATH_IMAGE025
下的尺度参数
Figure 466820DEST_PATH_IMAGE026
表示为：
Figure 551319DEST_PATH_IMAGE027
；
式中，
Figure 140564DEST_PATH_IMAGE028
为待定系数；
建立变点分布参数与加速温度应力之间的矩阵方程为：
Figure 641340DEST_PATH_IMAGE029
；
求解该矩阵方程得到最小二乘估计值
Figure 708653DEST_PATH_IMAGE030
，根据最小二乘估计值外 推出正常温度应力T 0下两阶段退化过程中的形状参数m 0、尺度参数η 0的参数值为：
Figure 991867DEST_PATH_IMAGE031
；
正常温度应力T 0下变点的概率密度函数
Figure 266859DEST_PATH_IMAGE032
表示为：
Figure 205996DEST_PATH_IMAGE033
。
优选的，所述步骤S4包括：
在Wiener过程第一阶段中，第k个加速温度应力T k 下的漂移参数
Figure 326268DEST_PATH_IMAGE034
表示为：
Figure 932830DEST_PATH_IMAGE035
；
第k个加速温度应力
Figure 378723DEST_PATH_IMAGE036
下的扩散参数
Figure 805157DEST_PATH_IMAGE037
表示为：
Figure 794366DEST_PATH_IMAGE038
：
式中，
Figure 989855DEST_PATH_IMAGE039
为待定系数；
满足加速因子不变原则，则
Figure 888541DEST_PATH_IMAGE040
，漂移参数
Figure 51538DEST_PATH_IMAGE041
和扩散参数
Figure 326661DEST_PATH_IMAGE042
在Wiener过 程第一阶段的加速退化模型为：
Figure 625924DEST_PATH_IMAGE043
；
根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数为：
Figure 961091DEST_PATH_IMAGE044
；
式中，
Figure 362116DEST_PATH_IMAGE045
表示退化增量，
Figure 627881DEST_PATH_IMAGE046
表示时间增 量；
求解得到第一阶段的极大似然估计值
Figure 860280DEST_PATH_IMAGE047
。
优选的，所述步骤S4还包括：
在Wiener过程第二阶段中漂移参数
Figure 290648DEST_PATH_IMAGE048
和扩散参数
Figure 241287DEST_PATH_IMAGE049
加速退化模型为：
Figure 795896DEST_PATH_IMAGE050
；
根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数为：
Figure 132068DEST_PATH_IMAGE052
；
式中，
Figure 746720DEST_PATH_IMAGE053
表示退化增量，
Figure 106026DEST_PATH_IMAGE054
表示时间 增量；
求解得到第二阶段的极大似然估计值
Figure 526644DEST_PATH_IMAGE055
。
优选的，所述步骤S5包括：
电机绝缘材料在正常温度应力
Figure 405738DEST_PATH_IMAGE056
下，根据加速模型得到Wiener过程第一阶段在
Figure 706138DEST_PATH_IMAGE057
下的漂移参数
Figure 365789DEST_PATH_IMAGE058
和扩散参数
Figure 237101DEST_PATH_IMAGE059
，以及得到Wiener过程第二阶段在T 0下的漂移参数
Figure 954390DEST_PATH_IMAGE060
和扩散参数
Figure 176424DEST_PATH_IMAGE061
：
Figure 57792DEST_PATH_IMAGE062
；
Figure 272742DEST_PATH_IMAGE063
；
得到正常温度应力
Figure 844537DEST_PATH_IMAGE064
下电机绝缘材料在首次达到失效阈值D时产品寿命的累积 分布函数：
Figure 237472DEST_PATH_IMAGE066
；
优选的，所述步骤S6包括：
获取与试验样本同种型号的待测电机绝缘材料，测量得到在正常温度应力下的现 场退化数据
Figure 327176DEST_PATH_IMAGE067
；
根据正常温度应力下变点的概率密度函数
Figure 627707DEST_PATH_IMAGE068
，确定现场退化数据中变点出现 的置信度C；
Figure 335900DEST_PATH_IMAGE069
；
根据CUSUM控制图检测现场退化数据中是否存在变点：计算现场退化数据的均值
Figure 680163DEST_PATH_IMAGE070
，将现场退化数据与均值的差累加到累计和上；
Figure 739385DEST_PATH_IMAGE071
；
其中，
Figure 843608DEST_PATH_IMAGE072
，极差值定义为
Figure 921154DEST_PATH_IMAGE073
；
执行N次循环，对现场退化数据进行有放回重采样，重新计算新序列的累计和，记 为
Figure 452629DEST_PATH_IMAGE074
，计算极差值，记为
Figure 982837DEST_PATH_IMAGE075
，
Figure 828433DEST_PATH_IMAGE076
，比较
Figure 573535DEST_PATH_IMAGE075
与
Figure 465792DEST_PATH_IMAGE077
，记n为
Figure 296345DEST_PATH_IMAGE078
的次数，计算变点出现的置信度P为;
Figure 929321DEST_PATH_IMAGE079
;
将置信度P与置信度C进行比较，若P大于C，则表明现场退化数据中存在的变点的 置信度水平至少为C，变点的位置为
Figure 997771DEST_PATH_IMAGE080
。
优选的，所述步骤S7包括：
对两阶段Wiener过程，假设参数
Figure 871049DEST_PATH_IMAGE081
服从共轭先验分布：
Figure 907007DEST_PATH_IMAGE082
；
式中，
Figure 94406DEST_PATH_IMAGE083
表示Gamma分布，
Figure 814100DEST_PATH_IMAGE084
表示Normal分布，a,b,c,d为超参数；
若现场退化数据在
Figure 310809DEST_PATH_IMAGE085
处出现变点，根据加速退化试验数据作为先验数据，利用EM 算法，获取的超参数的先验估计值
Figure 788058DEST_PATH_IMAGE086
和
Figure 828083DEST_PATH_IMAGE087
；
结合变点前的现场退化数据
Figure 339967DEST_PATH_IMAGE088
计算第一阶 段超参数的后验估计值，以及变点后的现场退化数据
Figure 555047DEST_PATH_IMAGE089
计算第二阶段超参数的后验估计值，其 中，超参数的的后验估计值可由Bayes公式推导出：
Figure 300018DEST_PATH_IMAGE090
；
其中，
Figure 563641DEST_PATH_IMAGE091
、
Figure 507195DEST_PATH_IMAGE092
为似然函数，
Figure 830860DEST_PATH_IMAGE093
、
Figure 345018DEST_PATH_IMAGE094
为联合先验密度函数，
Figure 661598DEST_PATH_IMAGE095
、
Figure 413654DEST_PATH_IMAGE096
为联合后验 密度函数；
随机参数的共轭先验分布与其后验分布具有相同的形式，推导出超参数的后验估 计值
Figure 234453DEST_PATH_IMAGE097
和
Figure 173590DEST_PATH_IMAGE098
，进一步推导出两阶 段随机参数
Figure 106911DEST_PATH_IMAGE099
、
Figure 697161DEST_PATH_IMAGE100
，
Figure 628208DEST_PATH_IMAGE101
，
Figure 366226DEST_PATH_IMAGE102
的后验期望值
Figure 40921DEST_PATH_IMAGE103
；
所述步骤S8包括：
更新后的正常温度应力下产品寿命的累积分布函数：
Figure 298727DEST_PATH_IMAGE105
。
优选的，所述步骤S7还包括：
现场退化数据还未出现变点，根据加速退化试验数据作为先验数据，利用EM算法 获取的超参数的先验估计值
Figure 384364DEST_PATH_IMAGE106
，结合现场退化数据
Figure 360410DEST_PATH_IMAGE107
，超参数的的后验估计值可由Bayes公式推导出：
Figure 90993DEST_PATH_IMAGE108
；
式中，
Figure 140989DEST_PATH_IMAGE109
为似然函数，
Figure 476155DEST_PATH_IMAGE110
为联合先验密度函数，
Figure 860869DEST_PATH_IMAGE111
为联合后验密度函数，其中，
Figure 939683DEST_PATH_IMAGE112
；
随机参数的共轭先验分布与其后验分布具有相同的形式，计算得到超参数的后验 估计值
Figure 93453DEST_PATH_IMAGE113
，进一步得到随机参数
Figure 802783DEST_PATH_IMAGE114
、
Figure 940372DEST_PATH_IMAGE115
的后验期望值
Figure 557298DEST_PATH_IMAGE116
；
所述步骤S8包括：
更新后的正常温度应力下产品寿命的累积分布函数为：
Figure 847465DEST_PATH_IMAGE118
。
与现有技术相比，本发明提供的一种融合现场数据和两阶段加速退化数据的寿命预测方法，所带来的有益效果为：根据电机绝缘材料试验样本的加速退化数据，利用SIC和残差平方和最小准则确定加速温度应力下退化过程的变点位置，构建两阶段Wiener过程退化模型，再根据Arrhenius方程外推产品在正常应力水平下退化过程中变点的分布及退化模型。利用待测电机绝缘材料的现场退化数据，用CUSUM控制图检测发生变点的置信度，结合先验变点概率密度函数，确定变点发生的概率，将加速退化数据作为先验信息获得超参数的先验估计值，再基于Bayes方法根据现场退化数据得到超参数的后验估计值，推导出随机参数的后验估计值，进而得到更新后的产品累积分布函数，实现对电机绝缘材料产品的寿命预测。利用加速退化数据研究产品的可靠性大大提升了试验效率，通过变点检测对退化过程分阶段建模，提高了产品可靠性评估的准确性，对产品的寿命预测与健康管理有很好的工程应用价值。
附图说明
图1是根据本发明的一个实施例的融合现场数据和两阶段加速退化数据的寿命预测方法的流程示意图；
图2是根据本发明的一个具体实施例的外推正常应力下退化过程变点分布的流程示意图；
图3是根据本发明的一个实施例的外推正常应力下产品寿命的累积分布函数的流程示意图；
图4是根据本发明的一个实施例的融合现场退化数据的正常应力下的产品寿命的累积分布函数的流程示意图。
具体实施方式
以下将结合附图所示的具体实施方式对本发明进行详细描述，但这些实施方式并不限制本发明，本领域的普通技术人员根据这些实施方式所做出的结构、方法、或功能上的变换均包含在本发明的保护范围内。
如图1所示的本发明的一个实施例，本发明提供一种融合现场数据和两阶段加速退化数据的寿命预测方法，该方法包括步骤：
S1、选取N个电机绝缘材料作为试验样本进行恒应力加速退化试验，设定n个加速温度应力，将N个电机绝缘材料分为n组，记录第k个加速温度应力下对应n k 个试验样本，在每一个加速温度应力下对相应的试验样本进行恒应力加速退化试验，得到N组加速退化数据，每一组加速退化数据对应为该电机绝缘材料在第k个加速温度应力下并基于各个测量时刻下所得到的加速退化数据序列，k取值为1~n；
S2、根据每一组加速退化数据建立两阶段Wiener过程退化模型，通过SIC确定每一个电机绝缘材料退化过程中变点的区间，基于残差平方和最小准则得到变点的估计值，进而得到了第k个加速温度应力下n k 个电机绝缘材料的n k 个变点位置，假设每一个加速温度应力下的变点位置服从两参数威布尔分布，对第k个加速温度应力下的n k 个变点位置利用极大似然估计法得到该加速温度应力下的变点分布参数，进而得到n组变点分布参数；
S3、根据变点分布参数和加速温度应力之间的变化规律构建变点分布参数的加速模型，建立变点分布参数与加速温度应力之间的矩阵方程，通过求解矩阵方程得到最小二乘估计值，并根据所述最小二乘估计值得到在正常温度应力下两阶段退化过程中的变点分布参数值，并得到正常温度应力下变点的概率密度函数；
S4、对N组加速退化数据构建两阶段Wiener过程退化模型, 在Wiener过程第一阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数，求解得到第一阶段的极大似然估计值，在Wiener过程第二阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数，求解得到第二阶段的极大似然估计值；
S5、根据第一阶段的极大似然估计值和第二阶段的极大似然估计值，得到正常温度应力下两阶段Wiener过程退化模型中两阶段的漂移参数和扩散参数，以及得到正常温度应力下变点的概率密度函数，并根据预设的性能参数失效阈值，得到正常温度应力下产品寿命的累积分布函数；
S6、获取待测电机绝缘材料在正常温度应力下测量得到的现场退化数据，根据正常温度应力下变点的概率密度函数计算现场退化数据出现变点的置信度，以及根据CUSUM控制图辨识退化过程变点，判断现场退化数据中是否出现变点；
S7、若现场退化数据中出现变点，根据两阶段维纳过程退化模型中漂移参数和扩散参数服从共轭先验分布，将加速退化数据作为先验信息，利用变点前的现场退化数据计算第一阶段的漂移参数和扩散参数的后验期望值，利用变点后的现场退化数据计算第二阶段的漂移参数和扩散参数的后验期望值；
S8、根据第一阶段的漂移参数和扩散参数的后验期望值、以及第二阶段的漂移参数和扩散参数的后验期望值更新正常温度应力下两阶段维纳过程退化模型中两阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品寿命的累积分布函数。
选取N个电机绝缘材料作为试验样本进行恒应力加速退化试验，加速退化试验各 个不同等级的应力下进行的实验。设定n个加速温度应力，将N个电机绝缘材料分为n组，对 每一个加速温度应力下的试验样本进行试验。记录第k个加速恒温应力下对应n k 个试验样 本，例如设定4个加速温度应力，那么将N分为4组，第一个加速温度应力下有n 1 个样本，第二 个加速温度应力下有n 2 个样本，第三个加速温度应力下有n 3 个样本，第四个加速温度应力下 有n 4 个样本。在每一个加速温度应力下对相应的试验样本进行恒应力加速退化试验，即在 第k个加速温度应力下对n k 个试验样本进行应力加速退化试验，并在各个测量时刻下测量 得到表征性能退化的退化数据，该退化数据是基于各个测量时刻的加速退化数据序列，得 到N组加速退化数据。用xijk表示第k个应力下第j个试验样本第i次测量得到的数据，t ijk 表 示第k个应力下第j个试验样本第i次测量时的测量时刻，
Figure 979894DEST_PATH_IMAGE119
表示 退化增量，
Figure 89932DEST_PATH_IMAGE120
表示时间增量，其中，i=1,2,3…,n jk , j=1,2,3…,n k ，k=1,2,3…,n，n jk 表示第k个加速温度应力第j个试样样本，n k 表示第k个加速温度应力下的 试样样本数，n为加速温度应力个数。
采用两阶段Wiener过程对加速退化数据建立退化模型，退化过程X(t)存在随机变 点
Figure 759817DEST_PATH_IMAGE121
，定义两阶段Wiener过程退化模型如下：
Figure 904490DEST_PATH_IMAGE122
（1）；
其中，
Figure 939311DEST_PATH_IMAGE123
与
Figure 598963DEST_PATH_IMAGE124
分别为两阶段Wiener过程的漂移参数，
Figure 495375DEST_PATH_IMAGE125
、
Figure 806139DEST_PATH_IMAGE126
为扩散系数，
Figure 497015DEST_PATH_IMAGE127
为标准Brownian运动，退化初始量为0，由Wiener过程性质可知，退化增量
Figure 99422DEST_PATH_IMAGE128
服从正态分布
Figure 861841DEST_PATH_IMAGE129
，对独立增量进行累 加，推导得出：
Figure 715528DEST_PATH_IMAGE130
（2）；
可以得到X(t)的概率密度函数为：
Figure 888889DEST_PATH_IMAGE132
（3）；
设D为性能参数的失效阈值，则产品的寿命
Figure 460816DEST_PATH_IMAGE133
可被定义为X(t)首次达到失效阈值 的时间
Figure 761347DEST_PATH_IMAGE134
，产品的累积分布函数为如下逆高斯分布：
Figure 249966DEST_PATH_IMAGE136
（4）；
其中，
Figure 548224DEST_PATH_IMAGE137
为标准正态累积分布函数。
根据加速退化数据建立两阶段Wiener过程退化模型，通过SIC、残差平方和最小准 则确定加速温度应力下变点位置，如图2所示的流程图。由于在加速退化试验中，只能离散 地测出产品的退化量，变点
Figure 856714DEST_PATH_IMAGE138
的位置未知，因此利用SIC（Schwarz information criterion，施瓦茨信息准则）确定变点
Figure 164199DEST_PATH_IMAGE138
的位置。根据加速退化数据建立两阶段Wiener过 程退化模型，并建立似然函数，利用SIC确定电机绝缘材料退化过程中变点的区间，SIC定义 为：
Figure 498535DEST_PATH_IMAGE139
（5）
其中，
Figure 233273DEST_PATH_IMAGE140
是退化模型的极大似然函数，
Figure 825797DEST_PATH_IMAGE141
是
Figure 671393DEST_PATH_IMAGE142
的极大似然估计，p是退化模型中 的参数个数，m是样本的大小，设在第k个加速温度应力下，第j个电机绝缘材料实验样本的 变点
Figure 869025DEST_PATH_IMAGE143
位于时间间隔
Figure 509085DEST_PATH_IMAGE144
内，其退化增量的似然函数为：
Figure 792168DEST_PATH_IMAGE145
（6）；
对上式取对数，得到对数似然函数为：
Figure 441455DEST_PATH_IMAGE146
（7）；
根据SIC假设检验，H 0 表示无变点，H 1 表示有变点，则
Figure 41064DEST_PATH_IMAGE147
Figure 369801DEST_PATH_IMAGE148
；
式中，
Figure 890913DEST_PATH_IMAGE149
。
原假设H 0 下的SIC(n jk )值为：
Figure 327579DEST_PATH_IMAGE150
（8）；
式中，
Figure 984956DEST_PATH_IMAGE151
是基于原假设H 0 的似然函数，
Figure 278403DEST_PATH_IMAGE152
。
备择假设H 1下的SIC(a jk )值为：
Figure 21231DEST_PATH_IMAGE154
（9）；
式中，
Figure 809059DEST_PATH_IMAGE155
是基于备择假设H 1的似然函数，其中，
Figure 570210DEST_PATH_IMAGE156
；
Figure 785291DEST_PATH_IMAGE157
；
如果
Figure 533192DEST_PATH_IMAGE158
，则拒绝原假设H 0，即接受H 1，认为存在变 点，变点的估计区间位于
Figure 796814DEST_PATH_IMAGE159
，
Figure 225521DEST_PATH_IMAGE160
。
第k个加速温度应力下，第j个产品的变点
Figure 860771DEST_PATH_IMAGE161
位于时间间隔
Figure 578191DEST_PATH_IMAGE159
内，对两 阶段Wiener过程退化模型利用残差平方和最小准则获取变点
Figure 894772DEST_PATH_IMAGE161
的估计值，
Figure 646827DEST_PATH_IMAGE163
（10）；
式中，
Figure 656240DEST_PATH_IMAGE164
，
Figure 657694DEST_PATH_IMAGE165
，
Figure 528698DEST_PATH_IMAGE166
，
Figure 184195DEST_PATH_IMAGE167
；
式（12）是关于
Figure 849663DEST_PATH_IMAGE161
的一元二次方程，并且
Figure 338413DEST_PATH_IMAGE168
，根据方程性质 即可求出当RSS j 取得最小值时，变点的取值为
Figure 996796DEST_PATH_IMAGE169
，即为变点的估计值。进而得到了第k个 加速温度应力下n k 个电机绝缘材料的n k 个变点位置。比如，由上述实施例可知，第一个加速 温度应力下n 1 个样本有n 1 个变点位置，第二个加速温度应力下n 2 个样本有n 2 个变点位置，第 三个加速温度应力下n3个样本有n 3 个变点位置，第四个加速温度应力下n 4 个样本有n 4 个变 点位置。因此可以得到每一个加速温度应力下的每一个电机绝缘材料的变点位置，共得到n*n k 个变点位置。
设
Figure 254602DEST_PATH_IMAGE170
为第个k加速温度应力，根据上述方案，可以求出
Figure 605818DEST_PATH_IMAGE170
下
Figure 581864DEST_PATH_IMAGE171
个电机绝缘材 料的变点估计值
Figure 529092DEST_PATH_IMAGE172
，假设变点服从两参数威布尔分布，则变点的 概率密度函数为：
Figure 156251DEST_PATH_IMAGE173
（11）；
式中，m为形状参数，η为尺度参数，利用极大似然法进行参数估计，则其似然函数为：
Figure 163521DEST_PATH_IMAGE174
（12）；
对上式取对数，得到对数似然函数为：
Figure 626864DEST_PATH_IMAGE175
（13）；
将上式分别对求m、η偏导，并令偏导数等于0，建立方程组：
Figure 641698DEST_PATH_IMAGE177
（14）；
求解上式可以得到形状参数m，尺度参数η的估计值
Figure 874096DEST_PATH_IMAGE178
，根据该方案，由上述 实施例可知，第一个加速温度应力下
Figure 567115DEST_PATH_IMAGE179
个样本的
Figure 252174DEST_PATH_IMAGE179
个变点位置估计得到一组变点的分 布参数
Figure 72362DEST_PATH_IMAGE180
，第二个加速温度应力下
Figure 346218DEST_PATH_IMAGE181
个样本的
Figure 23187DEST_PATH_IMAGE181
个变点位置估计得到一组变点 的分布参数
Figure 382493DEST_PATH_IMAGE182
，第三个加速温度应力下
Figure 537531DEST_PATH_IMAGE183
个样本
Figure 682204DEST_PATH_IMAGE183
个变点位置估计得到一组 变点的分布参数
Figure 782272DEST_PATH_IMAGE184
，第四个加速温度应力下有
Figure 379606DEST_PATH_IMAGE185
个样本
Figure 338335DEST_PATH_IMAGE185
个变点位置估计得到 一组变点的分布参数
Figure 321203DEST_PATH_IMAGE186
，可以得到n个加速温度应力下变点服从两参数威布尔的参 数估计值
Figure 339975DEST_PATH_IMAGE187
。
根据变点分布参数和加速温度应力之间的变化规律构建变点分布参数的加速模 型。变点分布参数与加速温度应力之间的变化规律符合Arrhenius模型，第k个加速温度应 力
Figure 673873DEST_PATH_IMAGE188
下的形状参数m表示为：
Figure 436293DEST_PATH_IMAGE189
（15）；
第k个加速温度应力T k 下的尺度参数η表示为：
Figure 289979DEST_PATH_IMAGE190
（16）；
式中，
Figure 463341DEST_PATH_IMAGE191
为待定系数，建立变点分布参数与加速温度应力之 间的矩阵方程为：
Figure 300847DEST_PATH_IMAGE193
（17）；
求解该矩阵方程得到最小二乘估计值
Figure 791258DEST_PATH_IMAGE194
，根据最小二乘估计值 由此可以外推出正常温度应力
Figure 827347DEST_PATH_IMAGE195
下，两阶段退化过程中的变点服从威布尔参数，即形状 参数
Figure 860026DEST_PATH_IMAGE196
、尺度参数
Figure 230833DEST_PATH_IMAGE197
的参数值为：
Figure 538317DEST_PATH_IMAGE198
（18）；
正常温度应力
Figure 350285DEST_PATH_IMAGE199
下变点的概率密度函数
Figure 881760DEST_PATH_IMAGE200
表示为：
Figure 162700DEST_PATH_IMAGE201
（19）；
对N组加速退化数据构建两阶段Wiener过程退化模型，对退化模型中的模型参数 进行估计，如图3所示的流程图。设定电机绝缘材料在加速温度应力
Figure 319880DEST_PATH_IMAGE202
、加速温度应力
Figure 2666DEST_PATH_IMAGE203
下，定义
Figure 705042DEST_PATH_IMAGE202
相当于
Figure 991055DEST_PATH_IMAGE203
的加速因子为
Figure 374763DEST_PATH_IMAGE204
，根据加速因子不变原则，两阶段Wiener过程退 化模型的漂移参数
Figure 223639DEST_PATH_IMAGE205
和扩散参数
Figure 34600DEST_PATH_IMAGE206
满足以下关系式：
Figure 352449DEST_PATH_IMAGE207
（20）；
在加速温度应力下，两阶段Wiener过程退化模型的漂移参数
Figure 523536DEST_PATH_IMAGE208
和扩散参数
Figure 243231DEST_PATH_IMAGE209
与加 速温度应力之间的变化规律用Arrhenius模型描述。
在Wiener过程第一阶段中，第k个加速温度应力
Figure 474361DEST_PATH_IMAGE210
下的漂移参数
Figure 279506DEST_PATH_IMAGE211
表示为：
Figure 5016DEST_PATH_IMAGE212
（21）；
第k个加速温度应力
Figure 819695DEST_PATH_IMAGE213
下的扩散参数
Figure 238038DEST_PATH_IMAGE214
表示为：
Figure 451851DEST_PATH_IMAGE215
（22）；
式中，
Figure 43369DEST_PATH_IMAGE216
为待定系数，要满足加速因子不变原则，则
Figure 675339DEST_PATH_IMAGE217
；因此漂 移参数
Figure 248272DEST_PATH_IMAGE218
和扩散参数
Figure 231271DEST_PATH_IMAGE219
在Wiener过程第一阶段的加速退化模型为：
Figure 344589DEST_PATH_IMAGE220
（23）；
根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数为：
Figure 627803DEST_PATH_IMAGE222
（24）；
求解得到第一阶段的极大似然估计值
Figure 653528DEST_PATH_IMAGE223
。
类似地，在Wiener过程第二阶段中漂移参数
Figure 844863DEST_PATH_IMAGE224
和扩散参数
Figure 778183DEST_PATH_IMAGE225
加速退化模型为：
Figure 102854DEST_PATH_IMAGE226
（25）；
根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数为：
Figure 96218DEST_PATH_IMAGE228
（26）；
求解得到第二阶段的极大似然估计值
Figure 522651DEST_PATH_IMAGE229
。
电机绝缘材料在正常温度应力
Figure 446614DEST_PATH_IMAGE230
下，根据加速模型得到Wiener过程第一阶段在T 0下的漂移参数
Figure 438841DEST_PATH_IMAGE231
和扩散参数
Figure 55636DEST_PATH_IMAGE232
，以及得到Wiener过程第二阶段在
Figure 766103DEST_PATH_IMAGE230
下的漂移参数
Figure 978909DEST_PATH_IMAGE233
和扩散参数
Figure 608998DEST_PATH_IMAGE234
：
Figure 616269DEST_PATH_IMAGE235
（27）；
Figure 79611DEST_PATH_IMAGE236
（28）；
根据上述步骤，得到正常温度应力
Figure 79797DEST_PATH_IMAGE237
下电机绝缘材料在首次达到失效阈值D时产 品寿命的累积分布函数：
Figure 312195DEST_PATH_IMAGE239
（29）；
可靠度函数为：
Figure 739634DEST_PATH_IMAGE240
（30）。
在正常温度应力下，基于加速退化先验信息对待测的电机绝缘材料进行寿命预 测。获取与试验样本同种型号的待测电机绝缘材料，测量得到在正常温度应力下的现场退 化数据
Figure 690273DEST_PATH_IMAGE241
。根据正常温度应力下变点的概率密度函数
Figure 244882DEST_PATH_IMAGE242
，确定现场退化数据中变点出现的置信度C：
Figure 581054DEST_PATH_IMAGE243
（31）；
根据CUSUM控制图辨识退化过程变点，判断现场退化数据中是否存在变点，具体 地，计算现场退化数据的均值
Figure 195706DEST_PATH_IMAGE244
，将现场退化数据与均值的差累加到累计和上：
Figure 368062DEST_PATH_IMAGE245
（32）；
其中，
Figure 712980DEST_PATH_IMAGE246
，极差值定义为
Figure 919970DEST_PATH_IMAGE247
，
执行N次循环，对现场退化数据进行有放回重采样，重新计算新序列的累计和，记 为
Figure 689212DEST_PATH_IMAGE248
，计算极差值，记为
Figure 348864DEST_PATH_IMAGE249
，
Figure 245275DEST_PATH_IMAGE250
，比较
Figure 556040DEST_PATH_IMAGE251
与
Figure 512495DEST_PATH_IMAGE252
，记n为
Figure 659442DEST_PATH_IMAGE253
的次数，计算变点出现的置信度P为：
Figure 343233DEST_PATH_IMAGE254
（33）；
将置信度P与置信度C进行比较，若P大于C，则表明现场退化数据中存在的变点的 置信度水平至少为C，变点的位置为
Figure 259237DEST_PATH_IMAGE255
。
不同产品个体的退化过程不可能完全一致，为描述不同产品个体之间的异质性， 可将两阶段Wiener过程退化模型的参数服从某种分布的随机参数。对两阶段Wiener过程退 化模型来说，假设
Figure 650509DEST_PATH_IMAGE256
服从如下的共轭先验分布：
Figure 284752DEST_PATH_IMAGE257
（34）；
式中，
Figure 522967DEST_PATH_IMAGE258
表示Gamma分布，
Figure 542744DEST_PATH_IMAGE259
表示Normal分布，a,b,c,d为超参数。
若现场退化数据在
Figure 841002DEST_PATH_IMAGE260
处出现变点，可基于Bayes方法对两阶段Wiener过程退化模 型的参数进行更新。根据两阶段维纳过程退化模型中漂移参数和扩散参数服从共轭先验分 布，将加速退化数据作为先验信息。根据加速退化试验数据作为先验数据，利用EM算法，获 取的超参数的先验估计值
Figure 962541DEST_PATH_IMAGE261
和
Figure 988135DEST_PATH_IMAGE262
，结合变点前的现场退化数据
Figure 878731DEST_PATH_IMAGE263
计算第一阶段超参数的后验估计值，利用变点 后的现场退化数据
Figure 597157DEST_PATH_IMAGE264
计算第二阶段超 参数的后验估计值。超参数的的后验估计值可由Bayes公式推导出：
Figure 940414DEST_PATH_IMAGE266
（35）；
式中，
Figure 786010DEST_PATH_IMAGE267
、
Figure 783309DEST_PATH_IMAGE268
为似然函数，
Figure 423369DEST_PATH_IMAGE269
、
Figure 253922DEST_PATH_IMAGE270
为联合先验密度函数，
Figure 621318DEST_PATH_IMAGE271
、
Figure 486506DEST_PATH_IMAGE272
为联合 后验密度函数。
因为随机参数的共轭先验分布与其后验分布具有相同的形式，推导出超参数的后 验估计值
Figure 546735DEST_PATH_IMAGE273
和
Figure 864584DEST_PATH_IMAGE274
，进一步推导出 两阶段随机参数
Figure 51982DEST_PATH_IMAGE275
、
Figure 958627DEST_PATH_IMAGE276
，
Figure 2807DEST_PATH_IMAGE277
，
Figure 997832DEST_PATH_IMAGE278
的后验期望值
Figure 785660DEST_PATH_IMAGE279
。
根据第一阶段的漂移参数和扩散参数的后验期望值、以及第二阶段的漂移参数和扩散参数的后验期望值更新正常温度应力下两阶段维纳过程退化模型中两阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品寿命的累积分布函数：
Figure 297543DEST_PATH_IMAGE281
（36） 。
若现场退化数据中变点未出现，则可基于Bayes方法利用现场退化数据对Wiener 过程第一阶段的参数值进行更新。根据加速退化试验数据作为先验数据，利用EM算法获取 的超参数的先验估计值
Figure 496312DEST_PATH_IMAGE282
，结合现场退化数据
Figure 726437DEST_PATH_IMAGE283
，超参数的的后验估计值可由Bayes公式推导出：
Figure 317955DEST_PATH_IMAGE285
（38）；
式中，
Figure 933613DEST_PATH_IMAGE286
为似然函数，
Figure 319595DEST_PATH_IMAGE287
为联合先验密度函数，
Figure 286283DEST_PATH_IMAGE288
为联合后验密度函数，
Figure 415913DEST_PATH_IMAGE289
（39）；
因为随机参数的共轭先验分布与其后验分布具有相同的形式，可推导出超参数的 后验估计值
Figure 636810DEST_PATH_IMAGE290
，进一步可推导出随机参数
Figure 711470DEST_PATH_IMAGE275
、
Figure 650607DEST_PATH_IMAGE276
的后验 期望值
Figure 318348DEST_PATH_IMAGE291
，更新后的正常温度应力下产品寿命的累积分布函数 为：
Figure 908599DEST_PATH_IMAGE293
（40）。
基于更新后的正常温度应力下产品寿命的累积分布函数，该累积分布函数用以计算产品寿命在某个时刻失效的概率，从而实现对电机绝缘材料产品的寿命预测。参考附图4的流程实现过程。
尽管为示例目的，已经公开了本发明的优选实施方式，但是本领域的普通技术人员将意识到，在不脱离由所附的权利要求书公开的本发明的范围和精神的情况下，各种改进、增加以及取代是可能的。

Claims (9)

1.一种融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述方法包括步骤：
S1、选取N个电机绝缘材料作为试验样本进行恒应力加速退化试验，设定n个加速温度应力，将N个电机绝缘材料分为n组，记录第k个加速温度应力下对应n k 个试验样本，在每一个加速温度应力下对相应的试验样本进行恒应力加速退化试验，得到N组加速退化数据，每一组加速退化数据对应为该电机绝缘材料在第k个加速温度应力下并基于各个测量时刻下所得到的加速退化数据序列，k取值为1~n；
S2、根据每一组加速退化数据建立两阶段Wiener过程退化模型，通过SIC确定每一个电机绝缘材料退化过程中变点的区间，基于残差平方和最小准则得到变点的估计值，进而得到了第k个加速温度应力下n k 个电机绝缘材料的n k 个变点位置，假设每一个加速温度应力下的变点位置服从两参数威布尔分布，对第k个加速温度应力下的n k 个变点位置利用极大似然估计法得到该加速温度应力下的变点分布参数，进而得到n组变点分布参数；
S3、根据变点分布参数和加速温度应力之间的变化规律构建变点分布参数的加速模型，建立变点分布参数与加速温度应力之间的矩阵方程，通过求解矩阵方程得到最小二乘估计值，并根据所述最小二乘估计值得到在正常温度应力下两阶段退化过程中的变点分布参数值，并得到正常温度应力下变点的概率密度函数；
S4、对N组加速退化数据构建两阶段Wiener过程退化模型, 在Wiener过程第一阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数，求解得到第一阶段的极大似然估计值，在Wiener过程第二阶段中，构建漂移参数和扩散参数与加速温度应力之间的加速退化模型，根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数，求解得到第二阶段的极大似然估计值；
S5、根据第一阶段的极大似然估计值和第二阶段的极大似然估计值，得到正常温度应力下两阶段Wiener过程退化模型中两阶段的漂移参数和扩散参数，结合正常温度应力下变点的概率密度函数，并根据预设的性能参数失效阈值，得到正常温度应力下产品寿命的累积分布函数；
S6、获取待测电机绝缘材料在正常温度应力下测量得到的现场退化数据，根据正常温度应力下变点的概率密度函数计算现场退化数据出现变点的置信度，以及根据CUSUM控制图辨识退化过程变点，判断现场退化数据中是否出现变点；
S7、若现场退化数据中出现变点，根据两阶段维纳过程退化模型中漂移参数和扩散参数服从共轭先验分布，将加速退化数据作为先验信息，利用变点前的现场退化数据计算第一阶段的漂移参数和扩散参数的后验期望值，利用变点后的现场退化数据计算第二阶段的漂移参数和扩散参数的后验期望值；
S8、根据第一阶段的漂移参数和扩散参数的后验期望值、以及第二阶段的漂移参数和扩散参数的后验期望值更新正常温度应力下两阶段维纳过程退化模型中两阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品寿命的累积分布函数;
其中，所述步骤S3包括：
设T k 为第k个加速温度应力，求出T k 下n k 个电机绝缘材料的变点估计值
Figure 105145DEST_PATH_IMAGE001
，假设变点服从两参数威布尔分布，则变点的概率密度函数为：
Figure 959968DEST_PATH_IMAGE002
；
式中，m为形状参数，η为尺度参数；
利用极大似然法进行参数估计，构建参数的似然函数，求解得到形状参数m，尺度参数η 的估计值
Figure 929586DEST_PATH_IMAGE003
，得到n个加速温度应力下变点服从两参数威布尔的参数估计值
Figure 528057DEST_PATH_IMAGE004
；
变点分布参数与加速温度应力之间的变化规律符合Arrhenius模型，第k个加速温度应力T k 下的形状参数m表示为：
Figure 529380DEST_PATH_IMAGE005
；
第k个加速温度应力下T k 的尺度参数η表示为：
Figure 617422DEST_PATH_IMAGE006
；
式中，a 1、a 2、a 3、a 4为待定系数；
建立变点分布参数与加速温度应力之间的矩阵方程为：
Figure 87718DEST_PATH_IMAGE007
；
求解该矩阵方程得到最小二乘估计值
Figure 473568DEST_PATH_IMAGE008
，根据最小二乘估计值外推出 正常温度应力T 0 下两阶段退化过程中的形状参数m 0、尺度参数η 0的参数值为：
Figure 80130DEST_PATH_IMAGE009
；
正常温度应力T 0下变点的概率密度函数
Figure 463707DEST_PATH_IMAGE010
表示为：
Figure 421299DEST_PATH_IMAGE011
。
2.如权利要求1所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S7包括：
若现场退化数据中变点未出现，利用现场退化数据计算第一阶段的漂移参数和扩散参数的后验期望值，更新正常温度应力下两阶段维纳过程退化模型第一阶段的漂移参数和扩散参数，并得到更新后的正常温度应力下产品的累积分布函数。
3.如权利要求2所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S2包括：
根据加速退化数据建立两阶段Wiener过程退化模型，并建立似然函数，利用SIC确定电机绝缘材料退化过程中变点的区间，SIC定义为：
Figure 336472DEST_PATH_IMAGE012
；
其中，
Figure 125437DEST_PATH_IMAGE013
是退化模型的极大似然函数，
Figure 758543DEST_PATH_IMAGE014
是
Figure 390382DEST_PATH_IMAGE015
的极大似然估计，p是退化模型中的 参数个数，m是样本的大小，设在第k个加速温度应力下，第j个电机绝缘材料实验样本的变 点t jk 位于时间间隔
Figure 868768DEST_PATH_IMAGE016
内，建立退化增量的似然函数；
根据SIC假设检验，H 0 表示无变点，H 1 表示有变点，则，
Figure 699189DEST_PATH_IMAGE017
Figure 299935DEST_PATH_IMAGE018
；
式中，
Figure 966540DEST_PATH_IMAGE019
；
计算原假设H 0 下的SIC(n jk )值以及备择假设H 1下的SIC(a jk )值；如果
Figure 497884DEST_PATH_IMAGE020
，则拒绝原假设H 0，即接受H 1，认为存在变点，变点的 估计区间位于
Figure 136807DEST_PATH_IMAGE016
，
Figure 98334DEST_PATH_IMAGE021
；
第k个加速温度应力下，第j个产品的变点t jk 位于时间间隔
Figure 252235DEST_PATH_IMAGE016
内，对两阶段 Wiener过程退化模型利用残差平方和最小准则获取变点t jk 的估计值，建立关于t jk 的一元 二次方程RSS j ，并且
Figure 400320DEST_PATH_IMAGE022
，根据方程性质求出当RSS j 取得最小值时，得到变点 的取值为
Figure 205333DEST_PATH_IMAGE023
。
4.如权利要求3所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S4包括：
在Wiener过程第一阶段中，第k个加速温度应力T k 下的漂移参数m 1k 表示为：
Figure 288827DEST_PATH_IMAGE024
；
第k个加速温度应力T k 下的扩散参数σ 1k 表示为：
Figure 179292DEST_PATH_IMAGE025
；
式中，
Figure 865488DEST_PATH_IMAGE026
为待定系数；
满足加速因子不变原则，则
Figure 275741DEST_PATH_IMAGE027
，漂移参数
Figure 841720DEST_PATH_IMAGE028
和扩散参数
Figure 907896DEST_PATH_IMAGE029
在Wiener过程 第一阶段的加速退化模型为：
Figure 587664DEST_PATH_IMAGE030
；
根据每一个电机绝缘材料变点前的所有加速退化数据建立第一阶段的似然函数为：
Figure 586844DEST_PATH_IMAGE031
；
式中，
Figure 136774DEST_PATH_IMAGE032
表示退化增量，
Figure 1830DEST_PATH_IMAGE033
表示时间增 量；
求解得到第一阶段的极大似然估计值
Figure 967512DEST_PATH_IMAGE034
。
5.如权利要求4所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S4还包括：
在Wiener过程第二阶段中漂移参数
Figure 539308DEST_PATH_IMAGE035
和扩散参数
Figure 932243DEST_PATH_IMAGE036
加速退化模型为：
Figure 97645DEST_PATH_IMAGE037
；
根据每一个电机绝缘材料变点后的所有加速退化数据建立第二阶段的似然函数为：
Figure 850707DEST_PATH_IMAGE039
；
式中，
Figure 558900DEST_PATH_IMAGE040
表示退化增量，
Figure 374933DEST_PATH_IMAGE041
表示时间增 量；
求解得到第二阶段的极大似然估计值
Figure 699735DEST_PATH_IMAGE042
。
6.如权利要求5所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S5包括：
电机绝缘材料在正常温度应力T 0下，根据加速模型得到Wiener过程第一阶段在T 0下的 漂移参数
Figure 335116DEST_PATH_IMAGE043
和扩散参数
Figure 678242DEST_PATH_IMAGE044
，以及得到Wiener过程第二阶段在T 0下的漂移参数
Figure 616242DEST_PATH_IMAGE045
和扩散参 数
Figure 677608DEST_PATH_IMAGE046
：
Figure 116679DEST_PATH_IMAGE047
；
Figure 330623DEST_PATH_IMAGE048
；
得到正常温度应力T 0下电机绝缘材料在首次达到失效阈值D时产品寿命的累积分布函数：
Figure 688792DEST_PATH_IMAGE049
。
7.如权利要求6所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S6包括：
获取与试验样本同种型号的待测电机绝缘材料，测量得到在正常温度应力下的现场退 化数据
Figure 988186DEST_PATH_IMAGE050
；
根据正常温度应力下变点的概率密度函数
Figure 965369DEST_PATH_IMAGE051
，确定现场退化数据中变点出现的置 信度C：
Figure 844656DEST_PATH_IMAGE052
；
根据CUSUM控制图检测现场退化数据中是否存在变点：计算现场退化数据的均值
Figure 186776DEST_PATH_IMAGE053
，将 现场退化数据与均值的差累加到累计和上；
Figure 160417DEST_PATH_IMAGE054
；
其中，
Figure 613395DEST_PATH_IMAGE055
，
Figure 785619DEST_PATH_IMAGE056
，极差值定义为
Figure 360957DEST_PATH_IMAGE057
，
执行N次循环，对现场退化数据进行有放回重采样，重新计算新序列的累计和，记为
Figure 369365DEST_PATH_IMAGE058
，计算极差值，记为
Figure 875301DEST_PATH_IMAGE059
，
Figure 121606DEST_PATH_IMAGE060
，比较
Figure 867845DEST_PATH_IMAGE061
与
Figure 350166DEST_PATH_IMAGE062
，记n为
Figure 410526DEST_PATH_IMAGE063
的次数，计算变点出现的置信度P为：
Figure 495026DEST_PATH_IMAGE064
；
将置信度P与置信度C进行比较，若P大于C，则表明现场退化数据中存在的变点的置信 度水平至少为C，变点的位置为
Figure 349849DEST_PATH_IMAGE065
。
8.如权利要求7所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S7包括：
对两阶段Wiener过程，假设参数
Figure 660745DEST_PATH_IMAGE066
、
Figure 508484DEST_PATH_IMAGE067
服从共轭先验分布；
Figure 994960DEST_PATH_IMAGE068
；
式中，
Figure 473215DEST_PATH_IMAGE069
表示Gamma分布，
Figure 943511DEST_PATH_IMAGE070
表示Normal分布，a,b,c,d为超参数；
若现场退化数据在
Figure 407990DEST_PATH_IMAGE071
处出现变点，根据加速退化试验数据作为先验数据，利用EM算法， 获取的超参数的先验估计值
Figure 1170DEST_PATH_IMAGE072
和
Figure 463375DEST_PATH_IMAGE073
；
结合变点前的现场退化数据
Figure 607918DEST_PATH_IMAGE074
计算第一阶段超参数的 后验估计值，以及变点后的现场退化数据
Figure 548192DEST_PATH_IMAGE075
计算第二阶段超参数的后验估计值，其中，超参数的后验估计值根据Bayes公式推导出：
Figure 337156DEST_PATH_IMAGE076
；
其中，
Figure 219530DEST_PATH_IMAGE077
、
Figure 133260DEST_PATH_IMAGE078
为似然函数，
Figure 64176DEST_PATH_IMAGE079
、
Figure 910909DEST_PATH_IMAGE080
为联合先验密度函数，
Figure 511654DEST_PATH_IMAGE081
、
Figure 696036DEST_PATH_IMAGE082
为联合 后验密度函数；
随机参数的共轭先验分布与其后验分布具有相同的形式，推导出超参数的后验估计值
Figure 915796DEST_PATH_IMAGE083
和
Figure 866303DEST_PATH_IMAGE084
，进一步推导出两阶段 随机参数
Figure 637950DEST_PATH_IMAGE085
、
Figure 57430DEST_PATH_IMAGE086
，
Figure 126886DEST_PATH_IMAGE087
，
Figure 417053DEST_PATH_IMAGE088
的后验期望值
Figure 749814DEST_PATH_IMAGE089
；
所述步骤S8包括：
更新后的正常温度应力下产品寿命的累积分布函数：
Figure 125432DEST_PATH_IMAGE091
。
9.如权利要求8所述的融合现场数据和两阶段加速退化数据的寿命预测方法，其特征在于，所述步骤S7还包括：
若现场退化数据还未出现变点，根据加速退化试验数据作为先验数据，利用EM算法获 取的超参数的先验估计值
Figure 989790DEST_PATH_IMAGE092
，结合现场退化数据
Figure 727939DEST_PATH_IMAGE093
，超参数的后验估计值根据Bayes公式推导出：
Figure 779071DEST_PATH_IMAGE095
；
式中，
Figure 94515DEST_PATH_IMAGE096
为似然函数，
Figure 256506DEST_PATH_IMAGE097
为联合先验密度函数，
Figure 36112DEST_PATH_IMAGE098
为联合后验密度函数，其中，
Figure 586042DEST_PATH_IMAGE099
；
随机参数的共轭先验分布与其后验分布具有相同的形式，计算得到超参数的后验估计 值
Figure 936252DEST_PATH_IMAGE100
，进一步得到随机参数
Figure 885622DEST_PATH_IMAGE101
、
Figure 473730DEST_PATH_IMAGE102
的后验期望值
Figure 118862DEST_PATH_IMAGE103
；
所述步骤S8包括：
更新后的正常温度应力下产品寿命的累积分布函数为：
Figure 221947DEST_PATH_IMAGE105
。
CN202210068042.3A 2022-01-20 2022-01-20 一种融合现场数据和两阶段加速退化数据的寿命预测方法 Active CN114091790B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210068042.3A CN114091790B (zh) 2022-01-20 2022-01-20 一种融合现场数据和两阶段加速退化数据的寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210068042.3A CN114091790B (zh) 2022-01-20 2022-01-20 一种融合现场数据和两阶段加速退化数据的寿命预测方法

Publications (2)

Publication Number Publication Date
CN114091790A CN114091790A (zh) 2022-02-25
CN114091790B true CN114091790B (zh) 2022-05-03

Family

ID=80308961

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210068042.3A Active CN114091790B (zh) 2022-01-20 2022-01-20 一种融合现场数据和两阶段加速退化数据的寿命预测方法

Country Status (1)

Country Link
CN (1) CN114091790B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114330148B (zh) * 2022-03-11 2022-06-28 浙江大学 基于加速退化数据的电机绝缘寿命预测方法
CN116227366B (zh) * 2023-05-08 2023-08-11 浙江大学 两阶段电机绝缘寿命预测方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111046564A (zh) * 2019-12-16 2020-04-21 重庆大学 两阶段退化产品的剩余寿命预测方法
CN112683535A (zh) * 2021-01-14 2021-04-20 大连理工大学 基于多阶段维纳过程的轴承寿命预测方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107480440B (zh) * 2017-08-04 2020-01-21 山东科技大学 一种基于两阶段随机退化建模的剩余寿命预测方法
CN108959676B (zh) * 2017-12-22 2019-09-20 北京航空航天大学 一种考虑有效冲击的退化建模与寿命预测方法
WO2020216530A1 (de) * 2019-04-23 2020-10-29 Volkswagen Aktiengesellschaft Verfahren zum bestimmen von restnutzungszyklen, restnutzungszyklusbestimmungsschaltung, restnutzungszyklusbestimmungsvorrichtung
CN110569623B (zh) * 2019-09-19 2023-03-24 西安工程大学 一种细纱机电机使用寿命预测方法
CN111753416A (zh) * 2020-06-17 2020-10-09 重庆大学 一种基于两阶段Wiener过程的锂离子电池RUL预测方法
CN111814331B (zh) * 2020-07-08 2023-10-20 重庆大学 复杂情况下多变点的设备剩余使用寿命预测方法
CN112765769B (zh) * 2020-12-23 2022-07-26 北京航空航天大学 一种多阶段变工况固体推进剂剩余贮存寿命预测方法
CN112883550A (zh) * 2021-01-19 2021-06-01 中国人民解放军火箭军工程大学 一种考虑多重不确定性的退化设备剩余寿命预测方法
CN112949060B (zh) * 2021-03-01 2023-07-07 中国人民解放军火箭军工程大学 一种时变温度下的锂电池健康状态估计与剩余寿命预测方法
CN113033015B (zh) * 2021-04-09 2024-05-14 中国人民解放军火箭军工程大学 考虑两阶段自适应Wiener过程的退化设备剩余寿命预测方法
CN113761751A (zh) * 2021-09-14 2021-12-07 上海交通大学 基于温度加速因子的锂离子电池剩余寿命预测方法及系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111046564A (zh) * 2019-12-16 2020-04-21 重庆大学 两阶段退化产品的剩余寿命预测方法
CN112683535A (zh) * 2021-01-14 2021-04-20 大连理工大学 基于多阶段维纳过程的轴承寿命预测方法

Also Published As

Publication number Publication date
CN114091790A (zh) 2022-02-25

Similar Documents

Publication Publication Date Title
CN114091790B (zh) 一种融合现场数据和两阶段加速退化数据的寿命预测方法
CN106779505B (zh) 一种基于大数据驱动的输电线路故障预警方法及系统
CN108037460B (zh) 一种批产锂离子电池容量实时评估方法
CN102262191B (zh) 一种固体钽电解电容器寿命预测方法
CN105930976B (zh) 基于加权理想点法的节点电压暂降严重程度综合评估方法
CN112926144B (zh) 一种多应力加速寿命试验耦合效应分析及寿命预测方法
CN107885928B (zh) 考虑测量误差的步进应力加速性能退化可靠性分析方法
CN110046453A (zh) 一种激光雷达的寿命预测方法
CN113484813A (zh) 多环境应力下的智能电表故障率预估方法及系统
CN111458661A (zh) 一种配电网线变关系诊断方法、装置及系统
CN110400231A (zh) 一种加权非线性贝叶斯的电能计量设备失效率预估方法
CN111881124A (zh) 一种基于改进算法的状态估计的数据处理方法及系统
CN112816881A (zh) 电池压差异常检测方法、装置及计算机存储介质
CN112485693A (zh) 一种基于温度概率密度函数的电池健康状态快速评估方法
CN110895626B (zh) 基于留一交叉验证的性能退化模型精度验证方法
CN112329272B (zh) 综合光伏组件退化和寿命数据的剩余寿命拟合预测方法
CN114036647A (zh) 一种基于实车数据的动力电池安全风险评估方法
CN113868947A (zh) 一种智能电能表寿命预测方法及装置
CN113408192A (zh) 基于ga-fsvr的智能电表误差预测方法
CN117665620A (zh) 基于新能源汽车数据的电池健康度评估方法
CN116879763A (zh) 一种基于卡尔曼滤波算法的电池故障预警方法
Zhang et al. Reliability life prediction of VFD by constant temperature stress accelerated life tests and maximum likelihood estimation
CN106646106A (zh) 基于变点探测技术的电网故障检测方法
CN113742929B (zh) 一种针对格点气象实况的数据质量评价方法
CN113486291B (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