CN116205127A - 一种基于时间序列分解与相似性度量的装备剩余寿命方法 - Google Patents

一种基于时间序列分解与相似性度量的装备剩余寿命方法 Download PDF

Info

Publication number
CN116205127A
CN116205127A CN202211538917.8A CN202211538917A CN116205127A CN 116205127 A CN116205127 A CN 116205127A CN 202211538917 A CN202211538917 A CN 202211538917A CN 116205127 A CN116205127 A CN 116205127A
Authority
CN
China
Prior art keywords
degradation
track
similarity
data
item
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.)
Pending
Application number
CN202211538917.8A
Other languages
English (en)
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
Beijing Aerospace Automatic Control Research Institute
Original Assignee
Beihang University
Beijing Aerospace Automatic Control Research Institute
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, Beijing Aerospace Automatic Control Research Institute filed Critical Beihang University
Priority to CN202211538917.8A priority Critical patent/CN116205127A/zh
Publication of CN116205127A publication Critical patent/CN116205127A/zh
Pending legal-status Critical Current

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
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Databases & Information Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于时间序列分解与相似性度量的装备剩余寿命方法，包括如下步骤：包括离散退化轨迹数据库构建和在线退化轨迹建模与剩余寿命预测部分，离散退化轨迹数据库构建包括预处理与健康因子构建两部分，在线退化轨迹建模与剩余寿命预测包括健康度预测、退化轨迹预测、剩余寿命预测三部分，使用模糊聚类方法划分装备的退化阶段，建立健康因子，描述退化过程；采用时间序列分解模型预测健康因子的退化过程；基于相似性比较的思想，通过实例学习比较时间序列的相似性，预测剩余使用寿命。本发明提高了装备剩余寿命预测的准确率和快速性，且预测结果偏向安全，为装备健康管理提供有力保障。

Description

