李鵬偉，戴士杰，張文華，母嘉恒，胡天榮

（河北工業(yè)大學(xué)，天津 300400）

隨著機器人智能控制與各個學(xué)科的交叉融合，在航空發(fā)動機葉片表面磨拋中高效率、高精度的機器人加工備受青睞。航空發(fā)動機葉片屬于典型的薄壁結(jié)構(gòu)零件，在機器人高速磨削加工中很容易發(fā)生顫振。機器人加工系統(tǒng)的顫振是發(fā)生在刀具和工件之間的相互振動，是一種復(fù)雜的彈性動力學(xué)現(xiàn)象，主要由機器人加工系統(tǒng)整體結(jié)構(gòu)從周期性切削力中吸收能量自激產(chǎn)生[1]。顫振不僅會造成磨削力波動，而且會嚴(yán)重降低加工表面質(zhì)量[2]。為了防止顫振對磨削系統(tǒng)破壞，早期、快速、準(zhǔn)確的顫振檢測至關(guān)重要。

顫振檢測的過程主要分為3部分：數(shù)據(jù)采集、特征提取、狀態(tài)識別[3]。其中，在顫振信號的數(shù)據(jù)采集過程中，可以利用各種加工信號，例如，聲音信號采集簡單可靠，但容易受到外界噪聲的影響；力信號對切削狀態(tài)的變化反應(yīng)迅速，但傳感器安裝位置的選擇比較困難；加速度信號對振動的表現(xiàn)最為直觀，但采集加速度信號需要昂貴的儀器[4-7]。在顫振信號的特征提取過程中，可以利用信號處理的方法構(gòu)造特征指標(biāo)，該指標(biāo)為后期定義顫振判據(jù)提供理論依據(jù)。Chen等[8]依靠數(shù)據(jù)處理算法，從提取的多個特征參數(shù)中自動選擇有效的顫振識別特征參數(shù)。Yao等[9]利用顫振信號的小波變換標(biāo)準(zhǔn)差和小波包變換能量比作為二維特征向量進行顫振檢測。在顫振的狀態(tài)識別過程中，實時性是顫振識別的一個重要要求，這意味著顫振應(yīng)該事先被識別出來。Chen等[10]提出一種基于EEMD和特征選擇的多通道顫振智能檢測策略，但煩冗的計算量會導(dǎo)致系統(tǒng)實時性較差。Li等[11]通過經(jīng)驗?zāi)B(tài)分解成功提取顫振識別的特征向量，該方法用時不足1 s。

目前，關(guān)于加工顫振檢測的研究主要圍繞銑削、車削等展開，而較少關(guān)注機器人磨削[12]，而且大部分檢測方法仍然不能在工件損壞之前識別顫振，達(dá)不到防止工件損壞的核心目標(biāo)。因此，本文提出一種基于排列熵的機器人磨削顫振檢測方法：首先，采集機器人末端的顫振加速度信號；然后，利用滑動窗口算法從原始信號中提取連續(xù)的排列熵值；最后，通過排列熵經(jīng)驗閾值判斷末端執(zhí)行器是否發(fā)生顫振。在此過程中，采用基于序數(shù)模式的排列熵算法，利用序數(shù)模式及其分布規(guī)律，簡化了排列熵計算過程，提高了排列熵特征的計算效率。仿真和試驗結(jié)果表明，排列熵值的驟降時刻遠(yuǎn)遠(yuǎn)早于顫振完全產(chǎn)生時刻，說明該方法能夠在顫振孕育過程的初始階段檢測出顫振，為后續(xù)采取抑制措施提供更多的反應(yīng)時間。

1 基于序數(shù)模式的排列熵

1.1 序數(shù)模式

序數(shù)模式描述了時間序列相鄰值之間的順序關(guān)系?？稍O(shè)時間序列 （xm，xm-t，…，xm-dt） 的d階序數(shù)模式為序數(shù)模式的計算表達(dá)式為

