基于瞬時頻率估計與Vold-Kalman濾波的銑削顫振識別

汪曉姍， 彭志科， 陳是扦

(上海交通大學 機械系統與振動國家重點實驗室，上海 200240)

顫振是金屬加工中出現的自激振動，嚴重制約了加工精度和加工效率。隨著機床加工柔性化的日益發展，要求加工能對不同的工件在不同的條件下進行，因此不能從根本上杜絕顫振。近年來，于英華等[1-2]很多國內外學者開始著力于研究顫振的在線監測和預報控制。

銑削顫振的形成是一個孕育的過程，在顫振早期，監測系統搶先把銑削過程即將發生顫振的預兆先通報給控制系統，是實現顫振預報的關鍵。各種各樣的傳感信號可以被用來采集顫振信號，比如動態銑削力信號[3-4]，加速度信號[5]，力矩信號[6]等。Kuljanic等[7]提出了利用多傳感器來監測識別顫振，這種方法提高了顫振監測的精確性和魯棒性；Tansel等[8]通過銑削位移信號的諧波分量來監測顫振；Tarng等[9]在發生顫振時切削力信號在的頻譜上的窄帶特性來預測顫振。

工程中，直接根據傳感信號難以識別顫振，因此需要從采集來的信號中獲得對故障敏感的特征量。時域中主要采用的是概率統計分析，如互相關系數[10]、標準差[11]等；頻域中，遲玉倫等[12]基于通過識別銑削力的譜峰頻率來監測顫振。然而，這些傳統的方法只適合分析平穩信號，而且無法準確預報顫振的發生。對于成分復雜的顫振信號，需要將其分解成單分量再進行分析。Cao等[13]應用EMD分解銑削力信號，提取顫振的特征，但是無法克服EMD分解的模態混疊。該方法易產生虛假的信號分量，干擾故障特征的診斷；Zhang等[14]提出用VMD分解銑削力信號，但是VMD本質上通過自適應的頻域濾波器組來提取信號，其分解銑削力信號得到的子分量沒有物理意義。Vold-Kalman濾波方法[15]基于數據方程誤差和最小化結構，有效地進行多分量信號分解。不僅可以避免EMD分解帶來的模式混淆，而且可以有效的從復雜的顫振信號中提取單分量。Vold-Kalman濾波提取單信號首先需要根據時頻脊線估計單分量信號的瞬時頻率。經典的瞬時頻率估計方法中，過零法和相位差分法較為簡單，但是誤差大。HHT繼承了EMD的缺點。在參數化時頻分析的基礎上，Yang等[16]提出了基于優化頻譜集中性指標，對信號瞬時頻率參數估計的方法，該方法具有抗噪性強、精確度高等特點。

顫振產生時，信號的頻率成分和能量分布發生改變。穩定銑削時，顫振的能量分布在噪聲、主軸的轉速頻率及其諧波處。發生顫振時，振動的能量逐漸轉移到顫振頻率處。因此本文提出基于頻域能量分布改變的顫振識別方法。針對復雜的銑削顫振信號，考慮Vold-Kalman濾波方法在單分量提取和優化頻譜集中性指標在瞬時頻率計算方面的獨特優勢，本文提出了基于兩者的多分量信號分解方法。采用分解后的子信號在頻域上的能量分布作為顫振識別的特征。多次銑削實驗驗證了本文所提出的方法的有效可行。

1 多分量信號分解

1.1 基于優化頻譜集中性指標估計信號瞬時頻率

該過程可以描述為，選取一個目標分量，采用合適的變換核參數構造對應的旋轉算子，將原信號旋轉變換到它的一個旋轉域。為了在上述旋轉域評估分量能量的集中性，提出頻譜集中性指標。旋轉域中能量越集中，集中性指標越大，可認為變換核參數與分量模型的相位參數越相近。

單分量解析信號可用多項式相位信號形式表示[17]

(1)

(2)

sd(t;C)=s(t)Φ(t;C)

(3)

sd(t;C)=a(t)exp(j(2πc0t+φ0))

(4)

此時能量集中于c0處， 為衡量頻譜集中性，定義頻譜集中性指標

(5)

式中：E(·)表示期望算子； F (·)表示傅里葉變換，以SCI最大為參數估計的指標。

(6)

本文，估計參數采用粒子群優化算法，通過找到解調信號的最大譜峰值對應的頻率來估計初始頻率參數

(7)

通過式(6)、式(7)估計瞬時頻率的所有參數(c0,…,ck)。

1.2 Vold-Kalman濾波