一种基于时间序列分解与相似性度量的装备剩余寿命方法
技术领域
本发明涉及装备的剩余寿命预测技术领域，尤其涉及一种基于时间序列分解与相似性度量的装备剩余寿命方法。
背景技术
装备测试参数众多，且在工作过程中必须承受极高的温度、压力和高频振动。在这种环境下，其各部件的功能随着时间的推移而继续恶化，剩余使用寿命也在下降。功能退化的装备很难很好地执行其任务，甚至造成严重后果。因此，长期以来，科研人员一直在研究和提出装备剩余使用寿命的预测方法。
发明内容
本发明目的是提供了一种基于时间序列分解与相似性度量的装备剩余寿命方法，以解决上述问题。
本发明解决技术问题采用如下技术方案：
一种基于时间序列分解与相似性度量的装备剩余寿命方法，包括如下步骤：
步骤S1、在离线阶段，对已有历史数据进行预处理，并进行特征选择和特征提取以降低维度；
步骤S2、通过模糊聚类对不同时刻特征提取融合得到的性能数据进行聚类，划分健康阶段、退化阶段、临近失效阶段，并得到各自的聚类中心点；
步骤S3、将聚类中心点作为评估基准，结合各时刻性能数据通过距离对健康度进行度量，设计健康因子转化模型，基于历史数据得到设备历史退化轨迹库；
步骤S4、对在线数据，经过步骤S1中的处理后，以历史数据的聚类中心进行健康阶段隶属度预测，完成健康因子转化；
步骤S5、通过时间序列分解，将退化轨迹分解为趋势项、周期项和随机项；
步骤S6、用ARIMA方法对趋势项建模，用核密度估计与拒绝接受采样完成随机项的生成，通过快速傅里叶分解得到周期项各频率的幅值、相角，选取模值较大的若干频段完成对周期项的建模；
步骤S7、分别对处理后的趋势项、周期项、随机项进行预测，再加和完成对退化轨迹的预测；
步骤S8、通过滑动平均滤波方法对退化轨迹进行平滑，通过度量轨迹窗口与历史轨迹滑动窗口的相似度比较，选取历史估计中与加入预测片段的在线轨迹最相近的若干片段，动态加权得到预测寿命。
进一步的，所述步骤S1具体包括：
对已有历史数据进行预处理，包括用最大最小归一化方法进行标准化；使用的特征选择方法是过滤式特征选择方法，以历史数据中各维特征的统计量和相关系数为评价指标，筛选出表现更优的特征子集，使用主成分分析法进行特征降维，降维后的特征是两维。
进一步的，主要筛选依据包括方差、Pearson系数和Spearman秩相关系数：方差：
Figure BDA0003976283480000021
其中，E表示随机变量X的数学期望，
Figure BDA0003976283480000031
是样本均值，n是样本容量，xi是随机变量X的n个样本值之一。
Pearson系数：
Figure BDA0003976283480000032
其中，xi和yi分别表示随机变量X和Y的n个样本值之一，
Figure BDA0003976283480000033
和
Figure BDA0003976283480000034
分别是X和Y的样本均值；
Spearman秩相关系数：设随机变量X和Y的2个随机样本为(x1,y1)，...，(xn,yn)，将x1，...，xn和y1，...，yn按升序方式进行排列，则X和Y的Spearman秩相关系数为
Figure BDA0003976283480000035
其中，R(xi)和R(yi)分别为xi和yi的秩。
进一步的，所述步骤S2具体包括：
步骤S1中得到的二维特征作为输入数据，设其为二维时间序列Z，
Figure BDA0003976283480000036
其中，N是二维时间序列的总长度，zi(i＝1,...,n)是二维时间序列数据，R2表示二维空间；
通过FCM模糊聚类算法对性能数据进行模糊聚类分析，设定聚类数目为三类，得到聚类中心矩阵C和隶属度矩阵U，
Figure BDA0003976283480000037
其中，C由3个聚类中心ci(i＝1，2，3)组成，U中uij代表了第j个数据对第i个聚类中心的隶属度。
进一步的，所述步骤S3具体包括：
根据各时刻性能数据对三阶段的隶属度，将各时刻与三个阶段对应起来，完成退化阶段划分，并以正常阶段的聚类中心作为基准记作
Figure BDA0003976283480000041
基准选取后，计算退化数据集中的每个数据与该基准的距离，每个发动机将得到一条健康状态变化曲线；在整个退化过程中，健康状态变化曲线呈上升趋势；假设发动机共有h个运行周期，则得到一条长度为h的距离变化轨迹：
Figure BDA0003976283480000042
为使健康因子满足相关定义，将得到的距离度量结果转化为由[0，1]区间内数值代表的健康度；对于第i个数据，它对基准的距离Di和它所表征的健康度之间的转化表达式如下：
Figure BDA0003976283480000043
其中Dmin表示Di中的最小值，Dmax表示Di中的最大值。
进一步的，所述步骤S6，具体步骤包括：
对于趋势项，使用ARMA(p，q)的扩展模型ARIMA(p，d，q)，该模型表示为：
Figure BDA0003976283480000044
其中p和q是模型的阶数；d是差分次数，d属于正整数；φi和θj为随机过程的系数(i＝1,2,…,p；j＝1,2,…,q)，特别地，假定φ0＝θ0＝-1；L是滞后算子；Xt是一个平稳的时间序列，εt为零均值白噪声序列。
通过快速傅里叶分解得到原始周期信号的频谱图，将能量较弱的频段滤除有效提取重要的周期信息，完成对周期项建模；
对随机项进行核密度估计得到其概率分布，完成对随机项的建模；在预测过程中，则通过拒绝接受采样从该分布中采样得到离散化的分布，生成随机项。
进一步的，所述步骤S8具体包括：
首先采用滑动平均滤波的方法平滑退化轨迹，同时采用边界对称延拓的方法解决边界平滑的问题；设有长度为m的待预测退化轨迹为L＝{l1,l2,...,lm}，以及n条历史退化轨迹{S1,...,Si,...,Sn}，其中
Figure BDA0003976283480000051
{S1,...,Si,...,Sn}长度分别为q1,…,qi,…,qn；
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L'，与各条历史轨迹上滑动窗口截取的片段进行相似性度量；对于待预测的退化轨迹根据其数据长度及预测结果的准确度，选择整体预测效果最佳的滑动窗口长度；对各历史退化轨迹，比较各位置滑动窗口截取片段的相似性，对各历史退化轨迹的最相似片段进行定位匹配；计算两段轨迹的欧氏距离来衡量其相似性；
通过对度量结果排序，根据预测效果动态筛选保留数目；同时设计三种加权策略，对寿命进行预测。
进一步的，计算两段轨迹的欧氏距离来衡量其相似性的方法如下：
对片段L'和Z'的距离度量表达式如下：
Figure BDA0003976283480000061
对于历史退化轨迹Si一共截取mi＝qi-l+1个片段，则对Si计算出mi个相似度；将计算出最小欧式距离的片段为最相似片段，设该最小距离设为di，该片段所对应的剩余寿命为ri；对退化模式库中的n条历史退化轨迹进行同样的运算，得到n个最小距离d1,d2,...,dn和对应的剩余寿命，组成集合RS＝{i:[di,ri]|i＝1,2,…,n}。
进一步的，三种加权策略包括：
策略A为朴素平均法，即将筛选得到的k个候选结果直接平均预测剩余寿命RULA：
Figure BDA0003976283480000062
其中，ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权计算剩余寿命RULB，其权重通过相似度计算得到：
Figure BDA0003976283480000063
Figure BDA0003976283480000064
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，wi表示其权重，ri为筛选出的第i个最相似片段所对应的剩余寿命；
策略C则是通过定义反比例加权函数进行加权计算得到剩余寿命RULC，其权重通过相似度计算得到：
Figure BDA0003976283480000071
Figure BDA0003976283480000072
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，wi表示其权重，ri为筛选出的第i个最相似片段所对应的剩余寿命。
有益效果：
本发明提供了一种基于时间序列分解与相似性度量的装备剩余寿命预测方法，采用时间序列分解模型预测健康因子的退化过程。基于相似性比较的思想，通过实例学习比较时间序列的相似性，准确、快速又偏向安全地预测装备剩余使用寿命。
附图说明
图1为本发明的方法流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
一种基于时间序列分解与相似性度量的装备剩余寿命方法，包括如下步骤：
步骤S1、在离线阶段，对已有历史数据进行预处理，并进行特征选择和特征提取以降低维度；
步骤S2、通过模糊聚类对不同时刻特征提取融合得到的性能数据进行聚类，划分健康阶段、退化阶段、临近失效阶段，并得到各自的聚类中心点；
步骤S3、将聚类中心点作为评估基准，结合各时刻性能数据通过距离对健康度进行度量，设计健康因子转化模型，基于历史数据得到设备历史退化轨迹库；
步骤S4、对在线数据，经过步骤S1中的处理后，以历史数据的聚类中心进行健康阶段隶属度预测，完成健康因子转化；
步骤S5、通过时间序列分解，将退化轨迹分解为趋势项、周期项和随机项；
步骤S6、用ARIMA方法对趋势项建模，用核密度估计与拒绝接受采样完成随机项的生成，通过快速傅里叶分解得到周期项各频率的幅值、相角，选取模值较大的若干频段完成对周期项的建模；
步骤S7、分别对处理后的趋势项、周期项、随机项进行预测，再加和完成对退化轨迹的预测；
步骤S8、通过滑动平均滤波方法对退化轨迹进行平滑，通过度量轨迹窗口与历史轨迹滑动窗口的相似度比较，选取历史估计中与加入预测片段的在线轨迹最相近的若干片段，动态加权得到预测寿命。
进一步的，所述步骤S1具体包括：
对已有历史数据进行预处理，包括用最大最小归一化方法进行标准化；使用的特征选择方法是过滤式特征选择方法，以历史数据中各维特征的统计量和相关系数为评价指标，筛选出表现更优的特征子集，主要筛选依据包括方差、Pearson系数和Spearman秩相关系数：
方差：
Figure BDA0003976283480000091
其中，E表示随机变量X的数学期望，
Figure BDA0003976283480000092
是样本均值，n是样本容量，xi是随机变量X的n个样本值之一，。
Pearson系数：
Figure BDA0003976283480000093
其中，xi和yi分别表示随机变量X和Y的n个样本值之一，
Figure BDA0003976283480000094
和
Figure BDA0003976283480000095
分别是X和Y的样本均值；
Spearman秩相关系数：设随机变量X和Y的2个随机样本为(x1,y1)，...，(xn,yn)，将x1，...，xn和y1，...，yn按升序方式进行排列，则X和Y的Spearman秩相关系数为
Figure BDA0003976283480000096
其中，R(xi)和R(yi)分别为xi和yi的秩。
初步筛选后的性能参数维度仍然较高，为便于进行模糊聚类，同时较好地观察聚类效果，使用主成分分析法(PCA)进行特征降维，降维后的特征是两维。
进一步的，所述步骤S2具体包括：
步骤S1中得到的二维特征作为输入数据，设其为二维时间序列Z，
Figure BDA0003976283480000101
其中，N是二维时间序列的总长度，zi(i＝1,...,n)是二维时间序列数据，R2表示二维空间。
通过FCM模糊聚类算法对性能数据进行模糊聚类分析，设定聚类数目为三类，得到聚类中心矩阵C和隶属度矩阵U，
Figure BDA0003976283480000102
其中，C由3个聚类中心ci(i＝1，2，3)组成，U中uij代表了第j个数据对第i个聚类中心的隶属度。
进一步的，所述步骤S3具体包括：
根据各时刻性能数据对三阶段的隶属度，将各时刻与三个阶段对应起来，完成退化阶段划分，并以正常阶段的聚类中心作为基准记作
Figure BDA0003976283480000104
基准选取后，计算退化数据集中的每个数据与该基准的距离，每个发动机将得到一条健康状态变化曲线；在整个退化过程中，健康状态变化曲线呈上升趋势；假设发动机共有h个运行周期，则得到一条长度为h的距离变化轨迹：
Figure BDA0003976283480000103
为使健康因子满足相关定义，将得到的距离度量结果转化为由[0，1]区间内数值代表的健康度；对于第i个数据，它对基准的距离Di和它所表征的健康度之间的转化表达式如下：
Figure BDA0003976283480000111
其中Dmin表示Di中的最小值，Dmax表示Di中的最大值。
进一步的，所述步骤S6，具体步骤包括：
对于趋势项，使用ARMA(p，q)的扩展模型ARIMA(p，d，q)，该模型表示为：
Figure BDA0003976283480000112
其中p和q是模型的阶数；d是差分次数，d属于正整数；φi和θj为随机过程的系数(i＝1,2,…,p；j＝1,2,…,q)，特别地，假定φ0＝θ0＝-1；L是滞后算子；Xt是一个平稳的时间序列，εt为零均值白噪声序列。
通过快速傅里叶分解得到原始周期信号的频谱图，将能量较弱的频段滤除有效提取重要的周期信息，完成对周期项建模；
对随机项进行核密度估计得到其概率分布，完成对随机项的建模；在预测过程中，则通过拒绝接受采样从该分布中采样得到离散化的分布，生成随机项。
进一步的，所述步骤S8具体包括：
首先采用滑动平均滤波的方法平滑退化轨迹，同时采用边界对称延拓的方法解决边界平滑的问题；设有长度为m的待预测退化轨迹为L＝{l1,l2,...,lm}，以及n条历史退化轨迹{S1,...,Si,...,Sn}，其中
Figure BDA0003976283480000122
{S1,...,Si,...,Sn}长度分别为q1,…,qi,…,qn；
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L'，与各条历史轨迹上滑动窗口截取的片段进行相似性度量；对于待预测的退化轨迹根据其数据长度及预测结果的准确度，选择整体预测效果最佳的滑动窗口长度；对各历史退化轨迹，比较各位置滑动窗口截取片段的相似性，对各历史退化轨迹的最相似片段进行定位匹配；计算两段轨迹的欧氏距离来衡量其相似性，对片段L'和Z'的距离度量表达式如下：
Figure BDA0003976283480000121
对于历史退化轨迹Si一共截取mi＝qi-l+1个片段，则对Si计算出mi个相似度；将计算出最小欧式距离的片段为最相似片段，设该最小距离设为di，该片段所对应的剩余寿命为ri；对退化模式库中的n条历史退化轨迹进行同样的运算，得到n个最小距离d1,d2,...,dn和对应的剩余寿命，组成集合RS＝{i:[di,ri]|i＝1,2,…,n}。
通过对度量结果排序，根据预测效果动态筛选保留数目；同时设计三种加权策略，对寿命进行预测。三种加权策略包括：
策略A为朴素平均法，即将筛选得到的k个候选结果直接平均预测剩余寿命：
Figure BDA0003976283480000131
其中，ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权，其权重通过相似度计算得到：
Figure BDA0003976283480000132
Figure BDA0003976283480000133
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略C则是通过定义反比例加权函数进行加权，其权重通过相似度计算得到：
Figure BDA0003976283480000134
Figure BDA0003976283480000135
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，ri为筛选出的第i个最相似片段所对应的剩余寿命。
实施例2
如图1所示，本发明提供了一种基于时间序列分解与相似性度量的剩余寿命预测方法，包括如下步骤：
步骤1、在离线阶段，对已有装备历史数据进行预处理，并进行特征选择和特征提取以降低特征的维度，特征选择方法采用主成分分析法(PCA)。装备的历史数据类型主要包括系统的电流、电压、姿态角、框架角、速度、位置、加速度等表征装备运行状态的信息。
具体地，主成分分析法的各步骤详细说明如下：
步骤1.1、对标准化后的数据计算相关系数矩阵；
步骤1.2、求出协方差矩阵的特征值和特征向量；
步骤1.3、确定主成分个数为2，选取主成分。
步骤2、通过模糊聚类对不同时刻特征提取融合得到的性能数据进行聚类，划分健康阶段、退化阶段、临近失效阶段，并得到各自的聚类中心点；
具体地，模糊聚类的各步骤详细说明如下：
步骤2.1、定义聚类对象：
设论域U＝[u1,u2,…,un]T是数为n的空间，每个ui样本有m个特征，即ui＝[xi1,xi2,…,xim](i＝1,2,…,n)，于是得到原始数据矩阵为
Figure BDA0003976283480000141
式中:xnm为第n个分类对象第m个特征的原始数据。
步骤2.2、数据标准化：
在实际问题中数据的量纲一般都不同，需要将数据进行标准化处理，从而相互之间可以比较。此处数据标准化是将数据变换到区间[0，1]，以满足模糊矩阵的要求。本发明中采用最大最小归一化方法进行标准化。计算式为：
Figure BDA0003976283480000142
式中k＝1,2,…,m,i＝1,2,…,n，xik表示第i个分类对象第k个特征的原始数据，
Figure BDA0003976283480000143
表示所有分类对象中第k个特征的原始数据中的最大值，
Figure BDA0003976283480000144
表示所有分类对象中第k个特征的原始数据中的最小值，显然有0≤xi'k≤1，即消除了量纲的影响。
步骤2.3、建立模糊相似矩阵：
由原始数据构造的模糊相似矩阵是后期聚类的重要依据，聚类的正确与否完全取决于该矩阵。本发明采用最大最小法来构造模糊相似矩阵，计算式为
Figure BDA0003976283480000151
式中：∧为“取小”运算；∨为“取大”运算，R为模糊相似矩阵，rij表示R的第i行第j列元素，xik表示第i个分类对象第k个特征的原始数据，xjk表示第j个分类对象第k个特征的原始数据。
步骤2.4、聚类：
本发明采用基于目标函数的模糊聚类分析方法，该算法实现过程中不需要人为的干预，是基于对目标函数的优化基础上的一种无监督的数据聚类方法。聚类结果是每一个数据点对聚类中心的隶属程度，该隶属程度用一个数值来表示。该算法允许同一样本属于多个不同的类。
目标函数为：
Figure BDA0003976283480000152
Figure BDA0003976283480000153
其中，C是聚类数，N是样本个数，U是隶属度矩阵，V是聚类中心，
Figure BDA0003976283480000154
表示第i个样本(即分类对象)的第k个特征所对应的隶属度矩阵；m是模糊因子用来决定聚类结果模糊度的权重指数，m∈[1,∞)，取其经验值范围为1.5≤m≤2.5；xi表示第i个样本，vk表示第k个特征的聚类中心，||xi-vk||表示取二者的欧氏距离。
步骤2.5、最佳阈值λ的确定：
在聚类中选择不同的阈值λ可得到不同的分类，λ取值越高分类数越多，反之越少.在实际问题中需要选择合适的λ确定样本的聚类数，一般情况下采用经验法和F统计量来确定最佳的阈值λ。
①经验法：由多位有经验的专家根据实际调整阈值λ，选择得到合适的分类数.
②F统计量：假设在阈值为λ时分类数为r，第j类的样本数为nj，其样本标记为
Figure BDA0003976283480000155
第j类的聚类中心为：
Figure BDA0003976283480000161
其中，
Figure BDA0003976283480000162
xik (j)表示第j类中第i个样本的第k个特征的原始数据。
则F统计量服从自由度为r-1、n-r的F分布，F统计量计算式为
Figure BDA0003976283480000163
式中,
Figure BDA0003976283480000164
分子表征不同类之间的距离；分母表征类内样本之间的距离。F统计量值越大，说明类内差异小、类间差异大，即分类效果越好。在显著水平α下，如果F>Fα(r-1，n-r)，根据统计学原理说明在该显著水平下分类是合理的，即不同类之间差异是显著的。
步骤3、将聚类中心点作为评估基准，结合各时刻性能数据通过距离对健康度进行度量，设计健康因子转化模型，基于历史数据得到设备历史退化轨迹库；
步骤4、对在线数据，经过步骤1中的处理后，以历史数据的聚类中心进行健康阶段隶属度预测，完成健康因子转化；
步骤5、通过时间序列分解，将退化轨迹分解为趋势项、周期项和随机项；
具体地，时间序列分解的各步骤详细说明如下：
对于一个给定的时间序列Y(t)，STL可以将其分解成三个相加的成分：趋势项T(t)、周期项S(t)、和随机项R(t)，即:
Y(t)＝T(t)+S(t)+R(t)
其中，周期项S(t)是指的具有周期性的成分(如，12个月，4个季度等)。STL包含内循环和外循环这两个嵌套的循环，在内循环中，主要是通过LOESS提取趋势项T(t)和周期项S(t)。在外循环中，计算随机项R(t)，并根据随机项，计算样本鲁棒权重Pr，来减少异常值对LOESS平滑结果的影响。
首先，假设内循环的循环次数为i＝1，2，…，I，初始化趋势项：T(0)(t)＝0，对于每次内循环，都有以下6个步骤:
步骤①：去趋势项。通过用原始序列Y(t)减去上个循环得到的趋势项T(i-1)(t)，去趋势的时间序列T(i)(t)：
Figure BDA0003976283480000171
步骤②：周期子序列平滑。从去趋势序列T(i)(t)中提取周期子序列，然后用LOESS(K1)每个周期子序列，同时向前后各延伸一个周期L；再将这些周期子序列按照时间先后的顺序排列，组成序列C(i)(t)，因此这里时间序列长度从T扩展到了T+2L。
步骤③：低通道滤波。对序列C(i)(t)依次做三个长度分别为L(即一个周期长度)、L、以及3的移动平均，最后再用LOESS(K2)平滑，得到序列E(i)(t),t＝1,2,...,N。与步骤②有所不同的是，这里时间序列长度从C(i)(t)的T+2L变回了E(i)(t)的T。
步骤④:计算周期项。根据下式计算周期项S(i)(t):
S(i)(t)＝C(i)(t)-E(i)(t),t＝1,2,...,T
步骤⑤:去周期项。从原始序列Y(t)中减去周期项S(i)(t)，得到去周期项的时间序列S(i)(t):
S(i)(t)＝Y(t)-S(i)(t),t＝1,2,...,T
步骤⑥:计算趋势项。使用LOESS(K3)对去周期项S(i)(t)平滑，得到趋势项T(i)(t)。
当上述6个步骤内循环完成时，可以得到趋势项T(t)和周期项S(t)。
在外循环中，首先根据下式计算随机项R(t):
R(t)＝Y(t)-T(t)-S(t)
随机项R(t)是用来反映原始数据的异常偏离情况，当t时刻随机项较大时，表示该样本异常情况较为严重，会影响内循环过程中LOESS的平滑效果。为了解决这个问题，引入样本鲁棒权重ρt，由于修正异常值的影响，ρt的计算公式如下：
Figure BDA0003976283480000181
其中，
Figure BDA0003976283480000182
median(|R(t)|)表示随机项R(t)绝对值的中位数。当求出所有样本对应的样本鲁棒权重后，在下一次的内循环中，步骤②至步骤⑥中所有的LOESS平滑过程，在寻找到邻近点后，将邻近点的数值再乘以样本鲁棒权重ρt。从式中可以看出，当样本t时刻的随机项R(t)较大时，样本获得的鲁棒权ρt就越小，当R(t)6median(|R(t)|)时，样本的鲁棒权重ρt＝0，表明LOESS平滑过程中，该点不会被考虑。
当外循环也结束时，可以得到原始序列Y(t)的三个分解成分：趋势项T(t)、周期项S(t)、和随机项R(t)。
步骤6、用ARIMA方法对趋势项建模；
具体地，ARIMA方法的各步骤详细说明如下：
步骤6.1、通过一阶差分法，对非平稳时间序列数据进行平稳操作；
步骤6.2、确定ARIMA模型的阶数，p和q的下界均为1，上界均为3，采用ACF/PACF图，找到最适合的p和q的取值；
步骤6.3、使用AICc等标准对其效果进行评估，选出最优模型；
步骤6.4、对模型进行残差分析，若模型的残差为白噪声，则建模结束，否则步骤6.2到步骤6.4。
步骤7、分别对处理后的趋势项、周期项、随机项进行预测，再完成对退化轨迹的预测；
步骤8、通过滑动平均滤波方法对退化轨迹进行平滑，通过度量轨迹窗口与历史轨迹滑动窗口的相似度比较，选取历史估计中与加入预测片段的在线轨迹最相近的若干片段，通过三种加权策略加和取平均值得到预测寿命
本发明公开提供了一种基于时间序列分解与相似性度量的剩余寿命预测方法，包括离散退化轨迹数据库构建和在线退化轨迹建模与剩余寿命预测两大模块，离线模块包括预处理与健康因子构建两部分，在线模块包括健康度预测、退化轨迹预测、剩余寿命预测三部分。基于模糊聚类方法划分装备的退化阶段，建立健康因子，描述退化过程。基于时间序列分解模型预测健康因子的退化过程。基于相似性比较的思想，通过实例学习比较时间序列的相似性，预测剩余使用寿命。计算快速，预测结果准确且偏向安全，为装备的安全运行提供了保障。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，包括如下步骤：
步骤S1、在离线阶段，对已有历史数据进行预处理，并进行特征选择和特征提取以降低维度；
步骤S2、通过模糊聚类对不同时刻特征提取融合得到的性能数据进行聚类，划分健康阶段、退化阶段、临近失效阶段，并得到各自的聚类中心点；
步骤S3、将聚类中心点作为评估基准，结合各时刻性能数据通过距离对健康度进行度量，设计健康因子转化模型，基于历史数据得到设备历史退化轨迹库；
步骤S4、对在线数据，经过步骤S1中的处理后，以历史数据的聚类中心进行健康阶段隶属度预测，完成健康因子转化；
步骤S5、通过时间序列分解，将退化轨迹分解为趋势项、周期项和随机项；
步骤S6、用ARIMA方法对趋势项建模，用核密度估计与拒绝接受采样完成随机项的生成，通过快速傅里叶分解得到周期项各频率的幅值、相角，选取模值较大的若干频段完成对周期项的建模；
步骤S7、分别对处理后的趋势项、周期项、随机项进行预测，再加和完成对退化轨迹的预测；
步骤S8、通过滑动平均滤波方法对退化轨迹进行平滑，通过度量轨迹窗口与历史轨迹滑动窗口的相似度比较，选取历史估计中与加入预测片段的在线轨迹最相近的若干片段，动态加权得到预测寿命。
2.根据权利要求1所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，所述步骤S1具体包括：
对已有历史数据进行预处理，包括用最大最小归一化方法进行标准化；使用的特征选择方法是过滤式特征选择方法，以历史数据中各维特征的统计量和相关系数为评价指标，筛选出表现更优的特征子集，使用主成分分析法进行特征降维，降维后的特征是两维。
3.根据权利要求2所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，主要筛选依据包括方差、Pearson系数和Spearman秩相关系数：
方差：
Figure FDA0003976283470000021
其中，E表示随机变量X的数学期望，X是样本均值，n是样本容量，xi是随机变量X的n个样本值之一；
Pearson系数：
Figure FDA0003976283470000022
其中，xi和yi分别表示随机变量X和Y的n个样本值之一，x和y分别是X和Y的样本均值；
Spearman秩相关系数：设随机变量X和Y的2个随机样本为(x1,y1)，...，(xn,yn)，将x1，...，xn和y1，...，yn按升序方式进行排列，则X和Y的Spearman秩相关系数为
Figure FDA0003976283470000023
其中，R(xi)和R(yi)分别为xi和yi的秩。
4.根据权利要求3所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，所述步骤S2具体包括：
步骤S1中得到的二维特征作为输入数据，设其为二维时间序列Z，
Figure FDA0003976283470000031
其中，N是二维时间序列的总长度，zi(i＝1,...,n)是二维时间序列数据，R2表示二维空间；
通过FCM模糊聚类算法对性能数据进行模糊聚类分析，设定聚类数目为三类，得到聚类中心矩阵C和隶属度矩阵U，
Figure FDA0003976283470000032
其中，C由3个聚类中心ci(i＝1，2，3)组成，U中uij代表了第j个数据对第i个聚类中心的隶属度。
5.根据权利要求4所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，所述步骤S3具体包括：
根据各时刻性能数据对三阶段的隶属度，将各时刻与三个阶段对应起来，完成退化阶段划分，并以正常阶段的聚类中心作为基准记作
Figure FDA0003976283470000033
基准选取后，计算退化数据集中的每个数据与该基准的距离，每个发动机将得到一条健康状态变化曲线；在整个退化过程中，健康状态变化曲线呈上升趋势；假设发动机共有h个运行周期，则得到一条长度为h的距离变化轨迹：
Figure FDA0003976283470000041
为使健康因子满足相关定义，将得到的距离度量结果转化为由[0，1]区间内数值代表的健康度；对于第i个数据，它对基准的距离Di和它所表征的健康度之间的转化表达式如下：
Figure FDA0003976283470000042
其中Dmin表示Di中的最小值，Dmax表示Di中的最大值。
6.根据权利要求5所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，所述步骤S6，具体步骤包括：
对于趋势项，使用ARMA(p，q)的扩展模型ARIMA(p，d，q)，该模型表示为：
Figure FDA0003976283470000043
其中p和q是模型的阶数；d是差分次数，d属于正整数；φi和θj为随机过程的系数(i＝1,2,…,p；j＝1,2,…,q)，特别地，假定φ0＝θ0＝-1；L是滞后算子；Xt是一个平稳的时间序列，εt为零均值白噪声序列；
通过快速傅里叶分解得到原始周期信号的频谱图，将能量较弱的频段滤除有效提取重要的周期信息，完成对周期项建模；
对随机项进行核密度估计得到其概率分布，完成对随机项的建模；在预测过程中，则通过拒绝接受采样从该分布中采样得到离散化的分布，生成随机项。
7.根据权利要求6所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，所述步骤S8具体包括：
首先采用滑动平均滤波的方法平滑退化轨迹，同时采用边界对称延拓的方法解决边界平滑的问题；设有长度为m的待预测退化轨迹为L＝{l1,l2,…,lm}，以及n条历史退化轨迹{S1,…,Si,…,Sn}，其中
Figure FDA0003976283470000051
其长度分别为q1,…,qi,…,qn；
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L'，与各条历史轨迹上滑动窗口截取的片段进行相似性度量；对于待预测的退化轨迹根据其数据长度及预测结果的准确度，选择整体预测效果最佳的滑动窗口长度；对各历史退化轨迹，比较各位置滑动窗口截取片段的相似性，对各历史退化轨迹的最相似片段进行定位匹配；计算两段轨迹的欧氏距离来衡量其相似性；
通过对度量结果排序，根据预测效果动态筛选保留数目；同时设计三种加权策略，对寿命进行预测。
8.根据权利要求7所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，计算两段轨迹的欧氏距离来衡量其相似性的方法如下：
对片段L'和Z'的距离度量表达式如下：
Figure FDA0003976283470000061
对于历史退化轨迹Si一共截取mi＝qi-l+1个片段，则对Si计算出mi个相似度；将计算出最小欧式距离的片段为最相似片段，设该最小距离设为di，该片段所对应的剩余寿命为ri；对退化模式库中的n条历史退化轨迹进行同样的运算，得到n个最小距离d1,d2,…,dn和对应的剩余寿命，组成集合RS＝{i:[di,ri]|i＝1,2,…,n}。
9.根据权利要求7所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法，其特征在于，三种加权策略包括：
策略A为朴素平均法，即将筛选得到的k个候选结果直接平均预测剩余寿命RULA：
Figure FDA0003976283470000062
其中，ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权计算剩余寿命RULB，其权重通过相似度计算得到：
Figure FDA0003976283470000063
Figure FDA0003976283470000064
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，wi表示其权重，ri为筛选出的第i个最相似片段所对应的剩余寿命；
策略C则是通过定义反比例加权函数进行加权计算得到剩余寿命RULC，其权重通过相似度计算得到：
Figure FDA0003976283470000071
Figure FDA0003976283470000072
其中，di表示筛选出的k个候选结果中第i个最小欧氏距离，wi表示其权重，ri为筛选出的第i个最相似片段所对应的剩余寿命。
CN202211538917.8A 2022-12-01 2022-12-01 一种基于时间序列分解与相似性度量的装备剩余寿命方法 Pending CN116205127A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211538917.8A CN116205127A (zh) 2022-12-01 2022-12-01 一种基于时间序列分解与相似性度量的装备剩余寿命方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211538917.8A CN116205127A (zh) 2022-12-01 2022-12-01 一种基于时间序列分解与相似性度量的装备剩余寿命方法