式中，m為時間序列中數(shù)據(jù)的數(shù)量；h=1，2，…，d；t為時間延遲，代表相鄰序數(shù)模式的中點間的距離；d為階數(shù)，代表一個序數(shù)模式中含有（d+ 1）個數(shù)據(jù)。

通過給每一種序數(shù)模式分配代號，可以更加方便地觀察序數(shù)模式的分布情況，其中代號按式 （2）進行一對一的分配。

例如，當(dāng)階數(shù)d= 2時，會產(chǎn)生6種序數(shù)模式類型，分別對應(yīng)數(shù)字代號“0、1、2、3、4、5”，如圖1所示。

圖1 d = 2時序數(shù)模式的類型Fig.1 Types of ordinal patterns at d = 2

1.2 基于序數(shù)模式的排列熵計算

排列熵是一種衡量時間序列復(fù)雜程度的特征，能夠準(zhǔn)確地定位系統(tǒng)發(fā)生突變的時刻，可以對突變信號的微小變化時刻進行放大。

滑動窗口內(nèi)的時間序列 （xm，xm-1，…，xm-M-dt+1） 在時間延遲t和階數(shù)d下的排列熵可以根據(jù)式 （3）進行計算。

式中，M為窗口中存在的序數(shù)模式的數(shù)量；qj= #[k∈{m，m-1，…，m-M+ 1}|（xk，xk-t，xk-dt）存在序數(shù)模式j(luò)]。

設(shè)由10個值組成的時間序列為 （x1，x2，…，x10），此時間序列中時間延遲t= 2、階數(shù)d= 2的序數(shù)模式如圖2所示。

圖2 滑動窗口中的序數(shù)模式Fig.2 Ordinal patterns in the sliding window

由圖2可知，兩個窗口存在4種序數(shù)模式的“重疊”部分。根據(jù)式 （3）可在窗口1和窗口2中相應(yīng)地求得排列熵（9）和（10），即

1.3 高效計算連續(xù)序數(shù)模式

圖3 序數(shù)模式“0”以及其連續(xù)序數(shù)模式的可能類型Fig.3 Ordinal pattern “0” and its possible types of consecutive ordinal patterns

由圖3可知，對于序數(shù)模式“0”的下一個值有3個可能的位置，則下一個連續(xù)序數(shù)模式有3種可能的類型：“0、3、4”。

同理可得，可以對其他階數(shù)的給定序數(shù)模式的下一個連續(xù)序數(shù)模式的可能類型進行分析。當(dāng)前序數(shù)模式已知，通過下一個數(shù)據(jù)與當(dāng)前序數(shù)模式中數(shù)據(jù)的大小比較，可求得下一個序數(shù)模式的類型，如式 （4）所示。

式中，ωd表示根據(jù)當(dāng)前序數(shù)模式來確定下一個序數(shù)模式的函數(shù)；I表示 （xm，xm-t，…，xm-(d–1)t）中有多少個值大于或等于xm+t。

根據(jù)式 （4）和 （5），可以求得所有階數(shù)序數(shù)模式的連續(xù)序數(shù)模式。d=2階的序數(shù)模式及其連續(xù)模式如圖4所示。例如，對于序數(shù)模式“0、1、2”，下一個值有3個可能的位置“I= 0、1、2”，也就有3個可能的連續(xù)序數(shù)模式“0，3，4”。對于序數(shù)模式“3、4、5”，下一個值有3個可能的位置“I= 0，1，2”，也就有3個可能的連續(xù)序數(shù)模式“1，2，5”。同理可得，每個代號為 （d+1）（k-1），…，（d+1）k-1的序數(shù)模式具有相同的連續(xù)序數(shù)模式，其中k=1，2，…，d!。根據(jù)上述方法，可以求得所有的序數(shù)模式以及其連續(xù)序數(shù)模式。

圖4 階數(shù)d =2的所有序數(shù)模式以及其連續(xù)序數(shù)模式Fig.4 All ordinal patterns of order d =2 and their consecutive ordinal patterns

