基于均方頻率與EMD的切削顫振特征提取方法*

邢諾貝 劉福軍 周 超 許黎明

(上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

顫振是金屬加工中由于刀具與工件之間的自激效應(yīng)而自發(fā)產(chǎn)生的一種振動(dòng)現(xiàn)象。顫振對(duì)工件的表面質(zhì)量、刀具的壽命、乃至機(jī)床的使用壽命都有很大的影響。切削顫振的產(chǎn)生機(jī)制于1946年以來(lái)已經(jīng)被眾多學(xué)者作為研究的熱點(diǎn)，研究發(fā)現(xiàn)切削顫振主要由再生[1]、摩擦[2]以及模態(tài)耦合[3]的影響而引起，其中再生效應(yīng)被廣泛用于解釋加工過(guò)程中的顫振機(jī)理。為了提高機(jī)床加工的穩(wěn)定性，穩(wěn)定性Lobe圖[4-5]于1964年被提出，通過(guò)建立穩(wěn)定切削模型并選擇穩(wěn)定區(qū)域的切削參數(shù)以避免顫振的產(chǎn)生。由于機(jī)床動(dòng)力學(xué)是非常復(fù)雜的系統(tǒng)，動(dòng)力學(xué)特性會(huì)隨著切削參數(shù)的改變產(chǎn)生瞬態(tài)變化，實(shí)際模型與理論模型往往還存在一定的差距，無(wú)法在工程上廣泛應(yīng)用。因此，利用傳感器信號(hào)對(duì)顫振現(xiàn)象進(jìn)行在線監(jiān)測(cè)與預(yù)報(bào)，具有非常重要的現(xiàn)實(shí)意義。

計(jì)算機(jī)與傳感器技術(shù)的快速發(fā)展，為在線監(jiān)控顫振現(xiàn)象創(chuàng)造了良好的條件。常見(jiàn)的顫振監(jiān)測(cè)傳感信號(hào)包括切削力信號(hào)[6-7]、振動(dòng)加速度信號(hào)[8]、聲壓信號(hào)[9]等。不同的傳感信號(hào)對(duì)不同切削顫振的敏感程度也不同，一般采集的原始信號(hào)也無(wú)法有效地識(shí)別顫振。研究人員通過(guò)在原始信號(hào)中提取能敏感反應(yīng)顫振的特征量，應(yīng)用分類(lèi)技術(shù)，進(jìn)行切削顫振的在線監(jiān)測(cè)研究。西安交通大學(xué)Yao Liu[10]提出了一種基于伺服進(jìn)給電機(jī)電流和振動(dòng)加速度信號(hào)的小波包變換提取顫振特征的在線監(jiān)測(cè)方法，將特征頻帶信號(hào)的標(biāo)準(zhǔn)差和能量比作為顫振監(jiān)測(cè)指標(biāo)。清華大學(xué)張智等人[11]采用加速度傳感器，獲取主軸振動(dòng)信號(hào)，以小波包能量熵值為指標(biāo)，對(duì)銑削加工的穩(wěn)定狀態(tài)進(jìn)行識(shí)別，并取得良好效果。大連理工大學(xué)的周帥[12]獲取聲壓信號(hào)，選取有效值、顫振頻段占能比以及反映能量集中頻段轉(zhuǎn)移特性的特征值作為判斷顫振的特征量，并選取合適的閾值，進(jìn)行多重判斷。Somkia[13]等對(duì)動(dòng)態(tài)切削力信號(hào)求解平均方差并提取相應(yīng)特征，實(shí)現(xiàn)了加工過(guò)程中顫振的在線監(jiān)測(cè)。泰國(guó)朱拉隆功大學(xué) S. Tangjitsitcharoen等[14]基于模式識(shí)別技術(shù)，將切削力的平均變化率作為顫振判別指數(shù)，進(jìn)行車(chē)削顫振的監(jiān)測(cè)和預(yù)報(bào)。

然而，目前研究中并沒(méi)有考慮特征值能否同時(shí)兼顧顫振提前預(yù)警和精確報(bào)警。因此，本文提出兩種新的切削顫振特征提取方法，特征分別著重于顫振早期預(yù)報(bào)的快速性以及狀態(tài)變化時(shí)的強(qiáng)烈敏感性，兩者有效結(jié)合互補(bǔ)，實(shí)現(xiàn)對(duì)切削顫振過(guò)程更快速、精確的監(jiān)測(cè)。

1 研究方法

本文基于均方頻率與經(jīng)驗(yàn)?zāi)B(tài)分解，提出一種切削顫振的特征提取方法。首先，對(duì)銑削試驗(yàn)得到的切削力信號(hào)進(jìn)行快速傅里葉變換，確定顫振發(fā)生的頻域范圍。對(duì)信號(hào)進(jìn)行離散小波包分解與重構(gòu)，提取特征頻段的均方頻率作為特征量T1；其次，對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解得到相應(yīng)的IMF(intrinsic mode function)，通過(guò)計(jì)算IMF與原始信號(hào)的相關(guān)系數(shù)選擇相應(yīng)的分量，最終計(jì)算IMF能量熵作為顫振特征量T2。最后，組成特征向量[T1,T2]作為檢測(cè)指標(biāo)。

