模式識(shí)別實(shí)驗(yàn)報(bào)告

C(徐清瑤,儀器科學(xué)與技術(shù)：針對(duì)滾動(dòng)軸承振動(dòng)信號(hào)的非平穩(wěn)特性，本文采用基于小波包與模糊C均值聚類的模式識(shí)別方法法對(duì)該特征向量進(jìn)行降維，將降維后的主成分矩陣作為故障特征向量，最后以模糊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í)別數(shù)據(jù)來故障數(shù)據(jù)來源于凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心本文選取的是采樣頻率為12kHz時(shí)的正4SKF62051797r/min29.95HS5162.2H107.4H70.6H1024（a）正 幅值幅值幅值幅值-- 樣本序號(hào)

-- 樣本序號(hào) 1基于時(shí)域特征和線性分類器的簡單模式識(shí)滾動(dòng)軸承振動(dòng)信號(hào)進(jìn)行時(shí)域特征提取。有效值與峭度的計(jì)算分別為x2x2

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

T Tw1類的輸出為y11,w2的輸出為y21。但實(shí)際的輸2所示：計(jì)算w2點(diǎn)，將其分為200組樣本，每組樣本512個(gè)樣本點(diǎn)。對(duì)每組樣本計(jì)算其有效值和峭度值，從而進(jìn)行時(shí)域特征提取，得到200x2的矩陣，從中選取150組樣本作為訓(xùn)練樣本，構(gòu)成150x2的訓(xùn)不同故障有效值概率密度圖外圈故障外圈故障內(nèi)圈故障滾動(dòng)體故障正常軸承特征樣本數(shù)特征樣本數(shù)0 樣本有效值1504

50

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

53db小波基，將到的軸承振動(dòng)信號(hào)S進(jìn)行3層分解，從而得到第3層從低頻到高頻的8個(gè)子S

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

NN2|2

n,i

式中，n為小波包分解層數(shù)，m表示分解節(jié)點(diǎn)的序號(hào)，N表示每個(gè)節(jié)點(diǎn)數(shù)據(jù)的長度，xn,iSn,mi個(gè)離散點(diǎn)的幅值。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個(gè)主成分，但是，由于各個(gè)主成分的方差是遞減的，包含的信息量某個(gè)特征值占全部特征值合計(jì)的。即貢獻(xiàn)率

ppi

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

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

i1

ij這里uij介于0和1之間；ci為模糊組I的聚類中心，dij=||ci-xj||為第I個(gè)聚類中心與第j個(gè)數(shù)據(jù)點(diǎn)間的歐幾里德距離；且m1,是一個(gè)模 指數(shù)，這樣便在C均值聚類中引入 ..,

,1,...,n

..,

nj1

cc

jj

2

uij

i1

j

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

nunuuij

j12cdij k1dkj 由上述兩個(gè)必要條件可知，模糊C模糊C均值聚類用下列步驟確定聚類中心ci和隸屬度矩陣0、1U，使其滿足式（10）用式（13）計(jì)算c個(gè)聚類中心ci，i=1,…,c用（14）計(jì)算新的U2模糊CCU不適U進(jìn)行硬化處理。在最大隸屬度法中，U矩陣中每列的最大隸屬度值被被置為1，該列的其他項(xiàng)則被置為0，則得到的Uhard矩陣就可以很直觀地看到樣本的分基于小波包與模糊C均值聚類的軸承故障模式識(shí)C620409680895%以C均值聚類的輸入。1010樣向聚提提輸使用FCM聚類算法對(duì)診斷樣本集進(jìn)行故障分類，判斷待識(shí)別故障屬于哪一類已知樣向聚提提輸向向樣

