模式识别实验报告

基于包与模糊C均值聚类的滚动轴承故障诊(,仪器科学与技术：针对滚动轴承振动信号的非平稳特性，本文采用基于包与模糊C均值聚类的模式识别方法法对该特征向量进行降维，将降维后的主成分矩阵作为故障特征向量，最后以模糊C：滚动轴承；包；模糊C均值聚类；主成分分RollingBearingfaultdiagnosisbasedonwaveletpacketandFuzzyC(Qing-yaoXu,InstrumentScienceand:Aimingatthenon-stationaryfeatureofavibrationsingalofrollingbearings,apatternrecognitionmethodbasedonwaveletpacket(WPT)andFuzzyCMeansClustering(FCM)wasusedtodiagnosethefaultsofrollingbearings.Firstly,thearticleysedthedisadvantagesofsimplepatternrecognitionmethodbasedontime-characteristicsandlinearclassifier.Secondly,waveletpacket positionofvibrationsignalsofbearingswasdone,andthewaveletpacketnodenormalizedenergyvaluewasextractedasthefeaturevectors.Thenthemethodofprincipalcomponentysis(PCA)wasusedtoreducethedimensionsofthefeaturevectors,andthecomponentmatrixwasusedasfaultfeaturevectors.Finally,FCMwasusedasfaultclassifiertorealizetheidentificationofrollingbearingfaulttypes.Resultsshowedthatthismethodcaneffectivelycarryoutthefaultdiagnosisofrolling:Rollingbearings;Waveletpacket;FuzzyCMeansClustering;Principalcomponent，C均值聚类(FuzzyCMeansClustering,FCM)刚好解决了“故障渐变的模糊性”问题。为此，本文将包变换与主成分分析相结合提取的特征向量作为模糊C均值聚类的输入，来识别故障数据来源于西储大学轴承数据中心本文选取的是采样频率为12kHz时的正4SKF62051024（a）正 -- 样本序号

-- 样本序号 1滚动轴承振动信号进行时域特征提取。有效值与峭度的计算分别为x2x2

4p(x)dx1T4(t)dt1T()4

T Tw1y11,w2y21。但实际的输2所示：2不同故障有效值概率密度图外圈故障外圈故障内圈故障滚动体故障正常轴承特征样本数特征样本数0 样本有效值1504

50

0 基于包的特征提包变换通过对振动信号进行包分解，可以得到每一频带内振动信号的变化规律，包的每次分解都是利用高、低通滤波器将信号分解到高、低两个频带范围内，相比于分解包分解在保留信号时频局部化特性的基础之上在时频连续区提取暂态信号的特征，本文采用db4基对振动信号进行3层包分解，分别提取第3层从低频到高频的，个频率成分信号包分解的树结构图如图5所示。第3层中各节点重构信号的频带范围，别为：节点(3,0)-[0,750]Hz、节点(3,1)-[750,1500]Hz、节点(3,2)-[2250,3000]Hz、节点基于包的特征提取步

图5包分解树结构本文在用包对原始振动信号进行特征提取时，分解层数取为3层，选用db基，将到的轴承振动信号S进行3层分解，从而得到第3层从低频到高频的8个子S

式中，Sn,mnm计算各个频带信号的总能量。设包分解后第n层第m个频带的重构信号

NN2|2

n,i

式中，n为包分解层数，m表示分解节点的序号，N表示每个节点数据的长度，xn,iSn,mi个离散点的幅值。求取包的能量分布。包3层分解各频带的总能量E3等于各频带的能量之77E3则归一化 包各频带的能量分布表示如

TE3,0,E3,1,E3,2,E3,3,E3,4,E3,5,E3,6,E3,7

E 388080x8

r1p rR

2p

r p ppX nrij xtixtjtn

(i,j1,2,,

某个特征值占全部特征值合计的。即贡献率

ppi

贡献率越大，说明该主成分所包含的原始变量的信息越强。主成分个数k的选取，主要根CC>1m，它是一个控制算法的柔性m过大，则聚类效果会很次。CC*N的一个模糊划分矩阵，这个矩阵表示的是每个（FCMISODATA，是用隶属度确定每个数据点属于某个聚类的程度的一种聚类算法。标的价值函数达到最小。FCMC均值聚类（HCM）FCM用模糊划分，使U01间的元素。不过，加上归一化规定，一个数据集的隶1：ccuij1,j1,...,那么，FCM的价值函数（或目标函数）

J(U,c,...,c)Jumd

i1

ij这里uij介于0和1之间；ci为模糊组I的聚类中心，dij=||ci-xj||为第I个聚类中心与第j个 指数，这样便在C均值聚类中引入

)J(U,c1,...,

)j)

cc

umd2

ij i1 uum ci

nunuuij

j12cdij k1dkj CCci0、1U，使其满足式（10）用式（13）cci，i=1,…,c用（14）U2CCU不适U进行硬化处理。在最大隶属度法中，U矩阵中每列的最大隶属度基于包与模糊C均值聚类的轴承故障模式识别流程如图6所示。具体诊断步骤如下204096对80组样本中的每组都进行包分解得到一系列具有不同特征尺度的平稳信号，求最后一层各节点的包归一化能量值作为8维初始特征向量。95%以C均值聚类的输入。1010

图6基于包与模糊C均值聚类的轴承故障模式识12维累积方差贡献率98.45%2C均值聚类。1104040组为待识别故障样本。将已知故障样本和待识别故障样本合并在一起，组成故障诊断117所示。从图中可以看2103所示。根据硬化后的结果，相同类型的44（内圈）3行，第二类（外圈）1行，第三类（动体）2行，第四类（正常）4行，跟已知故障样内圈内圈外圈滚动体正常中心0 0 723U456.2.1主成分分析降维的数据分析结当振动信号过主成分分析降维时，可通过对振动信号经过特征提取后的各频段包能量分布情况进行分析来提取能量分布集中的包归一化能量值各频段包能量分布情7所示。 0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)