排列熵是在具有最大重疊的大小為M的滑動窗口中計算的，即連續(xù)窗口的第1個值是前一個窗口的第2個值。連續(xù)的窗口 （xm-1，xm-2，…，xm-M-dt）和 （xm，xm-1，…，xm-M-dt+1）在值xm和xm-M-dt上有所不同，因此窗口中的序數(shù)分布在序數(shù)模式（m）和（m-M）的出現(xiàn)頻率上有所不同。所以在已知排列熵（m-1）的情況下，排列熵（m）的計算方法可表示為

1.4 基于序數(shù)模式的排列熵的計算效率

本文提出的基于序數(shù)模式的排列熵計算方法通過比較信號數(shù)據(jù)中相鄰數(shù)值的大小確定序數(shù)模式的類型（序數(shù)模式類型的定義已經(jīng)給出）。利用圖4的序數(shù)模式分布規(guī)律可快速獲得連續(xù)的序數(shù)模式分布情況。最后通過連續(xù)序數(shù)模式的類型，以及式 （1）～（6）計算基于序數(shù)模式的排列熵值，并利用滑動窗口算法得到連續(xù)的基于序數(shù)模式的排列熵值。

分別利用傳統(tǒng)排列熵算法和所提出的基于序數(shù)模式排列熵算法對不同信號長度以及不同階數(shù)的隨機信號進行分析，然后對兩種算法排列熵提取的計算時間進行比較 （計算時間是多次運行的平均值），結(jié)果如表1所示。

表1 不同方法計算排列熵值的執(zhí)行時間Table 1 Execution time for calculatating permutation entropy using different methods

由表1可知，當(dāng)原始信號長度加倍時，兩種計算方法的執(zhí)行時間都有所增加。然而，相比于傳統(tǒng)排列熵計算方法，基于序數(shù)模式排列熵計算方法的計算時間要低于一個數(shù)量級。因此，基于序數(shù)模式的排列熵計算方法效率更高，更有利于實現(xiàn)顫振的及時檢測。

基于序數(shù)模式的排列熵方法對上述有效計算進行總結(jié)，計算方案如圖5所示。其中，W為時間序列的長度。

圖5 基于序數(shù)模式的排列熵計算方案Fig.5 Permutation entropy calculation scheme based on ordinal pattern

2 仿真算例分析

顫振的發(fā)生一般有3個階段，即平穩(wěn)階段、過渡階段和顫振階段[13]。顫振信號通常表現(xiàn)為單次諧波特征，因此，當(dāng)顫振發(fā)生時，加速度信號可用單頻諧波函數(shù)描述。為了驗證所提出的磨削顫振檢測方法的有效性，現(xiàn)利用仿真軟件，直觀地表達(dá)出顫振從穩(wěn)定到不穩(wěn)定逐步發(fā)展的過程，且將存在平穩(wěn)、過渡和顫振3個階段的顫振信號模擬出來，顫振加速度模擬信號如圖6所示。

圖6 經(jīng)歷平穩(wěn)、過渡和顫振階段的加速度模擬信號Fig.6 Simulated acceleration signal through smooth, transition and chatter phases

本研究選取階數(shù)d= 3，時間延遲t= 1，滑動窗口大小M= 512對顫振原始信號進行窗口滑動計算，每次對每個窗口范圍的信號進行一次排列熵值計算，最終得到連續(xù)的排列熵值數(shù)據(jù)，可以更直觀地觀察加速度信號變化與排列熵值變化的差異對比情況?；瑒哟翱谒惴ㄓ嬎愠龌谛驍?shù)模式排列熵的動態(tài)變化情況如圖7所示。

圖7 顫振信號的排列熵特征Fig.7 Permutation entropy characteristics of chattering signal

由圖7可知，在平穩(wěn)階段的排列熵特征曲線大致保持在水平穩(wěn)定狀態(tài)。由于周期性振動的積累，顫振逐漸產(chǎn)生，排列熵特征曲線發(fā)生驟降，而且驟降發(fā)生在顫振過渡階段的前期。排列熵特征驟降可以作為顫振即將發(fā)生的征兆，而且排列熵驟降時刻要遠(yuǎn)遠(yuǎn)早于顫振完全產(chǎn)生時刻，因此可以利用排列熵驟降作為顫振檢測的特征。

