CN112101227A - 一种基于celmdan和sskfda的机械状态监测方法 - Google Patents
一种基于celmdan和sskfda的机械状态监测方法 Download PDFInfo
- Publication number
- CN112101227A CN112101227A CN202010972737.5A CN202010972737A CN112101227A CN 112101227 A CN112101227 A CN 112101227A CN 202010972737 A CN202010972737 A CN 202010972737A CN 112101227 A CN112101227 A CN 112101227A
- Authority
- CN
- China
- Prior art keywords
- samples
- matrix
- sample
- class
- representing
- 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.)
- Withdrawn
Links
- 238000000034 method Methods 0.000 title claims abstract description 113
- 238000012544 monitoring process Methods 0.000 title claims abstract description 50
- 230000006870 function Effects 0.000 claims abstract description 76
- 230000000737 periodic effect Effects 0.000 claims abstract description 16
- 230000009467 reduction Effects 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims description 102
- 238000004422 calculation algorithm Methods 0.000 claims description 53
- 239000013598 vector Substances 0.000 claims description 50
- 238000013507 mapping Methods 0.000 claims description 27
- 238000005457 optimization Methods 0.000 claims description 24
- 238000000354 decomposition reaction Methods 0.000 claims description 19
- 238000009826 distribution Methods 0.000 claims description 19
- 238000004458 analytical method Methods 0.000 claims description 18
- 230000001174 ascending effect Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 10
- 238000000926 separation method Methods 0.000 claims description 9
- 230000010355 oscillation Effects 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 5
- 238000011478 gradient descent method Methods 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 2
- 238000005314 correlation function Methods 0.000 claims 1
- 230000006911 nucleation Effects 0.000 claims 1
- 238000010899 nucleation Methods 0.000 claims 1
- 238000010187 selection method Methods 0.000 abstract description 3
- 230000002159 abnormal effect Effects 0.000 abstract 1
- 238000003745 diagnosis Methods 0.000 description 10
- 238000012360 testing method Methods 0.000 description 7
- 230000003044 adaptive effect Effects 0.000 description 6
- 230000008901 benefit Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 230000010354 integration Effects 0.000 description 6
- 238000012549 training Methods 0.000 description 6
- 238000000605 extraction Methods 0.000 description 5
- 238000010801 machine learning Methods 0.000 description 5
- 238000000513 principal component analysis Methods 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 239000000654 additive Substances 0.000 description 4
- 230000000996 additive effect Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000001514 detection method Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000005311 autocorrelation function Methods 0.000 description 2
- 238000003909 pattern recognition Methods 0.000 description 2
- 230000035939 shock Effects 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 238000003889 chemical engineering Methods 0.000 description 1
- 239000003638 chemical reducing agent Substances 0.000 description 1
- 239000012141 concentrate Substances 0.000 description 1
- 238000013527 convolutional neural network Methods 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 238000002790 cross-validation Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000012850 discrimination method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000011326 mechanical measurement Methods 0.000 description 1
- 238000005272 metallurgy Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000012567 pattern recognition method Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Artificial Intelligence (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于CELMDAN和SSKFDA的机械状态监测方法,采用以下方法:(1)利用CELMDAN方法把复杂振动信号分解为多个有物理意义的乘积函数。(2)提出周期调制强度PMI作为选择PFs准则的方法,能够准确选择有效的PFs。(3)提出SSKFDA维数约简方法,充分利用标签样本和未标签样本集几何信息,融合核方法、稀疏表示、流形学习以及FDA方法,更好地揭示嵌入在高维稀疏空间的低维子空间数据集,解决高维、稀疏和非线性数据维数约简问题。(4)提出一种快速SSKFDA模型选择方法。该方法通过基于类内局部结构测度最小、全局部结构测度最大准则求取最优模型参数。(5)提出全局监测统计量与Bayesian后验推理的机械已知与未知状态检测方法,较好解决大多数机械监测系统无法检测未知异常状态问题。
Description
技术领域
本发明涉及机械设备监测技术领域,特别是一种基于噪声自适应完全集成局部均值分解CELMDAN和半监督系数核Fisher鉴别分析SSKFDA的机械状态监测方法。
背景技术
旋转机械是飞机发动机,燃气轮机,减速器等设备的重要部分,广泛应用于电力、冶金、化工和机械等行业,其状态不仅导致局部破坏,影响产品质量和企业的经济效益,严重地威胁操作人员的人身安全。因此有必要深入研究机械运行状态和故障诊断,从而降低维护费用,保证生产安全稳定。振动信号含有丰富的机械设备部件状态信息,基于振动信号的旋转机械状态监测和故障诊断方法已经被证明是一种有效、主流方法。然而,机械工作环境恶劣、生产负荷和旋转速度时变性、设备零部件耦合以及故障冲击,振动信号往往表现出很强非平稳、非线性、信号调制和强背景噪声,给振动信号的故障信息抽取和故障诊断带来了严重的挑战。
为了准确对机械故障进行识别和诊断,需要从振动信号抽取与故障敏感的有用信息并进行故障识别。近年来,随着机器学习、信号处理以及人工智能的快速发展,融合先进信号处理算法和数据驱动的机器学习方法已经成为机械故障诊断的流行方法。时频分析方法已经成常用、有效的振动信号处理分析工具,比如短时傅里叶变换(short time Fouriertransform,STFT),稀疏表示、WVD(Wigner-Ville distribution,WVD)、小波变换(wavelettransformation,WT)、经验模式分解(empirical mode decomposition,EMD)、局部均值分解(local mean decomposition,LMD)等。EMD是一种新的自适应多尺度时频分析方法,但其性能往往受到模式混淆和端点效应的影响。旋转机械(比如转子,轴承或齿轮)产生故障时通常会产生高频冲击信号以及脉冲,振动信号被这些冲击信号所调制,导致故障特征被掩盖。由于旋转机械结构特性,测量故障信号往往具有典型的幅值和频率调制特性。LMD具有很强的信号解调分析能力,能够把复杂信号分解多个称为乘积信号(productionfunction,PF)的相对简单信号。与EMD相比,LMD无需使用Hilbert变换即可得到瞬态幅值和频率,不使用插值策略拟合上下包络。因此LMD不仅克服了EMD本身存在的问题,而且有用信息更加集中在少数分解成分中,非常适合非线性、非平稳信号分析,能够提高旋转机械信号故障诊断的性能。然而,LMD仍然存在受到端点效应影响导致分解结果失真而无法解释的问题。为此,人们提出一些改进算法,比如,为消除LMD端点效应问题,W.Y.Liu(2015)提出一种对原始信号左右两端进行自然扩展并保持信号的趋势不变的LMD积分扩展(Integralextension of LMD,IELM)算法;类似于EEMD思想,Y.Yang et al(2012)提出集成LMD(ensemble LMD,ELMD)算法。ELMD方法首先把多个有限幅值白噪声加到原始信号中,然后使用LMD分别把这些加噪信号分解为PFs集,并对每次实验PFs进行集成平均得到最终PFs。需要指出的是,对含噪原始振动信号,由ELMD得到含有感兴趣成分的PFs仍会存在噪声。为了消除或抑制感兴趣PFs噪声突出机械故障特征,Wang Lei(2018)使用kurtosis指标选取由ELMD求取表征故障特征信息的PFS,并使用快速kurtgram确定最优带宽对选取PFs滤波抽取冲击信号;Q.W.Gao(2018)利用多尺度熵(Multiscale entropy,ME)处理非线性信号以及抗噪(anti-noise)和抗干扰(anti-interference)能力,使用ME抽取不同尺度PFs的特征参数,准确地计算故障的复杂性。需要指出的是,ELMD仍然存在如下问题:(1)ELMD把相同尺度的信号成分以及加性噪声部分均被分解到相应的PFs,需要多次集成试验消除加性噪声部分,增加了算法的计算量;(2)施加噪声幅值对ELMD性能影响很大。施加噪声幅值太小不足以影响极点的分布,因而无法解决模式混淆问题。而噪声幅值太大导致模式包含相当大的残差,需要更多的集成试验次数;(3)不同高斯噪声加入原始信号,不同数量试验会产生不同的PFs,而不同PFs的物理意义也会有所区别,不利于系统运行状态的识别。为克服上述问题,Lei Wang(2018)提出噪声自适应完全ELMD(complete ELMD with adaptive noise,CELMDAN)算法。在每次试验中,该方法在每个分解阶段均加入一个特殊的自适应噪声,每次分离PF后会得到唯一的残差。然后把残差作为下一次分解的输入。CELMDAN较好地消除了模式混淆问题,求取的PFs具有明显的物理意义。
从振动信号中抽取包含故障信息特征是提高在线机械故障识别率的关键。一般说来,不同故障在信号的不同域中表现方式和程度不一致。因此,人们从测量信号时域和频域抽取机械故障特征参数,比如从原始信号时域抽取熵特征、Teager-Kaiser能量算子(Teager-Kaiser energy operator,TKEO),kurtosis,skewness,shape、peak等统计特征参数;从原始信号频域中抽取谱kurtosis,包络谱(envelope spectrum),谱分析、功率谱,频率均值,频率中心等。除了时域和频域特征外,多尺度时频特征,比如,小波能量特征、基于改进LMD的多尺度熵、多尺度高阶奇异谱熵、奇异值特征也成为获取原始信号的状态信息重要方法。机械多故障识别可以归结为模式识别问题,即把不同故障模态划分到多个不同子空间或者聚类。
高维数据往往包含大量的冗余信息,后续的故障识别算法往往面临“维数灾难”的问题。因此,人们往往使用维数约简算法抽取故障敏感特征,更好揭示机械状态的数据几何结构,进而提高故障诊断的精度和机械状态的理解性。常用的维数约简算法主要有主元分析(principal component analysis,PCA)、独立主元分析(Independent componentanalysis,ICA)、Fisher鉴别分析(Fisher discriminant analysis,FDA)、特征选择、流形学习算法等。流形学习方法能够利用样本数据集信息揭示数据类别隐藏几何结构,成为维数约简的研究热点,在故障诊断中得到极大的重视。Jianbo Yu.Hidden Markov modelscombining local and global information for nonlinear and multimodal processmonitoring.Journal of Process Control20:344-359,2010当中,提出利用数据集全局和局部几何结构鉴别信息的机械故障识别方法,更好地揭示高维数据中隐藏的低维故障几何信息;考虑到传统流形学习算法受到数据噪声影响以及近邻样本选择问题,Gang Chen,Fenglin Liu,Wei Huang.Sparse discriminant manifold projections for bearingfault diagnosis.Journal of Sound and vibration,399(2017)330-344中,提出一种基于稀疏保持投影和稀疏流形局部嵌入的稀疏鉴别流形投影算法,该方法通过稀疏表示自适应地选取近邻数据样本点,不仅较好地克服了传统流形学习算法中近邻点数据可能来自多个流形的问题、数据噪声影响,而且利用样本数据类别信息从高维数据抽取故障类别鉴别信息,算法还克服了out-of–sample问题。LDA以及其改进算法如直接LDA(direct LDA,DLDA)、正则化LDA(regularized LDA,RLDA)等被证明是一类有效的模式识别方法,广泛应用于人脸识别、基因识别、故障诊断等领域中。LDA通过最大化投影空间中类间散度与类内散度比例实现高维数据的维数约简。稀疏LDA方法把稀疏表示技术应用于LDA,优于传统依赖所有输入变量的LDA方法,成为重要的稀疏子空间学习方法。针对稀疏子空间学习方法没有有效利用数据鉴别信息问题,受到LDA思想启发,Lei Yu(2017)提出一种称为稀疏多最大散度差(sparse multiple maximum scatter difference,SMMSD)的稀疏子空间学习方法。该方法通过多个最大散度差准则使得低维约简空间上类内稀疏重构残差最小以及类间稀疏重构最大化,能够获得更多潜在分类鉴别信息而且避免了奇异值问题。SMMSD方法是一种新的监督稀疏子空间学习算法,其性能优于传统的维数约简算法。
发明内容
本发明所要解决的技术问题是克服现有技术的不足而提供一种基于CELMDAN和SSKFDA的机械状态监测方法,提高了非线性、不确定多机械状态识别的精度。
本发明为解决上述技术问题采用以下技术方案:
根据本发明提出的一种基于CELMDAN和SSKFDA的机械状态监测方法,包括离线建模阶段和在线监测机械状态阶段;具体步骤如下:
离线建模阶段:
初始化阶段:从旋转机械采集故障和正常状态的振动信号,类内、类间散度矩阵的正则化因子β,核函数参数初始值σ,半监督系数核Fisher鉴别分析SSKFDA算法约简后数据特征维数r,置信度水平γ',每个状态的先验概率;正常状态指非故障状态;
第1步、对振动信号归一化为[0,1]内,根据局部均值分解LMD算法以及噪声自适应完全集成局部均值分解CELMDAN算法对归一化后的振动信号进行分解,分为多个瞬时频率具有物理意义的乘积函数PF成分;
第2步、计算每个PF成分的周期调制强度PMI值,根据PMI值选取PF成分;
第3步、从第2步得到的所有PF成分中抽取PF奇异熵特征、PF时频熵特征、PF能量熵特征,分别从每一个选取的PF成分中抽取峰度Kurtosis、偏度Skewness、标准差和脉冲因子统计特征;根据PF奇异熵特征、PF时频熵特征、PF能量熵特征、抽取峰度Kurtosis、偏度Skewness、标准差和脉冲因子统计特征,生成标记样本数据集和未标记样本数据集;
第4步、基于标记样本数据集,首先根据特征空间中类内数据集尽可能靠近、类间数据集尽可能远离准则,求取最优墨西哥Mexico小波核函数参数;
第5步、根据所有样本数据集以及第4步求取的最优墨西哥Mexico小波核函数参数计算核矩阵K,根据核矩阵K计算中心化核矩阵利用扩展最优类内样本重构权重向量流形内相似性矩阵Gw∈RN×N计算类内散度Sw,利用流形间分离矩阵Gb∈RN×N、扩展最优类间样本重构权重向量计算类间散度Sb,利用总体相似矩阵Gt∈RN×N、标记和未标记样本对样本重构权重向量pi以及构建总体散度St;RN×1表示N×1的实数列向量,RN×N表示N×N的实数矩阵,N为标记样本和未标记样本的总数量;
第6步、基于第5步计算结果,构建正则化类内散度Srw和类间散度Srb,其中,Srw度量同类样本之间散度,Srb度量不同异类样本之间散度;
第7步、基于第5步和第6步计算结果,根据Srw最小、Srb最大准则,计算所示特征值和特征向量,选取前r个最大特征值对应的特征向量作为投影矩阵W=[q1,q2,…,qr];β∈(0,1)为正则化因子,q为特征向量,λ为特征值,I为N×N的单位矩阵,qp为第p个最大特征值对应的特征向量,p=1,2,…,r,r为选取特征向量的数量,Srb'表示特征空间不同类别数据之间的流形几何结构,Spt'表示特征空间全体数据之间的流形几何结构,Ssw'表示特征空间同类数据之间的流形几何结构;
第8步、保存投影矩阵、Mexico小波核函数参数、核矩阵作为求取的SSKFDA模型;
在线监测阶段:
步骤A、从运行的旋转机械采集振动信号;
步骤B、重复离线建模阶段第1步-第3步,生成新的样本数据xnew;
步骤C、根据第4步得到的Mexico小波核函数参数和第5步得到的核矩阵计算xnew对应的中心化核向量其中,K为第5步计算得到的核矩阵,表示与xnew相关的N×1列向量,k(xm,xnew)为Mexico小波核函数,xm就是SSKFDA模型学习时使用的所有标记和未标记样本集中第m个样本,m=1,2,…,N,I1×N=(1/N)1×N表示所有元素均为1/N的1×N行向量,IN×N=(1/N)N×N表示所有元素均为1/N的N×N矩阵,N为所有标记和未标记样本数量,上标T为转置;
根据第7步得到的投影矩阵,通过下面公式计算xnew鉴别特征向量ynew;
其中,φ(xnew)为xnew到高维特征空间的映射,满足k(xm,xnew)=φ(xnew)Tφ(xm)为Mexico小波核函数,φ(xm)为xm到高维特征空间的映射;
步骤E、判断BID是否超出预设阈值,如果超出阈值,则有新的未知故障产生并进行报警;否则机械运行处于已知状态,由下面规则判断目前机械所属已知类别状态;规则是:其中,表示ynew属于ωc的后验概率,p(ωc)为状态ωc的先验概率,c*为机械所属状态。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第2步具体如下:
乘积函数PF成分由一个包络信号和一个纯调频信号直接求出,PF成分包络信号的自相关函数Ra(·)第一个峰值Ra(A)对应的A被认为PF成分的周期,计算每个PF成分的周期调制强度PMI值,PMI=Ra(A)/(Ra(0)-Ra(A))表示包络信号能量与随机调制能量之比、用于度量有用的信号量;根据PMI值选取PF成分,其中Ra(0)表示PF成分的能量,Ra(A)被认为是PF成分中真实信号的能量。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第7步中,其中为第5步计算的流形间分离矩阵,表示第i个样本和第j个异类样本之间非相似性,Gb用于度量不同类别样本之间的区分性,为N×N的矩阵,其中,为第5步中的扩展最优类间样本重构权重向量,为对角矩阵,表示矩阵Gb中第i个行向量的重要程度,i=1,2,…,N;
其中为第5步得到的总体相似矩阵,是度量总体样本集中第i个样本和第j个样本之间相似性,矩阵Vpt=[p1,p2,…,pN]为N×N的矩阵,其中,pi为第5步中标记和未标记样本对样本重构权重向量,为对角距离,表示矩阵Gt中第i个行向量的重要程度,i=1,2,…,N;
表示特征空间同类数据之间的流形几何结构,其中为第5步得到的流形内相似性矩阵,是度量同类样本集中第i个样本和第j个样本之间相似性,为N×N的矩阵,其中,为第5步得到的扩展最优类内样本重构权重向量,为对角矩阵,表示矩阵Gw中第i个行向量的重要程度,i=1,2,…,N。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,步骤E中阈值设为0.95~0.99。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第1步中的CELMDAN算法具体如下:
初始化:PFs成分数量h=1,噪声数量L,rh-1=x表示从信号rh-1求取第h个PF成分,x为原始振动信号;
步骤①、首先使用局部均值分解LMD方法对L个白噪声信号w(l)进行分解,l=1,2,…,L,保存第h个PF成分Lh(w(l));
步骤②、原始振动信号与Lh(w(l))进行叠加生成L个新信号yl(t)=rh-1+AlLh(w(l)),其中幅值Al大小确定方法为Al=ε0std(x)/std(Lh(w(l))),std(x)表示原始振动信号x的标准差,std(Lh(w(l)))为Lh(w(l))的标准差,ε0∈[-1,1]为随机选取的加权系数;
步骤③、采用LMD算法对所有yl(t)进行分解,得到yl(t)的第1个PF及其相应的残差Mh(yl(t)),那么rh-1的第1个真实PF成分为PFh=rh-1-rh,其中为残差Mh(yl(t))的期望均值;判断rh是否为单频率震荡信号,如果rh是单频率震荡信号,执行步骤⑤;否则执行步骤④;
步骤④、h←h+1,以步骤③计算的rh-1作为需要进一步处理的信号,调用步骤①计算L个白噪声信号w(l)的第h个乘积函数PF成分Lh(w(l)),转到步骤②执行;
步骤⑤、保存对原始振动信号分解得到的PF成分PF1,PF2,…,PFh,其中,PFd′为第d′个PF成分,d′=1,2,…,h。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第4步中,求取最优墨西哥Mexico小波核函数参数的过程如下:
基于C类标记样本数据集X={X1,X2,…,XC},其中Xc表示第c类数据集,令X/Xc表示从X移出Xc后的数据集,对每个样本集{Xc,X/Xc},c=1,2,…,C;分别计算两个集合(Xc,Xc)、(X/Xc,X/Xc)以及(Xc,X/Xc)中样本之间距离,并按照升序方式分别对计算结果排序,令表示(Xc,Xc)样本之间距离升序的第αNc(Nc-1)个距离值、表示(X/Xc,X/Xc)样本之间距离升序的第α(N′-Nc)(N′-Nc-1)个距离值以及表示(Xc,X/Xc)中样本之间距离升序的第(1-α)Nc(N′-Nc)个距离值,比例系数α=η/8,η=1,2,…,7,N'为所有标记样本数量,Nc为第c类数据样本数量;
基于最小化样本集局部距离、最大化样本集非局部距离准则,达到同类样本之间尽可能靠近、异类样本之间尽可能远离,从而达到分开不同类别样本的目的,通过优化问题得到最优Mexico小波核函数参数σ*,
其中,为Mexico小波核函数,din是样本的维数,σ为核函数参数,函数diff(σ)表示样本集中数据样本之间的非局部距离与局部距离之差,表示类c样本与非类c之间距离的升序在α处距离值,表示类c样本之间距离的升序在α处距离值,表示非类c样本之间距离的升序在α处距离值,为特征空间中类c样本与其他类样本之间的距离,表示特征空间上异类样本之间的非局部距离;为特征空间上类c的样本之间的距离,表示特征空间上类c样本集的局部距离;为在特征空间上非类c样本集的局部距离,表示非类c样本集的局部距离;令γ=1/2σ2,γ为中间变量,由最优中间变量γ*很容易根据求出最优核函数参数σ*, 的含义分别与和的含义完全相同,那么上述优化问题使用下面梯度下降法求解;
梯度下降更新形式为
其中,τ为步长,γg+1、γg表示第g+1次和第g次迭代的γ值,g为迭代次数,diff(γg)第g次迭代时的值,为第g次迭代时特征空间中类c样本与其他类样本之间的距离,为第g次迭代时特征空间上类c样本集的局部距离,为第g次迭代时特征空间上非类c样本集的局部距离;为第g次迭代时在特征空间数据集的非局部距离与局部距离差值,
其中,E=(1/N)N×N为N×N的元素全部1/N的矩阵,核矩阵其中k(xm',xn')=(din-||xm'-xn'||2/σ2)exp(-||xm'-xn'||2/(2σ2))为样本xm'和xn'对应的Mexico小波核函数值,din为样本的维数,m',n′=1,2,…,N。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第5步中,
其中,φ(xm')表示与Mexico核函数相关的高维映射函数φ(·)对数据样本xm'的映射,为第5步中的流形内相似矩阵的元素,W为投影矩阵,上标T为转置,由进行填0扩展得到并满足 为使用同类样本对样本的重构向量,矩阵B=[B1,B2,…,BC,BU],其中B表示所有标记样本和未标记样本在高维特征空间映射组成的矩阵,表示类c所有样本在高维特征空间的映射矩阵,c=1,2,…,C,表示所有未标记样本在高维特征空间的映射所组成的矩阵,Nu为未标记样本数量,分别表示样本和到高维特征空间的映射,表示类c数据样本,b=1,2,…,Nc,Nc为类c数据样本的数量,为标记数据样本,b'=1,2,…,Nu,Nu为未标记数据样本的数量。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第5步中,
其中,表示流形间分离矩阵Gb∈RN×N的元素,φ(xm')表示与Mercer核函数对应的高维特征映射,rm′是由进行填0扩展得到的,并满足其中为使用非类c样本对样本的重构向量,表示非类c样本在高维特征空间的映射所组成的矩阵。
作为本发明所述的一种基于CELMDAN和SSKFDA的机械状态监测方法进一步优化方案,第6步,
Srb=(1-β)Sb+βSt
Srw=(1-β)Sw+βtr(WTW)
其中,总体数据集散度St表示所有样本在空间的分布,β∈[0,1]为正则化因子;当β=1时,SSKFDA退化为KPCA方法;当β=0时,SSKFDA退化为KFDA方法。
本发明采用以上技术方案与现有技术相比,具有以下技术效果:
(1)使用CELMDAN把含噪调制信号分解为多个单模态解调成分-乘积函数(productfunctions,PFs),消除传统LMD中模式混淆问题,PFs不受背景噪声干扰、能够更好地表征故障特征信息。
(2)提出周期调制强度(periodicmodulation intensity,PMI)选择PFs方法,该方法良好的鲁棒性、能量独立性、物理意义明确以及易于计算等优点,能够选择准确可靠的PFs;
(3)将核FDA方法扩展到半监督学习和压缩感知框架,提出SSKFDA(semi-supervised sparse kernel FDA,SSKFDA)算法,该方法充分利用标签样本和未标签集几何信息以及稀疏性把感兴趣的PFs多域特征约简到低维空间,更好地揭示嵌入在高维稀疏空间的低维子空间数据集,更好的刻画复杂数据集分布;
(4)提出周期调制强度(periodic modulation intensity,PMI)作为选择PFs准则的方法,解决PFs数量多于实际信号成分时导致部分PFs的有效信息量过少的问题;PMI利用故障振动信号具有周期调制特点,具有良好的鲁棒性、明确的物理意义以及尺度不变性,能够准确选择有效的PFs;
(5)提出一种快速SSKFDA模型选择方法,该方法通过基于类内局部结构测度最小、全局部结构测度最大准则求取最优模型参数;
(6)提出全局监测统计量与Bayesian后验推理的已知与未知机械状态检测方法,克服现有故障检测方法无法区分已知或未知故障的问题。
附图说明
图1是ELMD算法实现图。
图2是CELMDAN算法实现图。
图3是Bayesian推理的CELMDAN与SSKFDA机械状态监测流程图。
具体实施方式
下面结合附图对本发明的技术方案做进一步的详细说明:
本发明考虑实际标签样本数量往往远少于未标签样本数量,高维数据往往是稀疏的,以及CELMDAN对复杂调制信号良好的解调能力,本发明提出Bayesian框架下CELMDAN和半监督核稀疏FDA(semisupervised sparse kernel Fisher discriminant analysis,SSKFDA)的机械状态监测方法。该方法使用CELMDAN分为如下部分:(1)使用CELMDAN把含噪调制信号分解为多个单模态解调成分-乘积函数(product functions,PFs),消除传统LMD中模式混淆问题,PFs不受背景噪声干扰、能够更好地表征故障特征信息。(2)提出周期调制强度(periodic modulation intensity,PMI)选择PFs方法,该方法良好的鲁棒性、能量独立性、物理意义明确以及易于计算等优点,能够选择准确可靠的PFs;(3)在半监督学习和压缩感知框架对核Fisher鉴别方法KFDA方法进行扩展,提出SSKFDA(semi-supervisedsparse kernel FDA,SSKFDA)算法。该方法充分利用标签样本和未标签集信息把感兴趣的PFs多域特征约简到低维空间,有效地揭示嵌入在高维稀疏流形空间的低维子空间数据集,更好的刻画复杂数据类别分布;(3)提出周期调制强度(periodic modulation intensity,PMI)作为选择PFs准则的方法,解决PFs数量多于实际信号成分时导致出现部分冗余PFs问题。PMI利用故障振动信号具有周期调制特点,具有良好的鲁棒性、明确的物理意义以及尺度不变性,能够准确选择有效的PFs。(4)提出一种快速SSKFDA核函数参数选择方法。该方法通过基于类内局部结构测度最小、全局部结构测度最大准则求取最优模型参数。(5)提出基于贝叶斯推理的全局监测统计量的已知与未知状态检测方法,较好的考虑了运行状态的不确定性,克服现有机械故障检测方法无法区分已知或未知故障的问题。
相关技术
CELMDAN
LMD能够把复杂调制信号分解到一系列包络信号与纯频率调制信号乘积函数。LMD不仅实现幅值调制AM-频率调制FM信号的解调,PFs的瞬时频率和复制较好地描述原始信号复杂时频分布。对给定原始信号x'(t),LMD算法过程如下:
1)确定原始信号x'(t)的所有局部极值点nei′,计算相邻局部极值点nei′与nei′+1的平均值并使用直线连接起来,计算相应包络然后使用滑动平均法分别对所有mai′和eai′平滑估计局部均值函数ma′11(t)和包络估计函数ea′11(t)。
2)从信号x′(t)中分离出局部均值函数ma′11(t),可得残差rh11(t)
rh11(t)=x′(t)-ma′11(t) (1)
使用下式对残差信号rh11(t)进行解调,即
ds11(t)=rh11(t)/ea′11(t) (2)
3)判断解调信号是否为纯调频信号。如果包络函数ea′12(t)应满足条件|ea′12(t)-1|≤δ′,那么ds11(t)为纯调频信号,这里δ′为事先设定的阈值。如果ea12(t)不满足上述条件,重复步骤1)-2)将ds11(t)继续分解直至ds1F(t)成为满足-1≤ds1F(t)≤1的纯调频信号,即ds1F(t)的包络满足|ea′1(F+1)(t)-1|≤δ′。相关的残差信号rh1i(t)和纯调频信号ds1i(t)可以通过下式求取
rh11(t)=x′(t)-ma′11(t)
ds1i′(t)=rh1i′(t)/ea′1i′(t)
rh1(i′+1)(t)=ds1i′(t)-ma′1(i′+1)(t)i′=1,2,…,F
4)计算瞬时幅值函数ia1(t)
根据纯调频函数和上面瞬时幅值信号计算乘积函数PF1(t)
PF1(t)=ea1(t)ds1F(t) (4)
这里,PF1(t)是一个单成分AM-FM信号,包含原始信号的最高频率震荡信号。由上式可计算该成分的瞬时频率为
根据式(3)和(5)很容易求取调频分量PF1(t)的瞬时幅值和瞬时频率。
5)从原始信号x′(t)减去PF1(t)得
u1(t)=x(t)-PF1(t) (6)
对u1(t)信号重复步骤1)-4)求取单成分AM-FM信号PF2(t),然后对信号u2(t)=u1(t)-PF2(t)作为原始信号重复步骤1)-4)求取单成分AM-FM信号PF3(t)。以此类推,经过H次循环,直至uH(t)=uH-1(t)-PFH(t)成为单调函数时停止算法。这样,原始信号x′(t)可表示为B′个单分量PF与一个单调函数uB‘(t)之和,即
这样复杂信号的时频分布可以由所有PF成分的瞬时幅值和频率描述。
针对LMD存在的端点效应问题,目前主要有信号延拓法、波形匹配法、极点延拓法、回归方法等方法。针对LMD存在的模式混淆问题,ELMD采用类似EEMD思想,通过对原始信号施加多个不同幅度白噪声信号,然后使用LMD对这些含噪信号分解到多个调频成分,最后计算相应尺度上调频成分的均值作为最终PFs,ELMD实现如图1所示。
需要指出的是,ELMD仍然存在如下问题:(1)ELMD把相同尺度的信号成分以及加性噪声部分均被分解到相应的PFs,需要多次集成试验消除加性噪声部分,增加了算法的计算量;(2)施加噪声幅值对ELMD性能影响很大。施加噪声幅值太小不足以影响极点的分布,因而无法解决模式混淆问题。而噪声幅值太大导致模式包含相当大的残差,需要更多的集成试验次数;(3)不同高斯噪声加入原始信号,不同数量试验会产生不同的PFs,而不同PFs的物理意义也会有所区别,不利于系统运行状态的识别。为克服上述问题,Lei Wang(2017)提出噪声自适应完全ELMD(complete ELMD with adaptive noise,CELMDAN)算法。在每次试验中,该方法在每个分解阶段均加入一个特殊的自适应噪声,每次分离PF后会得到唯一的残差。然后把残差作为下一次分解的输入。
CELMDAN把真实PF定义为当前信号与残差均值之间的差值。CELMDAN在不同分解阶段把白噪声PF施加到信号,使得与施加噪声的相似尺度上信号成分被滤波为PF,因此残差含有很少的施加噪声。CELMDAN抽取的PFs具有更加合适的尺度,包含更少的残差噪声,并并且在每个阶段中降低筛选迭代次数,其实现如图2所示。当CELMDAN应用于复杂振动数据分析时,PFs的数量往往多余实际信号成分数量,导致一些PFs含有很少的信息量。基于相关分析和Kurtosis选取有效的PFs往往受到信号噪声的影响。正常机械振动往往呈现一定的周期性,即使机械设备发生故障时,机械测量信号也呈现周期幅值调制信号。针对这种现象,本发明使用周期调制强度(periodic modulation intensity,PMI)作为度量PFs信息量准则。PMI具有良好的鲁棒性、尺度不变性、物理意义明确等优点。
PMI
由上面信号仿真分析可知,CELMDAN能够对复杂信号进行有效分析。然而,当分析信号含有较强噪声时,由CELMDAN得到的部分PFs含有诊断信息量很少,需要移除这些冗余PFs。注意到机械设备发生故障时,机械振动信号能量将发生周期性变化,从而导致测量信号周期幅值调制(periodic amplitude modulation,PAM)现象。因此可以使用PAM评估机械健康程度。PAM一般由原始信号的包络进行辨识。假设原始信号为x'(t),其包络信号为
其中,i”表示虚部单位,算子H表示Hilbert变换。故障特征的检测性不仅依赖PAM本身能量,而且还依赖由测量噪声和其他干扰引起的随机调制信号能量。因此对应包络信号en(t)一般由PAM信号pm(t)和不相关随机调制rm(t)组成。周期调制强度(periodicmodulation intensity,PMI)定义为PAM本身能量Epm与随机调制能量Erm之比,可以用于作为度量信号的信息量,其形式如下
显然PMI值越大意味着信号特征越明显。注意到pm(t)与nm(t)的不相关性,很容易推导出自相关和Epm与Erm的关系
Ren(0)=Epm+Erm (10)
其中,Ren(·)表示en(t)的自相关函数,Ren(0)表示en(t)的能量。包络信号PAM部分一般具有周期性,Epm可以通过检测包络自相关信号相应峰值估计,即Ren(T)=Epm,T为周期。这样PMI可以由下式估计
PMI具有能量独立性,明确的物理意义、良好的鲁棒性以及易于计算等优点,可以使用PMI选择重要的PFs成分。
故障特征抽取
振动信号的复杂度和各个频段上能量分布能够在一定程度上反映振动信号的规则程度以及特征,本发明从含有诊断信息量大PFs中抽取如下的信息熵特征进行故障类型识别故障诊断的重要步骤,这些信息熵计算方法如下:
1)PF奇异熵(product function singular entropy,PFSE)
一般而言,工业过程数据的奇异值具有较强的鲁棒性,能够表征过程状态,因此该方法引入到机械故障诊断中。假设选取Mpf个PFs,这些PFs奇异值生成乘积函数奇异熵(product function singular entropy,PFSE)特征。令为CELMDAN对信号分解得到的第i个PF,为由构成的矩阵。对进行奇异值分解并有G个非负奇异值Λj”(j”=1,2,…,G),那么PFSE定义为
PFSE描述了信号的频率成分和分布特点。
2)PF时频熵(product function time-frequency entropy,PFTFE)
那么PFTFE定义为
3)PF能量熵(product functions energy entropy,PFEE)
乘积函数能量熵的定义为
1半监督稀疏核Fisher鉴别分析
在实际应用中,获取数据标签代价往往较为昂贵,导致未标签数据样本数量远远大于标签样本,有限标签样本数量往往不足以充分反映样本的总体分布,导致监督学习算法性能不能满意的问题。另外,样本数据具有非线性、高维、稀疏特点,高维数据集往往位于低维流形几何上。半监督学习能够充分利用未标签样本和标签样本数据信息,训练的学习机模型能够更好地刻画数据类别的几何分布,能够充分利用数据鉴别信息,因而能够取得较好的效果,受到人们的广泛关注。稀疏表示是计算机视觉和模式识别领域的一个研究热点,测试样本被表示成所有训练样本的线性组合,并且要求表示系数是稀疏的。基于核稀疏表示、基于局部保持的稀疏表示子空间特征抽取方法优于流形学习子空间特征抽取算法,具有更强的鉴别能力。把稀疏表示与流形学习、PCA、LDA等其他机器学习算法深度融合,提高传统机器学习算法的性能,成为机器学习研究的热点问题。
受到SMMSD算法、核稀疏鉴别分析和半监督学习的启发,本发明提出半监督核稀疏线性鉴别分析,把非线性高维、稀疏数据集投影到低维鉴别子空间中,提高不同机械状态数据的分离性。下面给出SSKFDA的详细描述。
设标签样本数据集XL和未标签样本数据集XU,全体数据集记为X=[XL,XU]。假设有C类数据集,第c类样本数据集记为其中,Nc为第c类样本数量,c=1,2,…,C,首先通过非线性映射φ(·)把样本数据映射到高维特征空间。令对每个训练样本把从Bc移出后的样本,然后使用Bc的其余样本集对进行线性表示。假设样本数据对应XL中的数据xm',那么重构权重通过下面基于l1-范数的最优重构问题求得
这里,求解上述优化问题后,对进行填0扩展使得其中,B=[B1,B2,…,Bc,BU],其中表示所有未标记样本在高维特征空间的映射所组成的矩阵,Nu为未标记样本数量,B表示所有标记样本和未标记样本在高维特征空间映射组成的矩阵。假设给定投影矩阵W,根据稀疏鉴别嵌入(sparse discriminant embedding,SDE)理论和流形学习理论,类内散度Sw被定义样本在原始高维空间到低维嵌入空间中被同类样本稀疏重构残差量,其定义如下
这里,l(xm')表示样本xm'的类别标签。对每个有标记训练样本设数据集表示从XL移出c-th类样本集后标签数据集,假设样本数据对应XL中的数据xm',那么由样本集对φ(xm')线性表示的最优权重可通过求解下面优化问题
这里, 这里p'=1,2,…,C,p'≠c。对进行填充0扩展为满足条件BL为由映射φ(·)把XL映射到高维特征空间中的数据集。令R=[R1,R2,…,RC,RU]∈RN×N,类间散度Sb度量在原始高维空间到低维嵌入空间中样本与异类样本重构误差,其定义形式为
这里,样本类标签为C+1时表示该样本为未标签样本。那么对样本φ(xm')∈B,使用φ(xm')从B移出后剩余样本线性表示φ(xm'),那么重构权重向量pi可以通过求解下面的优化问题求得
那么总体数据集散度St定义为
Sw=tr(WTBSsBTW) (26)
Sb=tr(WTBSrBTW) (27)
St=tr(WTBSpBTW) (28)
这里,Ss、Sr、Sp形式与Sv类似,在此不再赘述。半监督稀疏核LDA(semisupervisedsparse kernel LDA,SLDA)模型的正则化类间散度矩阵、类内散度矩阵定义为
Srb=(1-β)Sb+βSt (29)
Srw=(1-β)Sw+βtr(WTW) (30)
这里,It表示合适维数的单位矩阵;β∈[0,1]为正则化因子。当β=1时,SSKFDA退化为KPCA方法;当β=0时,SSKFDA退化为KFDA方法。类似于KFDA方法,半监督投影鉴别矩阵就在由下面的SSKFDA优化问题求取:
其最优解可以归结如下广义特征值问题;
((1-β)BSrBT+βBSpBT)w=λ((1-β)BSsBT+βI)w (32)
这里,λ是广义特征值,投影向量w为相应的特征向量。根据核技巧,投影向量w可以表示成w=Bq,那么式(27)转化为如下形式
BTB((1-β)Sr+βSp)BTBq=λBTB((1-β)Ss+βI)BTBq (33)
注意到K=BTB=[k(xm',xn')]m',n'=1,2,…,N,那么上式转化如下形式
((1-β)Sr+βSp)Kq=λ((1-β)Ss+βI)Kq (34)
假设式(34)对应的广义特征值降序排列记为λ1≥λ2≥…≥λm,相应的特征向量q1,q2,…,qm作为SSKFDA模型的方向。如果选择前面特征向量Qr=[q1,q2,…,qr]作为鉴别投影矩阵,则Wr=BQr,新的数据样本x通过下面公式计算鉴别向量
这里,kx=[k(x1,x),k(x2,x),…,k(xN,x)]T。
这里,I1×N=(1/N)1×N,IN×N=(1/N)N×N。
2 SSKFDA核函数参数选择
众所周知,核学习机性能受到核函数参数和惩罚系数的严重影响。目前,核学习机模型参数选择方法主要有网格搜索法、粒子群优化方法、遗传算法、交叉验证、留一法(leave one out,LOO)。这些方法虽然能够获得最优的模型参数,但是计算复杂度很高,不适合大数据量的核学习机训练情况,限制其应用范围。针对核学习机模型参数选择问题,Zhongben Xu(2009)基于高斯核函数具有数据结构保持性质以及多尺度高斯核函数能够描述特征空间全局和局部数据结构思想,提出基于局部结构测度最小、全局结构测度尽可能大准则的最优SVM核函数参数求取方法;基于支持向量值能够度量相应样本的分类正确程度,提出通过限制支持向量值迭代求取最优惩罚参数。Matthias Varewyck(2011)提出一种由维数确定高斯核函数参数方法,并能够使用2-3次方法确定合适的惩罚参数。
Mexico小波核函数ψ(xi,xj)=ψ(||xm'-xn'||)=(din-||xm'-xn'||2/σ2)exp(-||xm'-xn'||2/(2σ2))是一种满足Mercer条件的核函数,具有平移不变性质,也称为径向小波函数,这里din为输入变量的维数,σ为Mexico小波核函数参数。对于输入了空间xm'和xn'的欧氏距离为||xm'-xn'||,其在特征空间的欧氏距离为
径向基小波函数不仅具有更小的支撑区间,不同尺度参数的径向小波SVM能够地把数据集滤波到不同频段,在不同尺度上观测数据特性。比如,小尺度径向基小波SVM分类平面非常复杂,能够高度拟合训练数据样本标签,即每个训练数据样本可以正确分开。由于特征空间的高维映射数据集结构无法保持原始空间数据集局部结构,导致模型泛化性能变差。大尺度径向基小波SVM分类平面非常平滑,特征空间数据对距离趋于0,难以反映输入空间全局数据结构,导致逼近性能很差以及很强的分类器泛化性能(尽管泛化性没有意义)。综合以上分析以及FDA思想,选择合适的核函数参数使得映射后特征数据集仍然能够保持原始数据局部和全局(非局部)结构,同时全局(非局部)测度尽可能大、局部测度尽可能小,提高核函数的数据鉴别能力。基于上述思想,我们提出一种快速的小波SSKFDA模型选择方法。
假设给定C类数据样本集X={X1,X2,…,XC},Xc表示第c类数据集,第c类的数据量记为Nc,且表示标记样本数量,令X/Xc表示从X移出Xc后的数据集。对数据样本集{Xc,X/Xc},计算Xc与Xc中的样本对距离,X/Xc与X/Xc样本对距离,以及Xc与X/Xc数据样本对距离,并分别对上述计算结果按照升序方式进行排列(删除掉0-距离),选取αNc(Nc-1)位置的距离和α(N'-Nc)(N'-Nc-1)位置的距离作为局部结构测度,在(1-α)Nc(N'-Nc)位置的距离作为非局部测度,这里c=1,2,…,C,比例系数α=l/8,η=1,2,…,7。根据定义可知,表示第c类数据集的局部测度,表示非第c类数据集的局部测度,表示第c类数据与异类数据集之间非局部测度。显然在特征空间内,我们希望同类映射数据之间尽可能靠近、非同类数据之间尽可能远离,与此同时同类数据与异类数据之间尽可能远离。一般说来,我们不建议选取a值接近0。原因是如果a取值接近0,数据集的局部和非局部结构测度分别为最小、最大,这两类值差别最大,数据噪声或者在野点会影响基于局部测度与非局部测度的模型选择策略的有效性。最优径向基小波核函数参数σ*应该使得局部测度和非局部测度差别最大,即,即和尽可能小、尽可能大。根据上述思想,通过求取如下优化问题得到最优Mexico小波核函数参数σ*
其中,为Mexico小波核函数,din是样本的维数,σ为核函数参数,函数diff(σ)表示样本集中数据样本之间的非局部距离与局部距离之差。上述优化问题具有明确的物理意义,最优模型参数值能够使得特征空间中类内数据集尽可能靠近、异类数据集尽可能远离,因此相应的模型对数据具有更好的鉴别分析能力,提高数据维数约简性能。显然上述优化问题为非凸优化问题,我们使用梯度下降法求取优化问题的最优解。为简单起见,令中间变量γ=1/2σ2,由最优中间变量γ*很容易根据求出最优核函数参数σ*, 的含义分别与和的含义完全相同,那么上述优化问题使用下面梯度下降法求解;
梯度下降更新形式为
Mexico小波核函数参数初始值由下式确定,即
其中,dst表示所有样本之间距离的平均值,所有数据样本均进行了归一化处理。梯度下降算法终止条件设为|diff(γg)-diff(γg+1)|≤10-5。
3基于SSKFDA的Bayesian推理的机械状态监测
对新生成的样本数据xnew,新提出的故障分类方法第一步需要把该样本映射到高维特征空间φ(xnew),然后把φ(xnew)投影到W中的各个方向上得到鉴别向量
然后基于ynew使用监测统计量确定新样本所属过程状态类别。这里提出构建Bayesian推理的故障监测统计量实现对多个机械运行状态进行监测。
对低维约简空间的新样本数据ynew,其关于类ωc的条件概率记为p(ynew|ωc),先验概率记为p(ωc)=Nc/NL,基于Bayesian推理理论,ynew属于类ωc的后验概率为
一般情况下,SSKFDA的隐藏变量ynew往往是高斯分布。因此,每类样本的条件概率分布使用多变量高斯分布估计
这里和Ξc分别为属于类ωc的样本集均值与协方差,dout为ynew的维数。除了现有故障状态外,机械在运行过程中往往还会出现未知的故障状态,因此需要构建统计量判断新样本是否属于新的运行状态。对新的样本数据ynew,其对应的全局统计量定义为
其中Ξc由分别表示低维约简空间中第c类样本数据集的均值和协方差矩阵。于每类数据服从单一模态高斯分布,那么从一个高斯成分采样得到的样本的Mahalanobis距离服从分布。因此,BID控制上限在置信度水平γ'下服从自由度为fr的χ2分布。在实际应用中,γ'通常取95%或99%。当新样本出现时,其对应的监测统计量BID,状态判断准则:
基于以上所述,基于CELMDAN和SSKFDA的Bayesian推理的机械状态监测实现流程图如图3所示,本发明方法应用于机械状态监测分为在线监测阶段和离线建模阶段。
离线建模阶段:
初始化阶段:从旋转机械采集故障和正常状态的振动信号,类内、类间散度矩阵的正则化因子β,核函数参数初始值σ,SSKFDA算法约简后数据特征维数r,置信度水平γ'。
1、对振动信号归一化,根据上面介绍LMD理论以及图2给出的CELMDAN算法对归一化振动信号进行分解分为多个PFs成分;
2根据上面计算每个PF成分PMI指标,根据PMI值并选取重要的PF成分;
3从选取的PFs抽取上面介绍的特征(在上面的故障特征抽取中提到的),生成标记样本数据集和未标记数据样本集;
4基于标记样本数据集,根据上面介绍的SSKFDA模型选择算法求取最优Mexico小波函数参数;
6基于第5步计算结果,根据式(29)、(30)构建正则化类内散度矩阵和类间散度矩阵Srw,Srb;
7计算式(34)所示特征值和特征向量,选取前r个最大特征值对应的特征向量作为的投影矩阵;
8保存求取的SSKFDA模型;
在线监测阶段:
1从运行的旋转机械采集振动信号;
2重复离线建模阶段步1-步4生成样本数据;
3根据式(37)对中心化新样本对应的核向量,式(43)计算投影鉴别向量;
4根据式(46)计算监测统计量BID;
5判断BID是否超出阈值,如果超出阈值,则有新的未知故障产生并进行报警;否则机械运行处于已知状态,由式(47)确定目前机械所属状态。
图1是ELMD算法实现图,具体过程如下:
1、输入原始信号x,初始化施加白噪声次数l'=1;
2、施加不同幅值的白噪声nw后生成新的信号yl'=x+Bl'nw,幅值Bl'可以从正态分布中随机采样确定;
3、然后利用LMD算法把信号yl'分解成多个不同频段的乘积函数,记为PFIl',d',d'=1,2,…,D',l'=l'+1;
4、如果l'<D',转到第2步;否则对相频段的PFs信号进行平均:达到消除PF信号中的白噪声、提高LMD算法鲁棒性的目的。注意,基于matlab语言的LMD算法源程序已经公开,根据图1所示步骤很容易实现ELMD算法。
如图2的CELMDAN算法具体过程如下:
初始化:PFs成分数量h=1,噪声数量L,rj-1=x,x为原始信号
第1步、首先使用LMD程序对L个白噪声信号进w(l)进行分解,保存第h个PF成分Lh(w(l));
第2步、原始信号与变幅值Lh(w(l))进行叠加生成L个新信号yl(t),其中幅值大小确定方法为Al=ε0std(x)/std(Lh(w(l))),std(·)为标准差算子;
第4步、以第3步计算的残差均值rj作为原始信号,h←h+1,调用第1步计算白噪声信号的第h个PF成分Lh(w(l)),继续使用第2、3步得到第2个残差均值rj
第5步、按照第4步方法依次计算原始信号真实PF成分,直至残差为单频率震荡信号结束。
CELMDAN算法实现的基础仍然是LMD,因此使用LMD程序很容易实现CELMDAN。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围内。
Claims (10)
1.一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,包括离线建模阶段和在线监测机械状态阶段;具体步骤如下:
离线建模阶段:
初始化阶段:从旋转机械采集故障和正常状态的振动信号,类内、类间散度矩阵的正则化因子β,核函数参数初始值σ,半监督系数核Fisher鉴别分析SSKFDA算法约简后数据特征维数r,置信度水平γ',每个状态的先验概率;正常状态指非故障状态;
第1步、对振动信号归一化为[0,1]内,根据局部均值分解LMD算法以及噪声自适应完全集成局部均值分解CELMDAN算法对归一化后的振动信号进行分解,分为多个瞬时频率具有物理意义的乘积函数PF成分;
第2步、计算每个PF成分的周期调制强度PMI值,根据PMI值选取PF成分;
第3步、从第2步得到的所有PF成分中抽取PF奇异熵特征、PF时频熵特征、PF能量熵特征,分别从每一个选取的PF成分中抽取峰度Kurtosis、偏度Skewness、标准差和脉冲因子统计特征;根据PF奇异熵特征、PF时频熵特征、PF能量熵特征、抽取峰度Kurtosis、偏度Skewness、标准差和脉冲因子统计特征,生成标记样本数据集和未标记样本数据集;
第4步、基于标记样本数据集,首先根据特征空间中类内数据集尽可能靠近、类间数据集尽可能远离准则,求取最优墨西哥Mexico小波核函数参数;
第5步、根据所有样本数据集以及第4步求取的最优墨西哥Mexico小波核函数参数计算核矩阵K,根据核矩阵K计算中心化核矩阵利用扩展最优类内样本重构权重向量流形内相似性矩阵Gw∈RN×N计算类内散度Sw,利用流形间分离矩阵Gb∈RN×N、扩展最优类间样本重构权重向量计算类间散度Sb,利用总体相似矩阵Gt∈RN×N、标记和未标记样本对样本重构权重向量pi以及构建总体散度St;RN×1表示N×1的实数列向量,RN×N表示N×N的实数矩阵,N为标记样本和未标记样本的总数量;
第6步、基于第5步计算结果,构建正则化类内散度Srw和类间散度Srb,其中,Srw度量同类样本之间散度,Srb度量不同异类样本之间散度;
第7步、基于第5步和第6步计算结果,根据Srw最小、Srb最大准则,计算所示特征值和特征向量,选取前r个最大特征值对应的特征向量作为投影矩阵W=[q1,q2,…,qr];β∈(0,1)为正则化因子,q为特征向量,λ为特征值,I为N×N的单位矩阵,qp为第p个最大特征值对应的特征向量,p=1,2,…,r,r为选取特征向量的数量,Srb'表示特征空间不同类别数据之间的流形几何结构,Spt'表示特征空间全体数据之间的流形几何结构,Ssw'表示特征空间同类数据之间的流形几何结构;
第8步、保存投影矩阵、Mexico小波核函数参数、核矩阵作为求取的SSKFDA模型;
在线监测阶段:
步骤A、从运行的旋转机械采集振动信号;
步骤B、重复离线建模阶段第1步-第3步,生成新的样本数据xnew;
步骤C、根据第4步得到的Mexico小波核函数参数和第5步得到的核矩阵计算xnew对应的中心化核向量其中,K为第5步计算得到的核矩阵,kxnew=[k(x1,xnew),k(x2,xnew),…,k(xN,xnew)]T表示与xnew相关的N×1列向量,k(xm,xnew)为Mexico小波核函数,xm就是SSKFDA模型学习时使用的所有标记和未标记样本集中第m个样本,m=1,2,…,N,I1×N=(1/N)1×N表示所有元素均为1/N的1×N行向量,IN×N=(1/N)N×N表示所有元素均为1/N的N×N矩阵,N为所有标记和未标记样本数量,上标T为转置;
根据第7步得到的投影矩阵,通过下面公式计算xnew鉴别特征向量ynew;
其中,φ(xnew)为xnew到高维特征空间的映射,满足k(xm,xnew)=φ(xnew)Tφ(xm)为Mexico小波核函数,φ(xm)为xm到高维特征空间的映射;
2.根据权利要求1所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第2步具体如下:
乘积函数PF成分由一个包络信号和一个纯调频信号直接求出,PF成分包络信号的自相关函数Ra(·)第一个峰值Ra(A)对应的A被认为PF成分的周期,计算每个PF成分的周期调制强度PMI值,PMI=Ra(A)/(Ra(0)-Ra(A))表示包络信号能量与随机调制能量之比、用于度量有用的信号量;根据PMI值选取PF成分,其中Ra(0)表示PF成分的能量,Ra(A)被认为是PF成分中真实信号的能量。
3.根据权利要求1所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第7步中,其中为第5步计算的流形间分离矩阵,表示第i个样本和第j个异类样本之间非相似性,Gb用于度量不同类别样本之间的区分性,为N×N的矩阵,其中,为第5步中的扩展最优类间样本重构权重向量,为对角矩阵,表示矩阵Gb中第i个行向量的重要程度,i=1,2,…,N;
其中为第5步得到的总体相似矩阵,是度量总体样本集中第i个样本和第j个样本之间相似性,矩阵Vpt=[p1,p2,…,pN]为N×N的矩阵,其中,pi为第5步中标记和未标记样本对样本重构权重向量,为对角距离,表示矩阵Gt中第i个行向量的重要程度,i=1,2,…,N;
4.根据权利要求1所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,步骤E中阈值设为0.95~0.99。
5.根据权利要求1中所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第1步中的CELMDAN算法具体如下:
初始化:PFs成分数量h=1,噪声数量L,rh-1=x表示从信号rh-1求取第h个PF成分,x为原始振动信号;
步骤①、首先使用局部均值分解LMD方法对L个白噪声信号w(l)进行分解,l=1,2,…,L,保存第h个PF成分Lh(w(l));
步骤②、原始振动信号与Lh(w(l))进行叠加生成L个新信号yl(t)=rh-1+AlLh(w(l)),其中幅值Al大小确定方法为Al=ε0std(x)/std(Lh(w(l))),std(x)表示原始振动信号x的标准差,std(Lh(w(l)))为Lh(w(l))的标准差,ε0∈[-1,1]为随机选取的加权系数;
步骤③、采用LMD算法对所有yl(t)进行分解,得到yl(t)的第1个PF及其相应的残差Mh(yl(t)),那么rh-1的第1个真实PF成分为PFh=rh-1-rh,其中为残差Mh(yl(t))的期望均值;判断rh是否为单频率震荡信号,如果rh是单频率震荡信号,执行步骤⑤;否则执行步骤④;
步骤④、h←h+1,以步骤③计算的rh-1作为需要进一步处理的信号,调用步骤①计算L个白噪声信号w(l)的第h个乘积函数PF成分Lh(w(l)),转到步骤②执行;
步骤⑤、保存对原始振动信号分解得到的PF成分PF1,PF2,…,PFh,其中,PFd'为第d'个PF成分,d'=1,2,…,h。
6.根据权利要求1中所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第4步中,求取最优墨西哥Mexico小波核函数参数的过程如下:
基于C类标记样本数据集X={X1,X2,…,XC},其中Xc表示第c类数据集,令X/Xc表示从X移出Xc后的数据集,对每个样本集{Xc,X/Xc},c=1,2,…,C;分别计算两个集合(Xc,Xc)、(X/Xc,X/Xc)以及(Xc,X/Xc)中样本之间距离,并按照升序方式分别对计算结果排序,令表示(Xc,Xc)样本之间距离升序的第αNc(Nc-1)个距离值、表示(X/Xc,X/Xc)样本之间距离升序的第α(N'-Nc)(N'-Nc-1)个距离值以及表示(Xc,X/Xc)中样本之间距离升序的第(1-α)Nc(N'-Nc)个距离值,比例系数α=η/8,η=1,2,…,7,N'为所有标记样本数量,Nc为第c类数据样本数量;
基于最小化样本集局部距离、最大化样本集非局部距离准则,达到同类样本之间尽可能靠近、异类样本之间尽可能远离,从而达到分开不同类别样本的目的,通过优化问题得到最优Mexico小波核函数参数σ*,
其中,为Mexico小波核函数,din是样本的维数,σ为核函数参数,函数diff(σ)表示样本集中数据样本之间的非局部距离与局部距离之差,表示类c样本与非类c之间距离的升序在α处距离值,表示类c样本之间距离的升序在α处距离值,表示非类c样本之间距离的升序在α处距离值,为特征空间中类c样本与其他类样本之间的距离,表示特征空间上异类样本之间的非局部距离;为特征空间上类c的样本之间的距离,表示特征空间上类c样本集的局部距离;为在特征空间上非类c样本集的局部距离,表示非类c样本集的局部距离;令γ=1/2σ2,γ为中间变量,由最优中间变量γ*很容易根据求出最优核函数参数σ*, 的含义分别与和的含义完全相同,那么上述优化问题使用下面梯度下降法求解;
梯度下降更新形式为
其中,τ为步长,γg+1、γg表示第g+1次和第g次迭代的γ值,g为迭代次数,diff(γg)第g次迭代时的值,为第g次迭代时特征空间中类c样本与其他类样本之间的距离,为第g次迭代时特征空间上类c样本集的局部距离,为第g次迭代时特征空间上非类c样本集的局部距离;为第g次迭代时在特征空间数据集的非局部距离与局部距离差值,
8.根据权利要求1中所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第5步中,
其中,φ(xm′)表示与Mexico核函数相关的高维映射函数φ(·)对数据样本xm'的映射,为第5步中的流形内相似矩阵的元素,W为投影矩阵,上标T为转置,由进行填0扩展得到并满足 为使用同类样本对样本的重构向量,矩阵B=[B1,B2,…,BC,BU],其中B表示所有标记样本和未标记样本在高维特征空间映射组成的矩阵,表示类c所有样本在高维特征空间的映射矩阵,c=1,2,…,C,表示所有未标记样本在高维特征空间的映射所组成的矩阵,Nu为未标记样本数量,分别表示样本和到高维特征空间的映射,表示类c数据样本,b=1,2,…,Nc,Nc为类c数据样本的数量,为标记数据样本,b'=1,2,…,Nu,Nu为未标记数据样本的数量。
10.根据权利要求1中所述的一种基于CELMDAN和SSKFDA的机械状态监测方法,其特征在于,第6步,
Srb=(1-β)Sb+βSt
Srw=(1-β)Sw+βtr(WTW)
其中,总体数据集散度St表示所有样本在空间的分布,β∈[0,1]为正则化因子;当β=1时,SSKFDA退化为KPCA方法;当β=0时,SSKFDA退化为KFDA方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010972737.5A CN112101227A (zh) | 2020-09-16 | 2020-09-16 | 一种基于celmdan和sskfda的机械状态监测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010972737.5A CN112101227A (zh) | 2020-09-16 | 2020-09-16 | 一种基于celmdan和sskfda的机械状态监测方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN112101227A true CN112101227A (zh) | 2020-12-18 |
Family
ID=73759233
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010972737.5A Withdrawn CN112101227A (zh) | 2020-09-16 | 2020-09-16 | 一种基于celmdan和sskfda的机械状态监测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112101227A (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112949524A (zh) * | 2021-03-12 | 2021-06-11 | 中国民用航空飞行学院 | 一种基于经验模态分解与多核学习的发动机故障检测方法 |
CN113378191A (zh) * | 2021-06-01 | 2021-09-10 | 贵州大学 | 一种半诚实模型下基于信息熵的安全多方计算方案 |
CN113484336A (zh) * | 2021-07-06 | 2021-10-08 | 杭州电子科技大学 | 基于表面波和bp神经网络的亚表面裂纹尺寸测量方法 |
CN113610179A (zh) * | 2021-08-17 | 2021-11-05 | 安徽容知日新科技股份有限公司 | 一种设备故障检测分类器训练方法、计算设备及存储介质 |
CN114676736A (zh) * | 2022-05-13 | 2022-06-28 | 浙江工业大学 | 一种滚压工件表面质量预测方法 |
CN116299684A (zh) * | 2023-05-17 | 2023-06-23 | 成都理工大学 | 基于人工神经网络中双模态神经元的新型微震分类方法 |
-
2020
- 2020-09-16 CN CN202010972737.5A patent/CN112101227A/zh not_active Withdrawn
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112949524A (zh) * | 2021-03-12 | 2021-06-11 | 中国民用航空飞行学院 | 一种基于经验模态分解与多核学习的发动机故障检测方法 |
CN112949524B (zh) * | 2021-03-12 | 2022-08-26 | 中国民用航空飞行学院 | 一种基于经验模态分解与多核学习的发动机故障检测方法 |
CN113378191A (zh) * | 2021-06-01 | 2021-09-10 | 贵州大学 | 一种半诚实模型下基于信息熵的安全多方计算方案 |
CN113484336A (zh) * | 2021-07-06 | 2021-10-08 | 杭州电子科技大学 | 基于表面波和bp神经网络的亚表面裂纹尺寸测量方法 |
CN113610179A (zh) * | 2021-08-17 | 2021-11-05 | 安徽容知日新科技股份有限公司 | 一种设备故障检测分类器训练方法、计算设备及存储介质 |
CN114676736A (zh) * | 2022-05-13 | 2022-06-28 | 浙江工业大学 | 一种滚压工件表面质量预测方法 |
CN116299684A (zh) * | 2023-05-17 | 2023-06-23 | 成都理工大学 | 基于人工神经网络中双模态神经元的新型微震分类方法 |
CN116299684B (zh) * | 2023-05-17 | 2023-07-21 | 成都理工大学 | 基于人工神经网络中双模态神经元的新型微震分类方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112101227A (zh) | 一种基于celmdan和sskfda的机械状态监测方法 | |
CN110132598B (zh) | 旋转设备滚动轴承故障噪声诊断算法 | |
Grezmak et al. | Interpretable convolutional neural network through layer-wise relevance propagation for machine fault diagnosis | |
Mao et al. | A new online detection approach for rolling bearing incipient fault via self-adaptive deep feature matching | |
CN113255848B (zh) | 基于大数据学习的水轮机空化声信号辨识方法 | |
Zheng et al. | A review on non-model based diagnosis methodologies for PEM fuel cell stacks and systems | |
Wang et al. | CHMM for tool condition monitoring and remaining useful life prediction | |
Huang et al. | Memory residual regression autoencoder for bearing fault detection | |
CN105738109B (zh) | 基于稀疏表示与集成学习的轴承故障分类诊断方法 | |
Maurya et al. | Condition monitoring of machines using fused features from EMD-based local energy with DNN | |
Neupane et al. | Bearing fault detection using scalogram and switchable normalization-based CNN (SN-CNN) | |
US8630962B2 (en) | Error detection method and its system for early detection of errors in a planar or facilities | |
Liu et al. | A two-stage approach for predicting the remaining useful life of tools using bidirectional long short-term memory | |
CN105760839A (zh) | 基于多特征流形学习与支持向量机的轴承故障诊断方法 | |
CN105275833A (zh) | 一种基于CEEMD-STFT时频信息熵和multi-SVM的离心泵故障诊断方法 | |
Tang et al. | A wind turbine bearing fault diagnosis method based on fused depth features in time–frequency domain | |
Feng et al. | Finite-sensor fault-diagnosis simulation study of gas turbine engine using information entropy and deep belief networks | |
Lu et al. | Bearing fault diagnosis based on clustering and sparse representation in frequency domain | |
CN110701487B (zh) | 一种基于KPCA和Cas-SVDD的多工况管道泄漏检测方法 | |
Liu et al. | Composite multi-scale basic scale Entropy based on CEEMDAN and its application in hydraulic pump fault diagnosis | |
CN113159088B (zh) | 一种基于多特征融合和宽度学习的故障监测与诊断方法 | |
Wang et al. | A hybrid approach for identification of concurrent control chart patterns | |
CN113627539A (zh) | 滚动轴承复合故障诊断方法、系统、存储介质及计算设备 | |
Zheng et al. | Rolling element bearing fault diagnosis based on support vector machine | |
Deng et al. | Primary-auxiliary statistical local kernel principal component analysis and its application to incipient fault detection of nonlinear industrial processes |
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 | ||
WW01 | Invention patent application withdrawn after publication |
Application publication date: 20201218 |
|
WW01 | Invention patent application withdrawn after publication |