0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6) (a)正常轴 0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)

0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)7从图7中可以看出，正常状态下轴承振动信号的能量主要分布在节点(3,0)和(3,1)处，即0~750Hz750~1500Hz的频段范围内，这是由周期性振源引起的响应。当轴承出现故障时，2250~3000Hz3000~3750Hz的频段范围内，这是由于轴承的振动信号具有明显的调制特析之前取这四个节点的包归一化能量值作为特征直接输入FCM分类器中。61~40组为已101~10组、11~20组、21~30组、31~40组。从1~10组中的状态，51~6011~20组中的状态，61~70组对应已知样97.5%。6布集中频段的包能量特征虽然也能达到很高的分类正确率但却不能保证故障全部分类正针对滚动轴承故障振动信号的非平稳特征和故障征兆的模糊性，本文提出了基于包与模糊C均值聚类的滚动轴承故障诊断方法。该方法首先对滚动轴承的振动信号进行 周川,伍星,刘畅等.基于EMD和模糊C均值聚类的滚动轴承故障诊断[J].理工大学学，2009,34(6):34-黄建鸿.基于包分析的滚动轴承故障智能诊断[D].南昌：南昌大学 张威,,司爱威等.基于包-模糊C均值聚类算法诊断曲轴轴承故障[J].军事交通学王冬云,.基于包变换的滚动轴承故障诊断[J].中国机械工程,2013,23(3):295-线性分类器程序%patternloadIR028-0;%12kHz驱动端内圈故障loadOR021@6-0;%12kHz驱动端外圈故障loadB028-0;%12kHz驱动端滚动体故障loadNORMALBD-0;%12kHz驱动端正常样本xii=X056_DE_time(1:102400);%12kHz驱动端内圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHzOR021@6-0xbb=X048_DE_time(1:102400);%12kHz驱动端滚动体故障B028-0y=X097_DE_time(1:102400正常轴承xo=zeros(200,512即开始赋值0，后面再把特征数据装进去fori=1:200xi(i,:)=xii((512*(i-1)+1:512*i));%把故障数据个样本点分为组，每组512forj=1:200p2(j)=sqrt(sum(xi(j,:).^2)./N均方根值,又称有效值%q2(j)=sqrt(sum(xo(j,:).^2)./N均方根值,又称有效值q3(j)=sum(xo(j,:).^4)./N;

m2(j)=sqrt(sum(xb(j,:).^2)./N均方根值,又称有效值m3(j)=sum(xb(j,:).^4)./N;n2(j)=sqrt(sum(xn(j,:).^2)./N均方根值,又称有效值n3(j)=sum(xn(j,:).^4)./N;ksdensity(p2概率密度曲线图holdonksdensity(q2概率密度曲线图holdonholdonx1=[p2(1:150),p3(1:150)];%选前150组样本训练fori=1:150forholdony1ones(150,1w1类的期望输出为-1y2=ones(150,1);%w2类的期望输出为1x1(:,3)1;%考虑到不经过原点的超平面，对xx2(:,3)1;%使x'=[x1]，x为2维的，故加1扩为3维x=[x1;x2]';%使x矩阵化y[y1;y2];%使yR %E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个xyw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面for

fori=1:150forholdony3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,3)1;%考虑到不经过原点的超平面，对xx2(:,3)1;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化y[y3;y4];%使yR %E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个xyw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面for

fori=1:150forholdony3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,3)1;%考虑到不经过原点的超平面，对xx2(:,3)1;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化y[y3;y4];%使yR %E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个xyw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面for

fori=1:150for

holdony3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,3)1;%考虑到不经过原点的超平面，对xx2(:,3)1;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化y[y3;y4];%使yR %E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个xyw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面for

loadIR028-0;%12kHz驱动端内圈故障loadOR021@6-0;%12kHz驱动端外圈故障loadB028-0;%12kHz驱动端滚动体故障loadNORMALBD-0;%12kHz驱动端正常样本xii=X056_DE_time(1:102400);%12kHz驱动端内圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHzOR021@6-0xbb=X048_DE_time(1:102400);%12kHz驱动端滚动体故障B028-0y=X097_DE_time(1:102400);正常轴承fori=1:20forE0_total(i)=sum(E0(i,:));%求包分解总能量E_totalforj=1:8pi(i,j)=E0(i,j)/E0_total(i);%求归一化后包各频带的能量分forwpt1=wpdec(xo(i1,:),3,'db4');%对数据进行包分解forj1=1:8%对各层包分解系数进行重构E1_total(i1)=sum(E1(i1,:));%求包分解总能量E_totalforj1=1:8po(i1,j1)=E1(i1,j1)/E1_total(i1);%求归一化后包各频带的能量分forwpt2=wpdec(xb(i2,:),3,'db4');%对数据进行包分解forj2=1:8%对各层包分解系数进行重构E2_total(i2)=sum(E2(i2,:));%求包分解总能量E_totalforj2=1:8pb(i2,j2)=E2(i2,j2)/E2_total(i2);%求归一化后包各频带的能量分forwpt3=wpdec(xn(i3,:),3,'db4');%对数据进行包分解forj3=1:8%对各层包分解系数进行重构E3_total(i3)=sum(E3(i3,:));%求包分解总能量E_totalforj3=1:8pn(i3,j3)=E3(i3,j3)/E3_total(i3);%求归一化后包各频带的能量分 %score为y1latent为协方差矩阵的本征值[center,U, ]=fcm_xqy(pc,4);maxU=%Findthedatapointswithhighestgradeofmembershipincluster1index1=find(U(1,:)==maxU);%Findthedatapointswithhighestgradeofmembershipincluster2index2=find(U(2,:)==maxU);index3=find(U(3,:)==maxU);index4=find(U(4,:)==holdonholdonholdon%Plottheclustercentersholdonholdonholdonholdonfori=1:40forj=1:4ifUhard(j,i)==Uhard(b,i)

%score为y1latent为协方差矩阵的本征值[center2,U2, 2]=fcm_xqy(pc2,4);fori=1:80forj=1:4ifUhard2(j,i)==Uhard2(b,i)

模糊C均值聚类程function[center, ]=fcm_xqy(data,cluster_n, %%%nxm矩阵,表示n个样本,每个样本具有m4x1%%%1e-%%%U%% data= ]= plot(data(:,1), hold maxU= index1=find(U(1,:)== index2=find(U(2,:)== plot([center([12],1)],[center([1 holdifnargin~=2&nargin~=3, error('Toomanyortoofewinputarguments!');data_nsize(data1data的第一维(rowsin_nsize(data data的第二维(columnsdefault_options U ifnargin==options= iflength(options)4opition数少于4个那么其他用默认值;tmp=default_options;tmp(1:length(options))=options;options=tmp;options中是数的值为0(如NaN),不是数时为1nan_index=find(isnan(options)==