Vold-Kalman濾波方法可以從復雜多分量信號中提取有效的頻率成分。與傳統的濾波方法相比，該方法避免了由時域至頻域變換帶來的相位偏差，可以將復雜多分量信號分解為單分量。Vold-Kalman濾波的中心頻率能根據瞬時頻率進行自適應調節，因此能夠有效分離在時頻域內鄰近甚至交叉的信號分量[18]。

調制信號可以表示為

(8)

式中：Ak(t)為第k個分量的幅值包絡；Θk(t)為載波信號；k為階次。

(9)

振幅包絡線是一個低頻率調制的載波信號，可以表示為一個低階多項式。對于離散信號

(10)

Ak(n-1)-2Ak(n)+Ak(n+1)=εk(n)

(11)

對于實際信號，Ak(0)=0。 當n=1時，Ak(0)-2Ak(1)+Ak(2)=εk(1)； 當n=N時， 其中是信號的長度。 式(11)的矩陣的形式

(12)

矩陣形式為

MA=ε

(13)

實際測試信號

(14)

為噪聲或誤差項，式(14)可以寫成矩陣形式

X-BA=δ

(15)

優化頻譜集中性指標， 估計信號的瞬時頻率ωk(n)， 從而得到B， 根據(13)和(16)得到各信號分量的幅值包絡矩陣A， 由矩陣C和A可以重構各信號分量

S=AC

(16)

信號的非線性調頻分量分解的步驟如下：

(1) 優化頻譜集中性指標，估計分量的瞬時頻率；

(2) 根據估計到的瞬時頻率應用Vold-Kalman濾波器濾出分量ui；

(3) 從信號中減去復原的信號分量作為初始信號，重復以上步驟，得到信號的各個分量。

2 特征提取

為了顫振預報及時并且精確，正確提取特征是關鍵。所選擇的特征參數要充分反映銑削顫振過程的本質，同時還必須考慮信號采集和數據處理的簡易可行，使監測系統可以在幾秒時間內完成計算工作，向控制系統發出預報信號。

頻帶熵、信息熵、近似熵等的引用使信號復雜性表征成為可能，但是它們都存在缺點：在單一尺度上對時間序列進行分析、需要比較自身數據段且參數改變時結果一致性差。因此，本文采用多分量信號分解和熵結合的方法。在穩定銑削時，除了主軸轉速頻率及其諧波處，其他頻率處具有能量小、分布相對平均、不確定的特點。當發生顫振時，除了主軸轉速頻率及其諧波處，在相應的頻帶內就會出現相應共振頻率，能量便會集中在此頻率帶內，使能量分布的不確定性減少。因此可以通過計算分解得到的各子信號的能量分布以檢測顫振。信號經過分解可以得到n個單分量u1,u2,…,un， 計算出分別對應的能量

(17)

假設殘余分量可以忽略，由于本文的分解方法具有正交性，n個子信號的能量之應該等于原始信號的能量。由于各個子信號包含不同的頻率成分和能量，從而形成在銑削力信號在頻域的能量分布。由此引入能量熵銑削過程中的振動信號包括：由于銑刀旋轉引起的周期信號部分、由于再生效應引起的顫振和由于系統噪聲等引起的隨機擾動部分。從實驗信號中提取顫振分量對分析顫振有重要意義。本文提出的顫振監測方法如圖1(b)所示。首先，動態銑削力信號經過Vold-Kalman濾波器去除由于主軸旋轉產生的信號分量。然后優化頻譜集中性指標對信號進行瞬時頻率估計。Vold-Kalman濾波將復雜的多分量信號分解為單分量，從時域中提取對分析顫振有意義的信號分量。最后計算分量信號的能量熵作為識別顫振的特征。

(18)

3 實驗設置

實驗裝置圖1(a)所示，所有的銑削實驗都是在數控銑削機床上完成的，動態銑削力信號測量系統由Kistler 9272測力計、Kistler 5070A電荷放大器、5697數據采集卡和PC組成，采樣頻率為20 000 Hz。工件材料為400 mm←100 mm←100 mm的鋁合金， 刀具為四齒硬質合金球頭銑刀， 半徑4 mm， 懸伸長度為40 mm。 顫振的產生與銑削參數有直接的關系， 因此可以通過改變主軸轉速n、進給速度f、軸向銑削深度ap獲取感興趣的實驗數據。

本文進行了5次實驗，實驗條件如表1所示。 Test1、Test2、Test3、Test4構成第一組實驗； Test5分別在穩定、過渡、顫振時采樣，這3段采樣信號構成第二組實驗。

圖1 實驗裝置及方法流程圖Fig.1 The experimental platform and the flow chart

