CN116205127A - 一种基于时间序列分解与相似性度量的装备剩余寿命方法 - Google Patents
一种基于时间序列分解与相似性度量的装备剩余寿命方法 Download PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 36
- 238000005259 measurement Methods 0.000 title claims abstract description 26
- 238000006731 degradation reaction Methods 0.000 claims abstract description 67
- 230000015556 catabolic process Effects 0.000 claims abstract description 62
- 230000036541 health Effects 0.000 claims abstract description 41
- 239000012634 fragment Substances 0.000 claims description 18
- 239000011159 matrix material Substances 0.000 claims description 17
- 230000000737 periodic effect Effects 0.000 claims description 17
- 238000001914 filtration Methods 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000009499 grossing Methods 0.000 claims description 13
- 230000000694 effects Effects 0.000 claims description 12
- 238000006243 chemical reaction Methods 0.000 claims description 11
- 238000012216 screening Methods 0.000 claims description 10
- 238000003646 Spearman's rank correlation coefficient Methods 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 9
- YHXISWVBGDMDLQ-UHFFFAOYSA-N moclobemide Chemical compound C1=CC(Cl)=CC=C1C(=O)NCCN1CCOCC1 YHXISWVBGDMDLQ-UHFFFAOYSA-N 0.000 claims description 9
- 238000000605 extraction Methods 0.000 claims description 8
- 238000011156 evaluation Methods 0.000 claims description 7
- 230000006870 function Effects 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 7
- 238000010187 selection method Methods 0.000 claims description 7
- 230000001174 ascending effect Effects 0.000 claims description 6
- 230000009467 reduction Effects 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 5
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 238000010606 normalization Methods 0.000 claims description 5
- 230000004927 fusion Effects 0.000 claims description 4
- 241001123248 Arma Species 0.000 claims description 3
- 101100006960 Caenorhabditis elegans let-2 gene Proteins 0.000 claims description 3
- 238000007621 cluster analysis Methods 0.000 claims description 3
- 230000007850 degeneration Effects 0.000 claims description 3
- 238000012847 principal component analysis method Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 abstract description 5
- 238000000513 principal component analysis Methods 0.000 description 4
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000004836 empirical method Methods 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
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具体包括:
对已有历史数据进行预处理,包括用最大最小归一化方法进行标准化;使用的特征选择方法是过滤式特征选择方法,以历史数据中各维特征的统计量和相关系数为评价指标,筛选出表现更优的特征子集,使用主成分分析法进行特征降维,降维后的特征是两维。
Spearman秩相关系数:设随机变量X和Y的2个随机样本为(x1,y1),...,(xn,yn),将x1,...,xn和y1,...,yn按升序方式进行排列,则X和Y的Spearman秩相关系数为其中,R(xi)和R(yi)分别为xi和yi的秩。
进一步的,所述步骤S2具体包括:
步骤S1中得到的二维特征作为输入数据,设其为二维时间序列Z,
其中,N是二维时间序列的总长度,zi(i=1,...,n)是二维时间序列数据,R2表示二维空间;
通过FCM模糊聚类算法对性能数据进行模糊聚类分析,设定聚类数目为三类,得到聚类中心矩阵C和隶属度矩阵U,
其中,C由3个聚类中心ci(i=1,2,3)组成,U中uij代表了第j个数据对第i个聚类中心的隶属度。
进一步的,所述步骤S3具体包括:
基准选取后,计算退化数据集中的每个数据与该基准的距离,每个发动机将得到一条健康状态变化曲线;在整个退化过程中,健康状态变化曲线呈上升趋势;假设发动机共有h个运行周期,则得到一条长度为h的距离变化轨迹:
为使健康因子满足相关定义,将得到的距离度量结果转化为由[0,1]区间内数值代表的健康度;对于第i个数据,它对基准的距离Di和它所表征的健康度之间的转化表达式如下:
其中Dmin表示Di中的最小值,Dmax表示Di中的最大值。
进一步的,所述步骤S6,具体步骤包括:
对于趋势项,使用ARMA(p,q)的扩展模型ARIMA(p,d,q),该模型表示为:
其中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},其中{S1,...,Si,...,Sn}长度分别为q1,…,qi,…,qn;
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L',与各条历史轨迹上滑动窗口截取的片段进行相似性度量;对于待预测的退化轨迹根据其数据长度及预测结果的准确度,选择整体预测效果最佳的滑动窗口长度;对各历史退化轨迹,比较各位置滑动窗口截取片段的相似性,对各历史退化轨迹的最相似片段进行定位匹配;计算两段轨迹的欧氏距离来衡量其相似性;
通过对度量结果排序,根据预测效果动态筛选保留数目;同时设计三种加权策略,对寿命进行预测。
进一步的,计算两段轨迹的欧氏距离来衡量其相似性的方法如下:
对片段L'和Z'的距离度量表达式如下:
对于历史退化轨迹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:
其中,ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权计算剩余寿命RULB,其权重通过相似度计算得到:
其中,di表示筛选出的k个候选结果中第i个最小欧氏距离,wi表示其权重,ri为筛选出的第i个最相似片段所对应的剩余寿命;
策略C则是通过定义反比例加权函数进行加权计算得到剩余寿命RULC,其权重通过相似度计算得到:
其中,di表示筛选出的k个候选结果中第i个最小欧氏距离,wi表示其权重,ri为筛选出的第i个最相似片段所对应的剩余寿命。
有益效果:
本发明提供了一种基于时间序列分解与相似性度量的装备剩余寿命预测方法,采用时间序列分解模型预测健康因子的退化过程。基于相似性比较的思想,通过实例学习比较时间序列的相似性,准确、快速又偏向安全地预测装备剩余使用寿命。
附图说明
图1为本发明的方法流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例1
一种基于时间序列分解与相似性度量的装备剩余寿命方法,包括如下步骤:
步骤S1、在离线阶段,对已有历史数据进行预处理,并进行特征选择和特征提取以降低维度;
步骤S2、通过模糊聚类对不同时刻特征提取融合得到的性能数据进行聚类,划分健康阶段、退化阶段、临近失效阶段,并得到各自的聚类中心点;
步骤S3、将聚类中心点作为评估基准,结合各时刻性能数据通过距离对健康度进行度量,设计健康因子转化模型,基于历史数据得到设备历史退化轨迹库;
步骤S4、对在线数据,经过步骤S1中的处理后,以历史数据的聚类中心进行健康阶段隶属度预测,完成健康因子转化;
步骤S5、通过时间序列分解,将退化轨迹分解为趋势项、周期项和随机项;
步骤S6、用ARIMA方法对趋势项建模,用核密度估计与拒绝接受采样完成随机项的生成,通过快速傅里叶分解得到周期项各频率的幅值、相角,选取模值较大的若干频段完成对周期项的建模;
步骤S7、分别对处理后的趋势项、周期项、随机项进行预测,再加和完成对退化轨迹的预测;
步骤S8、通过滑动平均滤波方法对退化轨迹进行平滑,通过度量轨迹窗口与历史轨迹滑动窗口的相似度比较,选取历史估计中与加入预测片段的在线轨迹最相近的若干片段,动态加权得到预测寿命。
进一步的,所述步骤S1具体包括:
对已有历史数据进行预处理,包括用最大最小归一化方法进行标准化;使用的特征选择方法是过滤式特征选择方法,以历史数据中各维特征的统计量和相关系数为评价指标,筛选出表现更优的特征子集,主要筛选依据包括方差、Pearson系数和Spearman秩相关系数:
Spearman秩相关系数:设随机变量X和Y的2个随机样本为(x1,y1),...,(xn,yn),将x1,...,xn和y1,...,yn按升序方式进行排列,则X和Y的Spearman秩相关系数为其中,R(xi)和R(yi)分别为xi和yi的秩。
初步筛选后的性能参数维度仍然较高,为便于进行模糊聚类,同时较好地观察聚类效果,使用主成分分析法(PCA)进行特征降维,降维后的特征是两维。
进一步的,所述步骤S2具体包括:
步骤S1中得到的二维特征作为输入数据,设其为二维时间序列Z,
其中,N是二维时间序列的总长度,zi(i=1,...,n)是二维时间序列数据,R2表示二维空间。
通过FCM模糊聚类算法对性能数据进行模糊聚类分析,设定聚类数目为三类,得到聚类中心矩阵C和隶属度矩阵U,
其中,C由3个聚类中心ci(i=1,2,3)组成,U中uij代表了第j个数据对第i个聚类中心的隶属度。
进一步的,所述步骤S3具体包括:
基准选取后,计算退化数据集中的每个数据与该基准的距离,每个发动机将得到一条健康状态变化曲线;在整个退化过程中,健康状态变化曲线呈上升趋势;假设发动机共有h个运行周期,则得到一条长度为h的距离变化轨迹:
为使健康因子满足相关定义,将得到的距离度量结果转化为由[0,1]区间内数值代表的健康度;对于第i个数据,它对基准的距离Di和它所表征的健康度之间的转化表达式如下:
其中Dmin表示Di中的最小值,Dmax表示Di中的最大值。
进一步的,所述步骤S6,具体步骤包括:
对于趋势项,使用ARMA(p,q)的扩展模型ARIMA(p,d,q),该模型表示为:
其中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},其中{S1,...,Si,...,Sn}长度分别为q1,…,qi,…,qn;
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L',与各条历史轨迹上滑动窗口截取的片段进行相似性度量;对于待预测的退化轨迹根据其数据长度及预测结果的准确度,选择整体预测效果最佳的滑动窗口长度;对各历史退化轨迹,比较各位置滑动窗口截取片段的相似性,对各历史退化轨迹的最相似片段进行定位匹配;计算两段轨迹的欧氏距离来衡量其相似性,对片段L'和Z'的距离度量表达式如下:
对于历史退化轨迹Si一共截取mi=qi-l+1个片段,则对Si计算出mi个相似度;将计算出最小欧式距离的片段为最相似片段,设该最小距离设为di,该片段所对应的剩余寿命为ri;对退化模式库中的n条历史退化轨迹进行同样的运算,得到n个最小距离d1,d2,...,dn和对应的剩余寿命,组成集合RS={i:[di,ri]|i=1,2,…,n}。
通过对度量结果排序,根据预测效果动态筛选保留数目;同时设计三种加权策略,对寿命进行预测。三种加权策略包括:
策略A为朴素平均法,即将筛选得到的k个候选结果直接平均预测剩余寿命:
其中,ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权,其权重通过相似度计算得到:
其中,di表示筛选出的k个候选结果中第i个最小欧氏距离,ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略C则是通过定义反比例加权函数进行加权,其权重通过相似度计算得到:
其中,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),于是得到原始数据矩阵为
式中:xnm为第n个分类对象第m个特征的原始数据。
步骤2.2、数据标准化:
在实际问题中数据的量纲一般都不同,需要将数据进行标准化处理,从而相互之间可以比较。此处数据标准化是将数据变换到区间[0,1],以满足模糊矩阵的要求。本发明中采用最大最小归一化方法进行标准化。计算式为:
式中k=1,2,…,m,i=1,2,…,n,xik表示第i个分类对象第k个特征的原始数据,表示所有分类对象中第k个特征的原始数据中的最大值,表示所有分类对象中第k个特征的原始数据中的最小值,显然有0≤xi'k≤1,即消除了量纲的影响。
步骤2.3、建立模糊相似矩阵:
由原始数据构造的模糊相似矩阵是后期聚类的重要依据,聚类的正确与否完全取决于该矩阵。本发明采用最大最小法来构造模糊相似矩阵,计算式为
式中:∧为“取小”运算;∨为“取大”运算,R为模糊相似矩阵,rij表示R的第i行第j列元素,xik表示第i个分类对象第k个特征的原始数据,xjk表示第j个分类对象第k个特征的原始数据。
步骤2.4、聚类:
本发明采用基于目标函数的模糊聚类分析方法,该算法实现过程中不需要人为的干预,是基于对目标函数的优化基础上的一种无监督的数据聚类方法。聚类结果是每一个数据点对聚类中心的隶属程度,该隶属程度用一个数值来表示。该算法允许同一样本属于多个不同的类。
目标函数为:
其中,C是聚类数,N是样本个数,U是隶属度矩阵,V是聚类中心,表示第i个样本(即分类对象)的第k个特征所对应的隶属度矩阵;m是模糊因子用来决定聚类结果模糊度的权重指数,m∈[1,∞),取其经验值范围为1.5≤m≤2.5;xi表示第i个样本,vk表示第k个特征的聚类中心,||xi-vk||表示取二者的欧氏距离。
步骤2.5、最佳阈值λ的确定:
在聚类中选择不同的阈值λ可得到不同的分类,λ取值越高分类数越多,反之越少.在实际问题中需要选择合适的λ确定样本的聚类数,一般情况下采用经验法和F统计量来确定最佳的阈值λ。
①经验法:由多位有经验的专家根据实际调整阈值λ,选择得到合适的分类数.
则F统计量服从自由度为r-1、n-r的F分布,F统计量计算式为
式中,分子表征不同类之间的距离;分母表征类内样本之间的距离。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):
步骤②:周期子序列平滑。从去趋势序列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的计算公式如下:
其中,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秩相关系数:
5.根据权利要求4所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法,其特征在于,所述步骤S3具体包括:
基准选取后,计算退化数据集中的每个数据与该基准的距离,每个发动机将得到一条健康状态变化曲线;在整个退化过程中,健康状态变化曲线呈上升趋势;假设发动机共有h个运行周期,则得到一条长度为h的距离变化轨迹:
为使健康因子满足相关定义,将得到的距离度量结果转化为由[0,1]区间内数值代表的健康度;对于第i个数据,它对基准的距离Di和它所表征的健康度之间的转化表达式如下:
其中Dmin表示Di中的最小值,Dmax表示Di中的最大值。
6.根据权利要求5所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法,其特征在于,所述步骤S6,具体步骤包括:
对于趋势项,使用ARMA(p,q)的扩展模型ARIMA(p,d,q),该模型表示为:
其中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},其中其长度分别为q1,…,qi,…,qn;
选择长度为l的待预测退化轨迹末尾片段和时间序列分解建模得到的预测片段结合形成L',与各条历史轨迹上滑动窗口截取的片段进行相似性度量;对于待预测的退化轨迹根据其数据长度及预测结果的准确度,选择整体预测效果最佳的滑动窗口长度;对各历史退化轨迹,比较各位置滑动窗口截取片段的相似性,对各历史退化轨迹的最相似片段进行定位匹配;计算两段轨迹的欧氏距离来衡量其相似性;
通过对度量结果排序,根据预测效果动态筛选保留数目;同时设计三种加权策略,对寿命进行预测。
9.根据权利要求7所述的一种基于时间序列分解与相似性度量的装备剩余寿命方法,其特征在于,三种加权策略包括:
策略A为朴素平均法,即将筛选得到的k个候选结果直接平均预测剩余寿命RULA:
其中,ri为筛选出的第i个最相似片段所对应的剩余寿命。
策略B则是通过softmax进行加权计算剩余寿命RULB,其权重通过相似度计算得到:
其中,di表示筛选出的k个候选结果中第i个最小欧氏距离,wi表示其权重,ri为筛选出的第i个最相似片段所对应的剩余寿命;
策略C则是通过定义反比例加权函数进行加权计算得到剩余寿命RULC,其权重通过相似度计算得到:
其中,di表示筛选出的k个候选结果中第i个最小欧氏距离,wi表示其权重,ri为筛选出的第i个最相似片段所对应的剩余寿命。
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)
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 | 江西省水投江河信息技术有限公司 | 一种水利设施运行状态预测方法、系统及计算机 |
-
2022
- 2022-12-01 CN CN202211538917.8A patent/CN116205127A/zh active Pending
Cited By (7)
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 |