1.1 理論基礎(chǔ)

小波分解是一種窗口面積固定但其形狀可改變，即時(shí)間和頻率窗都可改變的時(shí)頻局部化分析方法，由于它在分解過(guò)程中只對(duì)低頻信號(hào)再分解，對(duì)高頻信號(hào)不再實(shí)施分解，使得它的頻率分辨率隨頻率升高而降低。

(1)小波包分解

小波包分解彌補(bǔ)了小波分解只對(duì)低頻信號(hào)再分解，對(duì)高頻信號(hào)不再分解的不足，能根據(jù)信號(hào)特性和分析要求自適應(yīng)地選擇相應(yīng)頻帶與頻譜進(jìn)行匹配。

小波包分解一般采用經(jīng)典的快速迭代算法，其分解算法為

(1)

(2)

重構(gòu)算法為

(3)

(2)經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition)

經(jīng)驗(yàn)?zāi)B(tài)分解是一種適用于處理非線性、非平穩(wěn)信號(hào)的時(shí)頻處理方法。它能夠?qū)?fù)雜的信號(hào)分解成有限個(gè)本征模函數(shù)(IMF)分量，每個(gè)IMF都包含了原始信號(hào)在不同時(shí)間尺度上的局部特征信號(hào)。IMF必須滿(mǎn)足以下兩個(gè)條件① 信號(hào)序列的極大值點(diǎn)與過(guò)零點(diǎn)，兩者的數(shù)目必須一致或只能相差一個(gè)；②在信號(hào)序列中，局部極大值與過(guò)零點(diǎn)所形成的平均包絡(luò)要為零。具體算法流程如下：

①對(duì)原始信號(hào)x(t)的所有極大值、極小值點(diǎn)，利用三次樣條插值函數(shù)分別擬合成上下包絡(luò)線，并求得平均包絡(luò)m1(t)。

②將原始信號(hào)減去平均包絡(luò)得到h1(t)，判斷所得序列是否滿(mǎn)足IMF條件。若滿(mǎn)足，則得到第一個(gè)IMF分量h1(t)，并將x(t)-h1(t)作為新的原始信號(hào)回到步驟1，進(jìn)入循環(huán)；若條件不滿(mǎn)足，則將h1(t)作為原始信號(hào)回到步驟1，再次篩選，直到滿(mǎn)足條件。

③重復(fù)執(zhí)行步驟①與步驟②，當(dāng)執(zhí)行到第n次，若第n個(gè)剩余量為無(wú)法再次分解的單調(diào)函數(shù)，則整個(gè)EMD分解過(guò)程完成，得到全部的IMF分量，如下式所示：

(4)

式中：x(t)為原始信號(hào)、ci為IMF分量、rs為殘余函數(shù)。

(3)均方頻率(mean square frequency)

均方頻率能夠反映信號(hào)在頻域中的變化，當(dāng)顫振發(fā)生時(shí)，頻域的變化是十分顯著的，它是信號(hào)頻率平方的加權(quán)平均，以功率譜的幅值為權(quán)。由于顫振現(xiàn)象的發(fā)生，信號(hào)的主頻帶下移，均方頻率的值也會(huì)急劇下降。均方頻率的計(jì)算公式為

(5)

式中：f表示頻率；S(f)表示功率譜幅值。

(4)能量熵

熵能夠反映一組時(shí)間序列波動(dòng)的規(guī)律性和復(fù)雜性，是序列混亂程度的一種度量。針對(duì)復(fù)雜信號(hào)，可通過(guò)計(jì)算信號(hào)的能量熵[15]來(lái)監(jiān)測(cè)顫振的發(fā)生。

假設(shè)一組信號(hào)經(jīng)過(guò)分解得到M個(gè)分量[S1,S2,…,SM]，計(jì)算每個(gè)分量的能量值：

(6)

式中：Si(t)是原信號(hào)分解的第i個(gè)分量。通過(guò)計(jì)算M個(gè)分量的能量值，得到信號(hào)分量向量[E1,E2,...,EM]。對(duì)能量值進(jìn)行歸一化：

(7)

式中：∑pi=1；∑Ei為總能量。

能量熵的計(jì)算公式如下：

(8)

2 試驗(yàn)研究