狀態n/(r·min-1)f/(mm·min-1)ap/mm第一組Test1穩定4 2005806.6Test2穩定3 6002502.6Test3顫振4 0005806.6Test4顫振4 0006406.7第二組Test5過渡4 0005806.3

4 結果分析

對于第一組實驗，分別對Test1、Test2、Test3、Test4進行多分量信號分解，計算子信號的能量在頻域上的分布。對比分析在不同實驗條件下，穩定和顫振的子信號的能量熵。

圖2(a)為Test5采集到的進給方向動態銑削力的時域波形，Test5經歷了從穩定銑削過渡到顫振的過程，實驗前0.6 s銑削穩定，0.6 s～0.8 s信號振幅急劇增大， 0.8 s開始顫振。 分別自0.2 s、0.8 s、2 s起， 往后截取5 000個采樣點長的信號記作Sample1(穩定)、Sample2(過渡)、Sample3(顫振)，作為第二組實驗。對比Sample1、Sample2、Sample3經過信號分解得到的子信號的能量分布，分析在實驗條件不變的情況下，隨著銑削從穩定過渡到顫振，子信號在頻域的能量熵的變化。

圖2 Test5及各采樣點時域波形Fig.2 Time domain waveform of test5 and sampling points

4.1 數據分析

第一組原始實驗數據的時域波形和頻譜如圖 3所示，Test1(圖3(a))、Test2(圖3(b))銑削穩定，頻譜由主軸轉速頻率及其諧波、噪聲組成；改變銑削條件，動態銑削力振幅急劇增加，發生顫振(圖3(c))，頻譜上的出現能量遠遠高于其他頻率的顫振峰值；繼續增加銑削深度，顫振的更加劇烈，如圖3(d)。

圖3 Test1、Test2、Test3、Test4的波形和FFTFig.3 Vibration signal in time domain andfrequency domain of Tests

第二組原始實驗數據的頻譜如圖4(a)、圖4(b)、圖4(c)所示，f表示主軸轉速頻率，fc表示顫振頻率。銑削穩定時(圖4(a))，頻譜由主軸轉速及其諧波成分和分布不確定的噪聲組成；過渡狀態(圖4(b))頻譜除了主軸轉速諧波成分和噪聲，還出現了顫振分量和它的各階諧波成份，此時工件表面還沒有明顯的顫振；隨著顫振程度增加(圖4(c))，顫振頻率及其諧波對應的峰值也急劇增長，信號的能量除了集中在主軸轉速頻率處，都匯集在fc及其諧波上，且顫振峰值還再繼續增加。

在分解信號、分析子信號能量頻域分布之前，為了避免主軸轉速及其諧波分量對顫振識別產生的干擾，先估計這些分量的瞬時頻率，用Vold-Kalman濾波提取此分量信號并從原始信號中剔除這一部分。對預處理過的銑削信號進行FFT，穩定銑削時頻譜(圖4(d))只剩下在頻率域上分布均勻且峰值較低的隨機振動；過渡到顫振時(圖4(e)、圖4(f))， 頻譜在fc及2倍的fc處出現能量很高的峰值，除此之外還有能量比較低的噪聲。

通常動態銑削力信號具有非線性、非平穩等特征，直接進行傅里葉變換無法揭示頻率分量隨時間的變化規律；短時傅里葉變換、小波變換等能夠建立信號在時間和頻率上的分布，然而它們是先驗性的自適應能力差。并且它們只能用來離線判斷銑削狀態，無法運用在實際中。預報銑削顫振要準確、及時，一旦出現顫振發生的預兆，計算機可以識別并且迅速的啟動控制系統，調整轉速、銑削用量等，將顫振控制在孕育階段。故障振動信號的能量分布會隨著顫振發生相應的改變，故采用分解后子信號在頻域分布的熵值的范圍可用來判斷是否顫振。

圖4 Sample1、Sample2、Sample3濾波前后的FFTFig.4 FFT of Sample1、Sample2、Sample3 beforeand after filtering

首先，Vold-Kalman濾波器波除去信號的轉速頻率及其諧波成分。優化頻譜集中性指標來估計發生顫振時預處理過的信號的瞬時頻率參數。用估計到的參數將非線性調頻信號解調為平穩信號。最后Vold-Kalman濾波器提取解調單分量信號。本文對實驗信號精細分解，按照分解后單分量子信號能量從高到低，得到這些分量集中度較好的頻譜表示。由于主要的故障信息集中在前幾個子信號中，剩余信號能量遠小于原始信號總能量，可以忽略，因此本文選用了前8個子信號，8個子信號的能量之和應該恒等于原始信號的總能量。由于篇幅限制，本文僅展示預處理后的第二組信號的Sample1、Sample3的分解結果。Sample1的經過分解后的各子信號按能量從高到低如圖5(a)所示，均為隨機分布的噪聲，能量分布相對平均和不確定；Sample3經過分解后的各子信號能量從高到低如圖5(b)所示，其中u1，u2能量遠高于其他各子信號能量，信號頻譜的幅值范圍均在0～4。u1的峰值正好在銑削系統的固有頻率附近，證明本實驗出現的是典型的再生顫振。u2為顫振頻率的2倍頻。剩余各階子信號是被齒頻調制的顫振分量或噪聲，其時域波形的振幅小，能量明顯低于u1和u2。

