模式識別實驗報告

C(徐清瑤,儀器科學與技術(shù)：針對滾動軸承振動信號的非平穩(wěn)特性，本文采用基于小波包與模糊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)剛好解決了“故障漸變的模糊性”問題。為此，本文將小波包變換與主成分分析相結(jié)合提取的特征向量作為模糊C均值聚類的輸入，來識別數(shù)據(jù)來故障數(shù)據(jù)來源于凱斯西儲大學軸承數(shù)據(jù)中心本文選取的是采樣頻率為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組樣本作為訓練樣本，構(gòu)成150x2的訓不同故障有效值概率密度圖外圈故障外圈故障內(nèi)圈故障滾動體故障正常軸承特征樣本數(shù)特征樣本數(shù)0 樣本有效值1504

50

0 圈故障與正常軸承的測試中有3組樣本測試測試準確率為94%；內(nèi)圈故障與外圈故障的測試中有17組樣本測試錯誤，測試準確率為66%。因此，線性分類器難以在不同故障類型中基于小波包的特征小波包變換通過對振動信號進行小波包分解，可以得到每一頻帶內(nèi)振動信號的變化規(guī)律，db433別為：節(jié)點(3,0)-[0,750]Hz、節(jié)點(3,1)-[750,1500]Hz、節(jié)點(3,2)-[2250,3000]Hz、節(jié)點

53db小波基，將到的軸承振動信號S進行3層分解，從而得到第3層從低頻到高頻的8個子S

式中，Sn,m表示第n層第m計算各個頻帶信號的總能量。設(shè)小波包分解后第nm個頻帶的重構(gòu)信號對應(yīng)的信號總能量為En,m，則

NN2|2

n,i

式中，n為小波包分解層數(shù)，m表示分解節(jié)點的序號，N表示每個節(jié)點數(shù)據(jù)的長度，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

貢獻率越大，說明該主成分所包含的原始變量的信息越強。主成分個數(shù)k85%模糊C均值聚類算模糊C均值聚類算法需要兩個參數(shù)：一個是聚類數(shù)目C，另一個是模糊指數(shù)m。一般CC>1m，它是一個控制算法的柔性的參數(shù)，如果m過大，則聚類效果會很次。CC*N的一個模糊劃分矩陣，這個矩陣表示的是每個模糊C均值聚類算法的步驟還是比較簡單的，模糊C均值聚類（FCM，即眾所周知的模糊ISODATA，是用隸屬度確定每個數(shù)據(jù)點屬于某個聚類的程度的一種聚類算法。FCM把n個向量xi（i=1,2,…,n）c個模糊組，并求每組的聚類中心，使得非相似性指標的價值函數(shù)達到最小。FCMC均值聚類（HCM）FCM用模糊劃分，使得每個給定數(shù)據(jù)點用值在0和1間的隸屬度來確定其屬于各個組的程度與引入模糊劃分相適U01間的元素。不過，加上歸一化規(guī)定，一個數(shù)據(jù)集的隸1：ccuij1,j1,...,那么，F(xiàn)CM的價值函數(shù)（或目標函數(shù)）

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

i1

ij這里uij介于0和1之間；ci為模糊組I的聚類中心，dij=||ci-xj||為第I個聚類中心與第j個數(shù)據(jù)點間的歐幾里德距離；且m1,是一個模 指數(shù)，這樣便在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聚類算法對診斷樣本集進行故障分類，判斷待識別故障屬于哪一類已知樣向聚提提輸向向樣

主成分分析降維后的數(shù)據(jù)分析結(jié)8012維累積方差貢獻率98.45%2維主成分變換矩陣作為降維后的特征樣本進行模糊C均值聚類。1104040組為待識別故障樣本。將已知故障樣本和待識別故障樣本合并在一起，組成故障診斷117所示。從圖中可以看2103所示。根據(jù)硬化后的結(jié)果，相同類型的10，這樣就可以很直觀地看出軸承屬于哪種故障44（內(nèi)圈）3行，第二類（外圈）1行，第三類（動體）2行，第四類（正常）4行，跟已知故障樣于小波包與模糊C均值聚類的故障振動方法的有效性。對于一組未知樣本，只需知道隸屬度內(nèi)圈內(nèi)圈外圈滾動體正常中心0 0 73U46.2.1主成分分析降維的數(shù)據(jù)分析結(jié)7所示。 歸一化能量歸一化能量歸一化能量 01

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)小波包節(jié)

01

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)小波包節(jié) 歸一化能量歸一化能量歸一化能量 0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)小波包節(jié)

0