Publications (1)

Publication Number Publication Date
CN116205127A true CN116205127A (zh) 2023-06-02

Family

ID=86511923

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211538917.8A Pending CN116205127A (zh) 2022-12-01 2022-12-01 一种基于时间序列分解与相似性度量的装备剩余寿命方法

Country Status (1)

Country Link
CN (1) CN116205127A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116563055A (zh) * 2023-07-07 2023-08-08 中国科学院地理科学与资源研究所 一种基于多源数据融合的风能潜力评估方法
CN117056692A (zh) * 2023-10-09 2023-11-14 山东芯赛思电子科技有限公司 一种SiC基电机驱动器件的老化预测方法
CN117310118A (zh) * 2023-11-28 2023-12-29 济南中安数码科技有限公司 一种地下水污染可视化监测方法
CN117349602A (zh) * 2023-12-06 2024-01-05 江西省水投江河信息技术有限公司 一种水利设施运行状态预测方法、系统及计算机

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116563055A (zh) * 2023-07-07 2023-08-08 中国科学院地理科学与资源研究所 一种基于多源数据融合的风能潜力评估方法
CN116563055B (zh) * 2023-07-07 2023-11-03 中国科学院地理科学与资源研究所 一种基于多源数据融合的风能潜力评估方法
CN117056692A (zh) * 2023-10-09 2023-11-14 山东芯赛思电子科技有限公司 一种SiC基电机驱动器件的老化预测方法
CN117056692B (zh) * 2023-10-09 2024-01-23 山东芯赛思电子科技有限公司 一种SiC基电机驱动器件的老化预测方法
CN117310118A (zh) * 2023-11-28 2023-12-29 济南中安数码科技有限公司 一种地下水污染可视化监测方法
CN117310118B (zh) * 2023-11-28 2024-03-08 济南中安数码科技有限公司 一种地下水污染可视化监测方法
CN117349602A (zh) * 2023-12-06 2024-01-05 江西省水投江河信息技术有限公司 一种水利设施运行状态预测方法、系统及计算机