經過預處理的動態銑削力信號經過分解后得到多個分量，且各個分量的瞬時頻率在頻譜上不相交，子信號的能量隨分解階數的增加而降低。在信噪比低的情況下，基于頻譜集中性指標估計分量的瞬時頻率較精確，抗噪性強。即使對于能量較小的單分量的瞬時頻率參數的估計也比較準確。

圖5 本文方法分解結果Fig.5 Decomposition results of this method

4.2 顫振識別

隨著顫振的發生，振動的能量發生轉移，所以振動的概率分布也發生了變化，這些都可以用熵來進行判斷。由于各階子信號含不同的頻率成分且具有不同的能量，分別求它們的能量E={E1,E2,…,En}。首先根據式(18)將各子信號的能量歸一化。分析第一、二組實驗經過分解后的能量分布，分別如圖6(a)，圖6(b)所示。對于第一組實驗，將分解得到的子信號歸一化，u1～u8的能量從高到低，分布相對均勻(圖6(a))，熵值較大；而Test3能量集中在u1和u2，所以熵值相對較小；Test4顫振更加劇烈，隨著顫振的程度增加，u3～u8所占能量幾乎為零，所以計算得到的熵值最小。計算Test1、Test2、Test3、Test4的能量熵結果如圖7(a)所示。對于第二組實驗，短暫的穩定銑削過程中，子信號能量分布較為均勻。隨著顫振的出現，能量急劇轉移到u1和u2，顫振劇烈時，u1和u2匯集了幾乎全部的能量。根據式(19)，計算Sample1、Sample2、Sample3的能量熵結果如圖7(b)所示，可以看出銑削從穩定過渡到顫振，熵的值急劇減小。第二組的分析更加證明了本文所提出的顫振識別方法的可靠性。

圖6 第一組、第二組歸一化后各子信號能量分布Fig.6 Normalized energy ration of the sub-signals

穩定銑削時，預處理后的銑削力信號經分解，得到的子信號能量分布較平均和不確定，能量熵較高。當出現顫振時，在相應的頻帶內就會出現相應的顫振頻率。此時，能量集中在顫振頻率上，使能量分布的不確定性減小，從而熵減小。通過第一組、第二組實驗的分析，可見本文提出的方法可以有效識別顫振。然而由于實驗條件的限制，本文沒有給出發生顫振時能量熵的閾值。

圖8所示為Test5從穩定銑削到顫振的過程，可見在0.6 s顫振發生時，能量熵明顯變小，證實了能量熵可以作為監測顫振的指標。本文的不足之處在于，由于實驗前后時間跨度大，無法提供工件表面加工質量前后對比照片。

圖7 第一組、第二組能量熵Fig.7 Energy entropy of the tests

圖8 Test5能量熵隨銑削力的變化Fig.8 The vibration signal and the relatedenergy entropy in Test 5

5 結 論

本文提出了一種簡單有效的銑削顫振識別方法。動態銑削力信號具有不平穩、頻率時變、噪聲大等特點，分解信號前先經過Vold-Kalman濾波器去除信號中的主軸轉速及其諧波分量，再通過優化頻譜集中性指標來估計信號瞬時頻率參數，用估計到的參數將非線性調頻信號解調為平穩信號，最后用Vold-Kalman濾波器提取單分量信號，并得到能量從高到低的各階子信號。該信號分解方法抗噪性好，對實驗信號分解精細，并得到這些分量集中度較好的頻譜表示。隨著銑削條件變化，經過分解得到的各子信號，其能量分布也呈現變化。因此，本文提出表征能量分布特征的能量熵作為顫振識別的特征。在穩定銑削時，實驗信號進行預處理后，其分解得到的各階子信號主要為噪聲，因此能量熵小。隨著顫振越來越劇烈，振動的能量集中到含有顫振信息的子信號上，因此能量分布從分散變得集中起來，能量熵變大。經過實驗驗證，該方法簡單有效，尤其在顫振發生的早期識別精度高，是實現顫振在線實時控制的關鍵。