模式识别实验报告

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时的正4SKF62051797r/min29.95HS5162.2H107.4H70.6H1024（a）正 幅值幅值幅值幅值-- 样本序号

-- 样本序号 1基于时域特征和线性分类器的简单模式识滚动轴承振动信号进行时域特征提取。有效值与峭度的计算分别为x2x2

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

T Tw1类的输出为y11,w2的输出为y21。但实际的输2所示：计算w2点，将其分为200组样本，每组样本512个样本点。对每组样本计算其有效值和峭度值，从而进行时域特征提取，得到200x2的矩阵，从中选取150组样本作为训练样本，构成150x2的训不同故障有效值概率密度图外圈故障外圈故障内圈故障滚动体故障正常轴承特征样本数特征样本数0 样本有效值1504

50

0 圈故障与正常轴承的测试中有3组样本测试测试准确率为94%；内圈故障与外圈故障的测试中有17组样本测试错误，测试准确率为66%。因此，线性分类器难以在不同故障类型中基于小波包的特征小波包变换通过对振动信号进行小波包分解，可以得到每一频带内振动信号的变化规律，db433别为：节点(3,0)-[0,750]Hz、节点(3,1)-[750,1500]Hz、节点(3,2)-[2250,3000]Hz、节点

53db小波基，将到的轴承振动信号S进行3层分解，从而得到第3层从低频到高频的8个子S

式中，Sn,m表示第n层第m计算各个频带信号的总能量。设小波包分解后第nm个频带的重构信号对应的信号总能量为En,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

1n1

(i,j1,2,,

R的特征值（12p）主成分分析可以得到p个主成分，但是，由于各个主成分的方差是递减的，包含的信息量某个特征值占全部特征值合计的。即贡献率

ppi

贡献率越大，说明该主成分所包含的原始变量的信息越强。主成分个数k85%模糊C均值聚类算模糊C均值聚类算法需要两个参数：一个是聚类数目C，另一个是模糊指数m。一般CC>1m，它是一个控制算法的柔性的参数，如果m过大，则聚类效果会很次。CC*N的一个模糊划分矩阵，这个矩阵表示的是每个模糊C均值聚类算法的步骤还是比较简单的，模糊C均值聚类（FCM，即众所周知的模糊ISODATA，是用隶属度确定每个数据点属于某个聚类的程度的一种聚类算法。FCM把n个向量xi（i=1,2,…,n）c个模糊组，并求每组的聚类中心，使得非相似性指标的价值函数达到最小。FCMC均值聚类（HCM）FCM用模糊划分，使得每个给定数据点用值在0和1间的隶属度来确定其属于各个组的程度与引入模糊划分相适U01间的元素。不过，加上归一化规定，一个数据集的隶1：ccuij1,j1,...,那么，FCM的价值函数（或目标函数）

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

i1

ij这里uij介于0和1之间；ci为模糊组I的聚类中心，dij=||ci-xj||为第I个聚类中心与第j个数据点间的欧几里德距离；且m1,是一个模 指数，这样便在C均值聚类中引入 ..,

,1,...,n

..,

nj1

cc

jj

2

uij

i1

j

这里j，j=1到n，是（10）nuum ci

nunuuij

j12cdij k1dkj 由上述两个必要条件可知，模糊C模糊C均值聚类用下列步骤确定聚类中心ci和隶属度矩阵0、1U，使其满足式（10）用式（13）计算c个聚类中心ci，i=1,…,c用（14）计算新的U2模糊CCU不适U进行硬化处理。在最大隶属度法中，U矩阵中每列的最大隶属度值被被置为1，该列的其他项则被置为0，则得到的Uhard矩阵就可以很直观地看到样本的分基于小波包与模糊C均值聚类的轴承故障模式识C620409680895%以C均值聚类的输入。1010样向聚提提输使用FCM聚类算法对诊断样本集进行故障分类，判断待识别故障属于哪一类已知样向聚提提输向向样

主成分分析降维后的数据分析结8012维累积方差贡献率98.45%2维主成分变换矩阵作为降维后的特征样本进行模糊C均值聚类。1104040组为待识别故障样本。将已知故障样本和待识别故障样本合并在一起，组成故障诊断117所示。从图中可以看2103所示。根据硬化后的结果，相同类型的10，这样就可以很直观地看出轴承属于哪种故障44（内圈）3行，第二类（外圈）1行，第三类（动体）2行，第四类（正常）4行，跟已知故障样于小波包与模糊C均值聚类的故障振动方法的有效性。对于一组未知样本，只需知道隶属度内圈内圈外圈滚动体正常中心0 0 73U46.2.1主成分分析降维的数据分析结7所示。 归一化能量归一化能量归一化能量 01

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)小波包节