Similar Documents

Publication Publication Date Title
CN116205127A (zh) 一种基于时间序列分解与相似性度量的装备剩余寿命方法
CN109472110B (zh) 一种基于lstm网络和arima模型的航空发动机剩余使用寿命预测方法
CN110610035B (zh) 一种基于gru神经网络的滚动轴承剩余寿命预测方法
CN108536971B (zh) 一种基于贝叶斯模型的结构损伤识别方法
CN111222290B (zh) 一种基于多参数特征融合的大型设备剩余使用寿命预测方法
Wang et al. Remaining useful life estimation using functional data analysis
CN105975443B (zh) 基于Lasso的网络异常行为检测方法及系统
CN107144428B (zh) 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法
CN103824137B (zh) 一种复杂机械设备多工况故障预测方法
CN109376892B (zh) 一种基于设备所处生命周期阶段的设备状态预测方法
Mo et al. Multi-head CNN-LSTM with prediction error analysis for remaining useful life prediction
CN109767043B (zh) 一种电力负荷时间序列大数据智能建模与预测方法
CN116380445B (zh) 基于振动波形的设备状态诊断方法及相关装置
CN112785092B (zh) 一种基于自适应深层特征提取的道岔剩余寿命预测方法
CN113012766B (zh) 一种基于在线选择性集成的自适应软测量建模方法
CN113554148A (zh) 一种基于贝叶斯优化的BiLSTM电压偏差预测方法
CN113539382B (zh) 一种亚磷酸二甲酯关键工艺参数的预警定位方法及系统
Chen et al. A deep learning feature fusion based health index construction method for prognostics using multiobjective optimization
CN114417699A (zh) 泵阀故障检测方法
Zhu et al. Res-HSA: Residual hybrid network with self-attention mechanism for RUL prediction of rotating machinery
CN109947076A (zh) 一种基于贝叶斯信息准则的工业过程故障诊断方法
CN109298633A (zh) 基于自适应分块非负矩阵分解的化工生产过程故障监测方法
US20240085274A1 (en) Hybrid bearing fault prognosis with fault detection and multiple model fusion
CN113052271A (zh) 基于深度神经网络的生物发酵数据预测方法
Xiao et al. A noise-boosted remaining useful life prediction method for rotating machines under different conditions

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