主成分分析降維后的數(shù)據(jù)分析結(jié)8012維累積方差貢獻(xiàn)率98.45%2維主成分變換矩陣作為降維后的特征樣本進(jìn)行模糊C均值聚類。1104040組為待識(shí)別故障樣本。將已知故障樣本和待識(shí)別故障樣本合并在一起，組成故障診斷117所示。從圖中可以看2103所示。根據(jù)硬化后的結(jié)果，相同類型的10，這樣就可以很直觀地看出軸承屬于哪種故障44（內(nèi)圈）3行，第二類（外圈）1行，第三類（動(dòng)體）2行，第四類（正常）4行，跟已知故障樣于小波包與模糊C均值聚類的故障振動(dòng)方法的有效性。對(duì)于一組未知樣本，只需知道隸屬度內(nèi)圈內(nèi)圈外圈滾動(dòng)體正常中心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)下軸承振動(dòng)信號(hào)的能量主要分布在節(jié)點(diǎn)(3,0)和(3,1)處，即0~750Hz750~1500Hz的頻段范圍內(nèi)，這是由周期性振源引起的響應(yīng)。當(dāng)軸承出現(xiàn)故障時(shí)，2250~3000Hz3000~3750Hz的頻段范圍內(nèi)，這是由于軸承的振動(dòng)信號(hào)具有明顯的調(diào)制特FCM分類器中。61~40組為已101~10組、11~20組、21~30組、31~40組。從再看第41~80組樣本的最大隸屬度值分布情況。第41~80組為待識(shí)別樣本。各類待識(shí)別樣本分別有10組，分別為41~50組、51~60組、61~70組、71~8041~50組對(duì)應(yīng)已知樣本1~10組中的狀態(tài)，51~6011~20組中的狀態(tài)，61~70組對(duì)應(yīng)已知樣97.5%。6故障診斷集經(jīng)過聚類分析后的隸屬度矩陣結(jié)與模糊C95%以上的主成分特征向量。最后以模糊C均值聚類作為故障分類器，對(duì)故障特征進(jìn)參考文 周川,,劉暢等.基于EMD和模糊C均值聚類的滾動(dòng)軸承故障診斷[J].理工大學(xué)學(xué)，2009,34(6):34-黃建鴻.基于小波包分析的滾動(dòng)軸承故障智能診斷[D]. 張威,肖云魁,司愛威等.基于小波包-模糊C均值聚類算法診斷曲軸軸承故障[J].軍事交通學(xué),張文志.基于小波包變換的滾動(dòng)軸承故障診斷[J].中國機(jī)械工程,2013,23(3):295-附線性分類器程序%lsmtest-kmeans-%patternloadIR028-0;%12kHz驅(qū)動(dòng)端內(nèi)圈故障loadOR021@6-012kHz驅(qū)動(dòng)端外圈故障loadB028-0;%12kHz驅(qū)動(dòng)端滾動(dòng)體故障loadNORMALBD-0;%12kHz驅(qū)動(dòng)端正常樣本xii=X056_DE_time(1:102400);%12kHz驅(qū)動(dòng)端內(nèi)圈故障IR028-0xoo=X234_DE_time(1:102400);%12kHz驅(qū)動(dòng)端外圈故障OR021@6-0xbb=X048_DE_time(1:102400);%12kHz驅(qū)動(dòng)端滾動(dòng)體故障B028-0y=X097_DE_time(1:102400正常軸承xo=zeros(200,512);%即開始賦值0，后面再把特征數(shù)據(jù)裝進(jìn)去fori=1:200xi(i,:)=xii((512*(i-1)+1:512*i));%把故障數(shù)據(jù)個(gè)樣本點(diǎn)分為組，每組512個(gè)數(shù)據(jù)，下面xo(i,:)=xoo((512*(i-xb(i,:)=xbb((512*(i-xn(i,:)=y((512*(i-p3=zeros(200,1);%峭度值初始化，賦值%計(jì)算時(shí)域特征，在那個(gè)時(shí)域統(tǒng)計(jì)特征的那個(gè)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組樣本訓(xùn)練fori=1:150forr3(j)=x2(j,1);holdon%初始化函y1ones(150,1w1類的期望輸出為-1y2=ones(150,1);%w2類的期望輸出為1x1(:,31;%考慮到不經(jīng)過原點(diǎn)的超平面，對(duì)x進(jìn)行擴(kuò)x2(:,31;%使x'=[x1]，x為2維的，故加1擴(kuò)為3維x=[x1;x2]';%使x矩陣化yy1;y2];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計(jì)值xlinspace(0,1,5000);%取5000個(gè)x的點(diǎn)作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)過原點(diǎn)的超平面，對(duì)x進(jìn)行擴(kuò)x2(:,31;%使x'=[x1]，x為2維的，故加1擴(kuò)為3維x=[x3;x2]';%使x矩陣化yy3;y4];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計(jì)值xlinspace(0,1,5000);%取5000個(gè)x的點(diǎn)作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)過原點(diǎn)的超平面，對(duì)x進(jìn)行擴(kuò)x2(:,31;%使x'=[x1]，x為2維的，故加1擴(kuò)為3維x=[x3;x2]';%使x矩陣化yy3;y4];%使y矩陣R= %求出自相關(guān)矩E=x*y; w=inv(R)*E;%求權(quán)向量估計(jì)值xlinspace(0,1,5000);%取5000個(gè)x的點(diǎn)作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('滾動(dòng)體故障fori=1:150for

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