銑削顫振試驗(yàn)裝置和方法如下[16]：銑削力采集系統(tǒng)由 Kistler9272測(cè)力計(jì)、5697數(shù)據(jù)采集卡、Kistler-5070A電荷放大器以及PC組成。工件材料為400 mm×100 mm×100 mm的鋁合金，刀具為四齒硬質(zhì)合金球頭銑刀，半徑4 mm，懸伸長(zhǎng)度為40 mm，銑削力采樣頻率為20 000 Hz。通過(guò)改變銑削加工的切削參數(shù)，如：主軸轉(zhuǎn)速n、進(jìn)給速度f(wàn)以及軸向切削深度ap，得到五組加工切削力信號(hào)，其中前三組均為穩(wěn)定狀態(tài)，后兩組出現(xiàn)顫振現(xiàn)象，如表1所示。

表1 實(shí)驗(yàn)條件設(shè)置

以第5組數(shù)據(jù)為例，圖1a為時(shí)長(zhǎng)7 s的加工信號(hào)時(shí)域圖，是從穩(wěn)定加工到顫振現(xiàn)象發(fā)生的一個(gè)過(guò)程。圖1b～d分別為其中穩(wěn)定、過(guò)渡、顫振3個(gè)階段的時(shí)域圖。可以看出，與穩(wěn)定狀態(tài)的切削力信號(hào)相比，顫振狀態(tài)信號(hào)幅值急劇增大，是一種破壞性極大的異常切削狀態(tài)。

對(duì)3個(gè)狀態(tài)的信號(hào)都進(jìn)行快速傅里葉變換，得到相應(yīng)的頻譜圖。由圖可見(jiàn)，圖2a的頻譜整體上幅值較小，且主要集中在低頻段，在65 Hz左右達(dá)到幅值最大。和穩(wěn)定狀態(tài)相比，圖2b在500～1 000 Hz之間幅值有所增加，這說(shuō)明在顫振過(guò)渡階段，切削力信號(hào)開(kāi)始比較明顯地出現(xiàn)能量頻移現(xiàn)象。從圖2c中可以發(fā)現(xiàn)，切削力信號(hào)在500～1 000 Hz頻段之間幅值急劇增大，峰值在788 Hz附近?；趯?duì)信號(hào)時(shí)域和頻域兩方面的分析，可以得到結(jié)論：當(dāng)發(fā)生顫振現(xiàn)象時(shí)，切削力信號(hào)在時(shí)域上的幅值會(huì)急劇增大，頻域上信號(hào)能量會(huì)發(fā)生頻移。

3 結(jié)果分析

3.1 基于均方頻率的顫振特征識(shí)別

通過(guò)上述研究發(fā)現(xiàn)，切削顫振的產(chǎn)生伴隨著能量頻移的現(xiàn)象，在788 Hz處達(dá)到峰值的最大。通過(guò)對(duì)第五組信號(hào)截取包含穩(wěn)定、過(guò)渡以及顫振狀態(tài)的信號(hào)點(diǎn)，通過(guò)小波包變換對(duì)其進(jìn)行4層分解，得到16個(gè)頻率分段,為d1～d16。其中，d2段的頻率范圍為625～1 250 Hz，d2段包含的顫振信息最為豐富，對(duì)d2段的小波包系數(shù)進(jìn)行重構(gòu)，計(jì)算該頻率段信號(hào)的均方頻率。

為了提高顫振特征的普適性，選擇顫振發(fā)生頻段切削力信號(hào)均方頻率MSF與顫振狀態(tài)下均方頻率最大值MSFmax的占比系數(shù)T1作為顫振特征。該特征是個(gè)無(wú)量綱參量，且數(shù)值在區(qū)間[0,1]內(nèi)，對(duì)機(jī)床加工的不同工藝以及不同的切削參數(shù)都具有適用性，能夠有效的對(duì)顫振進(jìn)行識(shí)別。T1的計(jì)算公式如下：

(9)

式中：rmax為均方頻率的最大值。

為了驗(yàn)證顫振特征T1是否能夠敏感反應(yīng)顫振的產(chǎn)生，計(jì)算信號(hào)從穩(wěn)定狀態(tài)過(guò)渡到顫振狀態(tài)的T1，并繪制隨信號(hào)變化的曲線圖，如圖3所示。

從圖3中可以看出，在穩(wěn)定切削狀態(tài)，均方頻率占比系數(shù)整體變化不是很大，主要在0.9～1。在0.4 s左右發(fā)生突變，此時(shí)加工進(jìn)入過(guò)渡階段，T1急劇下降到0.75附近，并且隨著顫振的完全發(fā)生，一直維持在0.75附近，且波動(dòng)不大，證實(shí)了T1的有效性，能夠作為監(jiān)測(cè)并預(yù)報(bào)顫振的指標(biāo)。

3.2 基于EMD能量熵的顫振特征識(shí)別