(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)小波包節(jié)7從圖7中可以看出，正常狀態(tài)下軸承振動信號的能量主要分布在節(jié)點(3,0)和(3,1)處，即0~750Hz750~1500Hz的頻段范圍內(nèi)，這是由周期性振源引起的響應(yīng)。當軸承出現(xiàn)故障時，2250~3000Hz3000~3750Hz的頻段范圍內(nèi)，這是由于軸承的振動信號具有明顯的調(diào)制特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組對應(yīng)已知樣本1~10組中的狀態(tài)，51~6011~20組中的狀態(tài)，61~70組對應(yīng)已知樣97.5%。6故障診斷集經(jīng)過聚類分析后的隸屬度矩陣結(jié)與模糊C95%以上的主成分特征向量。最后以模糊C均值聚類作為故障分類器，對故障特征進參考文 周川,,劉暢等.基于EMD和模糊C均值聚類的滾動軸承故障診斷[J].理工大學學，2009,34(6):34-黃建鴻.基于小波包分析的滾動軸承故障智能診斷[D]. 張威,肖云魁,司愛威等.基于小波包-模糊C均值聚類算法診斷曲軸軸承故障[J].軍事交通學,張文志.基于小波包變換的滾動軸承故障診斷[J].中國機械工程,2013,23(3):295-附線性分類器程序%lsmtest-kmeans-%patternloadIR028-0;%12kHz驅(qū)動端內(nèi)圈故障loadOR021@6-012kHz驅(qū)動端外圈故障loadB028-0;%12kHz驅(qū)動端滾動體故障loadNORMALBD-0;%12kHz驅(qū)動端正常樣本xii=X056_DE_time(1:102400);%12kHz驅(qū)動端內(nèi)圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHz驅(qū)動端外圈故障OR021@6-0xbb=X048_DE_time(1:102400);%12kHz驅(qū)動端滾動體故障B028-0y=X097_DE_time(1:102400正常軸承xo=zeros(200,512);%即開始賦值0，后面再把特征數(shù)據(jù)裝進去fori=1:200xi(i,:)=xii((512*(i-1)+1:512*i));%把故障數(shù)據(jù)個樣本點分為組，每組512個數(shù)據(jù)，下面xo(i,:)=xoo((512*(i-xb(i,:)=xbb((512*(i-xn(i,:)=y((512*(i-p3=zeros(200,1);%峭度值初始化，賦值%計算時域特征，在那個時域統(tǒng)計特征的那個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('特征樣本數(shù)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;%考慮到不經(jīng)過原點的超平面，對x進行擴x2(:,31;%使x'=[x1]，x為2維的，故加1擴為3維x=[x1;x2]';%使x矩陣化yy1;y2];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計值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('正常軸承與內(nèi)圈故障軸承分類結(jié)果for

disp('正常軸承fori=1:150forholdon%初始化函y3ones(150,1w1類的期望輸出為-1y4=ones(150,1);%w2類的期望輸出為1x3(:,31;%考慮到不經(jīng)過原點的超平面，對x進行擴x2(:,31;%使x'=[x1]，x為2維的，故加1擴為3維x=[x3;x2]';%使x矩陣化yy3;y4];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計值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('正常軸承與外圈故障軸承分類結(jié)果for

disp('正常軸承fori=1:150forholdon%初始化函y3ones(150,1w1類的期望輸出為-1y4=ones(150,1);%w2類的期望輸出為1x3(:,31;%考慮到不經(jīng)過原點的超平面，對x進行擴x2(:,31;%使x'=[x1]，x為2維的，故加1擴為3維x=[x3;x2]';%使x矩陣化yy3;y4];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計值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;%考慮到不經(jīng)過原點的超平面，對x進行擴x2(:,31;%使x'=[x1]，x為2維的，故加1擴為3維x=[x3;x2]';%使x矩陣化yy3;y4];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計值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驅(qū)動端內(nèi)圈故障loadOR021@6-012kHz驅(qū)動端外圈故障loadB028-0;%12kHz驅(qū)動端滾動體故障loadNORMALBD-0;%12kHz驅(qū)動端正常樣本xii=X056_DE_time(1:102400);%12kHz驅(qū)動端內(nèi)圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHz驅(qū)動端外圈故障OR021@6-0xbb=X048_DE_time(1:102400);%12kHz驅(qū)動端滾動體故障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對數(shù)據(jù)進行小波包分解forj=1:8%對各層小波包分解系數(shù)進行重構(gòu)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對數(shù)據(jù)進行小波包分解forj1=1:8%對各層小波包分解系數(shù)進行重構(gòu)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對數(shù)據(jù)進行小波包分解forj2=1:8%對各層小波包分解系數(shù)進行重構(gòu)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對數(shù)據(jù)進行小波包分解forj3=1:8%對各層小波包分解系數(shù)進行重構(gòu)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為協(xié)方差矩陣的本征值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為協(xié)方差矩陣的本征值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均值對數(shù)據(jù)集data聚為cluster_n%輸入%----nxm矩陣,表示n個樣本,每個樣本具有m維的特征%----標量,表示聚合中心數(shù)目,即類別%4x1矩陣，其 options(1):隸屬度矩陣U的指數(shù) (缺省值: options(2):最大迭代次 (缺省值: options(3):隸屬度最小變化量,迭代終止條 (缺省值:1e- options(4):每次迭代是否輸出信息標 (缺省值:%輸出%類中%U屬度矩%標函數(shù)% 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)數(shù),即樣本個in_nsize(data %求出data的第二維(columns)數(shù)，即特征值長認操作參default_options 屬度矩陣U的指 大迭代次1e- %隸屬度最小變化量,迭代終止條 次迭代是否輸出信息標ifnargin==options= %分析有options做參數(shù)時候的情%如果輸入?yún)?shù)個數(shù)是二那么就調(diào)用默認的iflength(options4如果用戶給的opition數(shù)少于4個那么其他用默認值;tmp=default_options;tmp(1:length(options))=options;options=tmp;回options中是數(shù)的值為0(如NaN),不是數(shù)時為1nan_index=find(isnan(options)==1);options(nan_index)=default_options(nan_index);ifoptions(1)<=1,%如果模糊矩陣的指數(shù)小于等于1error('Theexponen