3 試驗驗證

3.1 試驗平臺搭建

本節(jié)通過航空發(fā)動機葉片機器人磨削試驗來驗證本文所提出的磨削顫振檢測的可靠性，試驗平臺的搭建如圖8所示，加工工件是材質(zhì)為鎳基合金的航空發(fā)動機壓氣機葉片，磨具為粒度500目的柔性砂輪。

圖8 航空發(fā)動機葉片機器人磨拋試驗平臺Fig.8 Robot grinding experimental platform for aero-engine blade

機器人磨削顫振信號采集裝置如圖9所示，相比于末端執(zhí)行器，三軸加速度傳感器的體積及重量極小，在此可忽略加速度傳感器重量對末端顫振的影響。其他試驗設(shè)備及相關(guān)參數(shù)如表2所示。

表2 試驗設(shè)備及其主要參數(shù)Table 2 Experimental equipment and its main parameters

圖9 機器人磨削顫振信號采集裝置Fig.9 Acquisition device of chattering signal of robot grinding

3.2 顫振信號分析

通過改變磨削加工參數(shù)，采集機器人末端在不同加工狀態(tài)下的顫振信號，磨削加工參數(shù)如表3所示。

表3 試驗中的磨削加工參數(shù)Table 3 Grinding parameters in the experiment

試驗由安裝在磨削末端執(zhí)行器上的三軸加速度傳感器測量顫振加速度信號，通過HY-DAS數(shù)據(jù)系統(tǒng)采集顫振加速度信號，采用上位機對顫振加速度信號進行實時排列熵提取。其中，利用本文所提出基于序數(shù)模式的排列熵算法進行磨削顫振檢測，即選取階數(shù)d= 3、時間延遲t= 1、滑動窗口大小M= 512提取排列熵特征。根據(jù)表3中的加工參數(shù)進行航空發(fā)動機葉片機器人磨削，在試驗結(jié)果中選取幾組代表性數(shù)據(jù)進行分析，得到顫振原始信號及其連續(xù)排列熵特征，如圖10所示。

圖10 顫振原始信號和排列熵Fig.10 Chattering raw signal and permutation entropy

圖10（a）中排列熵值在0.23 s時低于0.95，早于顫振完全發(fā)生階段0.41s；圖10（b）中排列熵值在0.2 s時低于0.95，早于顫振完全發(fā)生階段0.5 s；圖10（c）中排列熵值在0.88 s時低于0.95，早于顫振完全發(fā)生階段0.52 s。由此可知，排列熵值低于0.95時刻早于顫振完全發(fā)生時刻約0.48 s。

由圖10可知，當(dāng)顫振未發(fā)生時，排列熵特征處于水平穩(wěn)定的狀態(tài)。隨著顫振的產(chǎn)生，原始信號幅值隨之增加，在過渡階段的初始位置，排列熵值發(fā)生驟降 （虛線為排列熵驟降時刻），此時即可判斷出磨削過程中發(fā)生了顫振，而且排列熵值驟降時刻遠(yuǎn)遠(yuǎn)早于顫振完全產(chǎn)生時刻，即可以在較早階段檢測出顫振。因此以排列熵值的驟降作為顫振檢測的標(biāo)準(zhǔn)，能夠在顫振完全產(chǎn)生前0.48 s檢測出顫振，為磨削系統(tǒng)采取顫振抑制措施提供更多的反應(yīng)時間。該試驗證明所提出的方法在機器人磨削顫振檢測方面是有效的。

3.3 特征對比

基于振動信號時域分析理論，選取方差、裕度因子作為顫振檢測的特征量來反映磨削顫振信號的變化趨勢，利用滑動窗口算法提取方差、裕度因子特征曲線，并與排列熵特征進行對比，結(jié)果如圖11所示。

圖11 顫振檢測特征：排列熵、方差、裕度因子Fig.11 Features of chatter detection: permutation entropy, variance,margin factor