對(duì)第五組信號(hào)做經(jīng)驗(yàn)?zāi)B(tài)分解，原信號(hào)自適應(yīng)地分解為11個(gè)IMF分量和殘余rs，如圖4所示。相關(guān)系數(shù)可以用來(lái)衡量不同信號(hào)間的關(guān)聯(lián)程度[17]，計(jì)算每個(gè)IMF與原始信號(hào)的相關(guān)系數(shù)。

如表2所示，信號(hào)的能量主要集中在前四階IMF中，而且與原信號(hào)的相關(guān)性也要比其余IMF要大一個(gè)數(shù)量級(jí)，說(shuō)明分解后信號(hào)的信息主要分布在前四階IMF中，IMF5～I(xiàn)MF11可以看作是虛假分量，忽略不計(jì)。

表2 IMF與原信號(hào)相關(guān)系數(shù)

5組加工信號(hào)EMD分解后，計(jì)算前四階IMF的能量熵T2，如圖5a所示，前3組為穩(wěn)定狀態(tài)，后2組為顫振狀態(tài)。發(fā)現(xiàn)穩(wěn)定銑削時(shí)，信號(hào)的能量熵值要遠(yuǎn)遠(yuǎn)大于顫振狀態(tài)下的能量熵值。計(jì)算第5組信號(hào)中穩(wěn)定、過(guò)渡、顫振狀態(tài)下的能量熵，結(jié)果如圖5b所示，可以看出銑削從穩(wěn)定過(guò)渡到顫振，熵的值急劇減小。第5組的分析更加證明了本文所提出的顫振識(shí)別方法的可靠性。

穩(wěn)定銑削時(shí)，得到的切削力信號(hào)經(jīng)過(guò)分解，能量分布比較平均，信息也比較不確定，能量熵高。當(dāng)發(fā)生顫振現(xiàn)象時(shí)，信號(hào)能量會(huì)出現(xiàn)頻移現(xiàn)象，向顫振頻率處集中，此時(shí)能量分布不再平均，信息不確定性變小，能量熵減小。通過(guò)圖5可得能量熵作為特征T2可以有效識(shí)別顫振。

圖6為第5組信號(hào)從穩(wěn)定銑削到顫振過(guò)程中特征T2的變化趨勢(shì)。大約在0.52 s處T2發(fā)生突變，進(jìn)入過(guò)渡狀態(tài)的T2開(kāi)始急劇下降。當(dāng)顫振已經(jīng)完全發(fā)生時(shí)，T2一直維持在0.2附近且波動(dòng)變化不大，證實(shí)了T2可以作為監(jiān)測(cè)與預(yù)報(bào)顫振的指標(biāo)。

經(jīng)過(guò)分析比較，T1、T2雖然都在顫振發(fā)生前有一個(gè)急劇的變化趨勢(shì)，但兩者的變化程度仍有一定差距。顫振發(fā)生時(shí)，T1的大小相較于穩(wěn)定狀態(tài)降低近20%；而T2的大小相較于穩(wěn)定狀態(tài)降低近80%，兩者的變化程度相差近4倍，在顫振預(yù)報(bào)的精確性上T2要更優(yōu)于T1。

3.3 方法總結(jié)

通過(guò)上述的分析，可以得知均方頻率占比系數(shù)能夠更快地實(shí)現(xiàn)顫振預(yù)警，在0.4 s左右特征值就已經(jīng)開(kāi)始發(fā)生突變，但是總體變化率只有20%，過(guò)小的變化率不能作為顫振報(bào)警的決定因素，單獨(dú)使用該種方法具有一定的局限性；EMD能量熵在0.52 s處值才開(kāi)始發(fā)生突變，在早期預(yù)警上略有不足，但是其總體變化率達(dá)到近80%，能夠?qū)崿F(xiàn)精確報(bào)警，與均方頻率形成有效的互補(bǔ)，對(duì)監(jiān)測(cè)切削顫振過(guò)程，反應(yīng)更加快速、精確。

4 結(jié)語(yǔ)

基于小波包分析與經(jīng)驗(yàn)?zāi)B(tài)分解，提出了兩種切削顫振的特征提取方法。該兩種方法分別將小波包顫振頻帶的均方頻率占比系數(shù)以及IMF分量能量熵作為特征T1、T2。研究表明T1、T2能夠有效識(shí)別銑削加工過(guò)程中穩(wěn)定、過(guò)渡以及顫振3種狀態(tài)，形成有效互補(bǔ)且通用性好的新方法。后續(xù)工作中可以基于特征向量[T1,T2]，通過(guò)人工智能分類(lèi)算法提高加工狀態(tài)的分類(lèi)識(shí)別精度，實(shí)現(xiàn)臥式加工中心銑削加工穩(wěn)定性的在線監(jiān)測(cè)和預(yù)報(bào)警。