01

(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)小波包节

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组。从再看第41~80组样本的最大隶属度值分布情况。第41~80组为待识别样本。各类待识别样本分别有10组，分别为41~50组、51~60组、61~70组、71~8041~50组对应已知样本1~10组中的状态，51~6011~20组中的状态，61~70组对应已知样97.5%。6故障诊断集经过聚类分析后的隶属度矩阵结与模糊C95%以上的主成分特征向量。最后以模糊C均值聚类作为故障分类器，对故障特征进参考文 周川,,刘畅等.基于EMD和模糊C均值聚类的滚动轴承故障诊断[J].理工大学学，2009,34(6):34-黄建鸿.基于小波包分析的滚动轴承故障智能诊断[D]. 张威,肖云魁,司爱威等.基于小波包-模糊C均值聚类算法诊断曲轴轴承故障[J].军事交通学,张文志.基于小波包变换的滚动轴承故障诊断[J].中国机械工程,2013,23(3):295-附线性分类器程序%lsmtest-kmeans-%patternloadIR028-0;%12kHz驱动端内圈故障loadOR021@6-012kHz驱动端外圈故障loadB028-0;%12kHz驱动端滚动体故障loadNORMALBD-0;%12kHz驱动端正常样本xii=X056_DE_time(1:102400);%12kHz驱动端内圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHz驱动端外圈故障OR021@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));%把故障数据个样本点分为组，每组512个数据，下面xo(i,:)=xoo((512*(i-xb(i,:)=xbb((512*(i-xn(i,:)=y((512*(i-p3=zeros(200,1);%峭度值初始化，赋值%计算时域特征，在那个时域统计特征的那个m文件里选一forj=1:200p2(j)=sqrt(sum(xi(j^2N);%均方根值,又称有效值p3(j)=sum(xi(j,:)^4)/N;%峭度q2(j)=sqrt(sum(xo(j^2N均方根值,又称有效值q3(j)=sum(xo(j,:)^4)/N;%峭度

m2(j)=sqrt(sum(xb(j,:^2N均方根值,又称有效值m3(j)=sum(xb(j,:)^4)/N;%峭度n2(j)=sqrt(sum(xn(j^2N);%均方根值,又称有效值n3(j)=sum(xn(j,:)^4)/N;%峭度%%%%%%%%%%%%%%%%%%有效值概率密度图ksdensity(p2概率密度曲线图holdonksdensity(q2概率密度曲线图holdonholdonxlabel('样本有效值ylabel('特征样本数title('不同故障有效值概率密度图%线性分类器x1=[p2(1:150),p3(1:150)];%选前150组样本训练fori=1:150forr3(j)=x2(j,1);holdon%初始化函y1ones(150,1w1类的期望输出为-1y2=ones(150,1);%w2类的期望输出为1x1(:,31;%考虑到不经过原点的超平面，对x进行扩x2(:,31;%使x'=[x1]，x为2维的，故加1扩为3维x=[x1;x2]';%使x矩阵化yy1;y2];%使y矩阵R= %求出自相关矩E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个x的点作yw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面title('正常轴承与内圈故障轴承分类结果for

disp('正常轴承fori=1:150forholdon%初始化函y3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,31;%考虑到不经过原点的超平面，对x进行扩x2(:,31;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化yy3;y4];%使y矩阵R= %求出自相关矩E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个x的点作yw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面title('正常轴承与外圈故障轴承分类结果for

disp('正常轴承fori=1:150forholdon%初始化函y3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,31;%考虑到不经过原点的超平面，对x进行扩x2(:,31;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化yy3;y4];%使y矩阵R= %求出自相关矩E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个x的点作yw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面%测试程for

disp('滚动体故障fori=1:150for

holdon%初始化函y3ones(150,1w1类的期望输出为-1y4=ones(150,1);%w2类的期望输出为1x3(:,31;%考虑到不经过原点的超平面，对x进行扩x2(:,31;%使x'=[x1]，x为2维的，故加1扩为3维x=[x3;x2]';%使x矩阵化yy3;y4];%使y矩阵R= %求出自相关矩E=x*y; w=inv(R)*E;%求权向量估计值xlinspace(0,1,5000);%取5000个x的点作yw(1)/w(2))*x-w(3)/w(2);%x*w1+y*w2+w0=0,w=[w1;w2;w0]plot(x,y,'r');%用红线画出分界面%测试程for

disp('外圈故障基于小波包与模糊C均值聚类的程loadIR028-0;%12kHz驱动端内圈故障loadOR021@6-012kHz驱动端外圈故障loadB028-0;%12kHz驱动端滚动体故障loadNORMALBD-0;%12kHz驱动端正常样本xii=X056_DE_time(1:102400);%12kHz驱动端内圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHz驱动端外圈故障OR021@6-0xbb=X048_DE_time(1:102400);%12kHz驱动端滚动体故障B028-0y=X097_DE_time(1:102400正常轴承fori=1:20xi(i,:)=xii((5120*(i-xo(i,:)=xoo((5120*(i-xb(i,:)=xbb((5120*(i-xn(i,:)=y((5120*(i-forxi(i,:)=xi(i,:)-wpt0=wpdec(xi(i,:),3,'db4对数据进行小波包分解forj=1:8%对各层小波包分解系数进行重构E0(i,j)=norm(wprcoef(wpt0,[3,j-1]),2)*norm(wprcoef(wpt0,[3,j-E0_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(i1,j1)=norm(wprcoef(wpt1,[3,j1-1]),2)*norm(wprcoef(wpt1,[3,j1-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(i2,j2)=norm(wprcoef(wpt2,[3,j2-1]),2)*norm(wprcoef(wpt2,[3,j2-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(i3,j3)=norm(wprcoef(wpt3,[3,j3-1]),2)*norm(wprcoef(wpt3,[3,j3-E3_total(i3)=sum(E3(i3求小波包分解总能量E_totalforj3=1:8pn(i3,j3)=E3(i3,j3)/E3_total(i3);%求归一化后小波包各频带的能量分 %coeff为变换空间中的基向%score为y1在变换空间中的主成分表latent为协方差矩阵的本征值con_rate=cumsum(latentsum(latent);[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)

%coeff为变换空间中的基向%score为y1在变换空间中的主成分表latent为协方差矩阵的本征值con_rate2=cumsum(latent2sum(latent2);[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,% 采用模糊C均值对数据集data聚为cluster_n%输入%----nxm矩阵,表示n个样本,每个样本具有m维的特征%----标量,表示聚合中心数目,即类别%4x1矩阵，其 options(1):隶属度矩阵U的指数 (缺省值: options(2):最大迭代次 (缺省值: options(3):隶属度最小变化量,迭代终止条 (缺省值:1e- options(4):每次迭代是否输出信息标 (缺省值:%输出%类中%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(data1出data的第一维(rows)数,即样本个in_nsize(data %求出data的第二维(columns)数，即特征值长认操作参default_options 属度矩阵U的指 大迭代次1e- %隶属度最小变化量,迭代终止条 次迭代是否输出信息标ifnargin==options= %分析有options做参数时候的情%如果输入参数个数是二那么就调用默认的iflength(options4如果用户给的opition数少于4个那么其他用默认值;tmp=default_options;tmp(1:length(options))=options;options=tmp;回options中是数的值为0(如NaN),不是数时为1nan_index=find(isnan(options)==1);options(nan_index)=default_options(nan_index);ifoptions(1)<=1,%如果模糊矩阵的指数小于等于1error('Theexponen