由圖11可知，方差特征曲線會隨著信號幅值的增減而增減，其變化受原始信號幅值的影響較大，而且方差特征曲線的變化幅度小于原始信號變化幅度，所以其檢測準(zhǔn)確性較差。裕度因子特征曲線變化不穩(wěn)定，曲線波動范圍較大，所以其檢測精度較差。相比于方差和裕度因子，小波包熵特征驟降時刻要更早 （圖11中的虛線為各個特征的轉(zhuǎn)折時間點），而且其波動較小、變化穩(wěn)定。

綜上所述，排列熵特征對顫振發(fā)生的反應(yīng)速度更快，因此利用排列熵特征能夠更早地檢測出顫振，為磨削系統(tǒng)采取顫振抑制措施提供更多反應(yīng)時間。

3.4 顫振狀態(tài)識別

顫振狀態(tài)識別最簡單的方法就是特征閾值法，通過對比某個特征量在平穩(wěn)階段與顫振征兆出現(xiàn)階段的變化情況，設(shè)定一個經(jīng)驗型的閾值，一旦超過或低于預(yù)先設(shè)定的閾值則判定為出現(xiàn)顫振，即可提醒系統(tǒng)采取顫振抑制措施。根據(jù)圖10中的排列熵特征的動態(tài)變化規(guī)律，可設(shè)置經(jīng)驗閾值為0.95，即當(dāng)排列熵低于0.95則判定此時加工系統(tǒng)中發(fā)生顫振。

綜上所述，機器人磨削末端顫振檢測方法的流程如圖12所示。首先設(shè)置磨削加工參數(shù)，利用機器人帶動末端執(zhí)行器進行航空發(fā)動機葉片磨削，并通過安裝在末端執(zhí)行器上的加速度傳感器對磨削顫振信號進行采集；然后利用序數(shù)模式進行排列熵特征連續(xù)提??；最后當(dāng)排列熵特征值小于0.95時，即判定此時機器人末端正在產(chǎn)生顫振，可為系統(tǒng)發(fā)出警告，提醒系統(tǒng)采取顫振抑制措施。

圖12 顫振檢測方法流程圖Fig.12 Flow chart of chatter detection method

試驗結(jié)果表明，當(dāng)顫振信號加速度幅值剛開始增加時排列熵值就會發(fā)生驟降，因此，在過渡階段的初始位置將顫振檢測出來 （排列熵值低于0.95時即可認(rèn)定此時的磨削過程中發(fā)生了顫振）為進行抑制操作提供了更多反應(yīng)時間。

4 結(jié)論

針對航空發(fā)動機葉片機器人磨削顫振檢測不及時導(dǎo)致顫振抑制滯后的問題，本文提出了一種基于排列熵的航空發(fā)動機葉片機器人磨削末端顫振檢測方法。該方法通過對原始信號進行排列熵特征提取，并利用該特征進行顫振檢測，當(dāng)排列熵值低于0.95時，判定機器人末端正在產(chǎn)生顫振，立即向系統(tǒng)發(fā)出警告，提醒系統(tǒng)采取顫振抑制措施。

針對實時性要求，采用基于序數(shù)模式的排列熵算法對顫振原始信號實時計算排列熵，該方法相較于傳統(tǒng)排列熵算法計算時間減小了一個數(shù)量級，當(dāng)信號長度為10000時，基于序數(shù)模式排列熵算法的計算時間約為0.25 s。

最后，試驗結(jié)果表明，通過排列熵特征提取算法的改進，基于排列熵的航空發(fā)動機葉片機器人磨削末端顫振檢測方法能夠在顫振爆發(fā)前約0.48 s檢測出顫振，可及時向系統(tǒng)發(fā)出警告，為后續(xù)采取措施提供了更多的反應(yīng)時間。

目前本文中排列熵特征閾值的應(yīng)用范圍有待進一步拓展，后續(xù)還需進一步研究狀態(tài)識別閾值的自適應(yīng)選擇標(biāo)準(zhǔn)。