劉麗娜，秦飛躍，2，宋欣鋼

(1.黃河交通學(xué)院機(jī)電工程學(xué)院，河南焦作 454950；2.河南理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院，河南焦作 454003)

0 前言

數(shù)控機(jī)床是目前最為熱門的加工設(shè)備，而銑頭是數(shù)控機(jī)床中最重要的設(shè)備之一。銑頭加工成本較高、生產(chǎn)周期較長(zhǎng)，因此在有限的環(huán)境和運(yùn)行數(shù)據(jù)下，建立銑頭的壽命預(yù)測(cè)模型具有重要的研究意義[1]。銑頭在運(yùn)行過(guò)程中因關(guān)鍵零部件磨損、結(jié)構(gòu)問(wèn)題以及本身質(zhì)量問(wèn)題，導(dǎo)致加工精度下降，進(jìn)而使得銑頭失效[2-3]。當(dāng)銑頭出現(xiàn)一定的磨損之后，最直接的表現(xiàn)就是溫度不斷上升，因此可以通過(guò)測(cè)定運(yùn)行過(guò)程中的溫升數(shù)據(jù)判斷銑頭的運(yùn)行狀態(tài)[4]。彭衛(wèi)文等[5]使用運(yùn)行過(guò)程中的溫升數(shù)據(jù)，應(yīng)用貝葉斯方法，結(jié)合故障機(jī)制構(gòu)建加速性能退化實(shí)驗(yàn)，建立了銑頭的可靠性評(píng)估預(yù)測(cè)模型。劉毅[6]分析了直角銑頭的常見(jiàn)故障，發(fā)現(xiàn)直角銑頭故障一般都是由臺(tái)階、螺旋傘齒輪齒面以及軸承磨損、殼體破壞等構(gòu)成，故障的主要原因還是銑頭內(nèi)部磨損。但是，目前有關(guān)銑頭的壽命預(yù)測(cè)模型研究相對(duì)較少，再者銑頭的壽命預(yù)測(cè)數(shù)據(jù)相對(duì)難以獲取，許多研究學(xué)者目前主要采用可靠性研究方法對(duì)數(shù)控機(jī)床及其關(guān)鍵零部件進(jìn)行壽命預(yù)測(cè)[7-8]。

銑頭運(yùn)行過(guò)程中的溫度升高是其內(nèi)部磨損、性能退化的直接體現(xiàn)。實(shí)際生產(chǎn)過(guò)程中，銑頭的退化狀態(tài)不可直接測(cè)量，本文作者提出基于狀態(tài)空間模型對(duì)銑頭運(yùn)行過(guò)程中的退化狀態(tài)進(jìn)行描述，并使用貝葉斯?fàn)顟B(tài)估計(jì)對(duì)建立的模型進(jìn)行求解，計(jì)算銑頭的失效時(shí)間分布，實(shí)現(xiàn)對(duì)銑頭實(shí)時(shí)壽命的預(yù)測(cè)[9]。

1 模型建立

在設(shè)備或者系統(tǒng)運(yùn)行過(guò)程中，會(huì)呈現(xiàn)出兩種狀態(tài)變量：一種是可以直接通過(guò)相應(yīng)的傳感器測(cè)量到的數(shù)據(jù)，這種狀態(tài)變量稱為觀測(cè)量，也可稱為可觀測(cè)序列；另外一種是無(wú)法通過(guò)相應(yīng)的傳感器測(cè)量，但是卻真實(shí)存在的變量，這種狀態(tài)變量稱為狀態(tài)向量。狀態(tài)向量能真實(shí)反映出設(shè)備或者系統(tǒng)運(yùn)行的內(nèi)在運(yùn)行狀態(tài)，它由一個(gè)馬爾科夫鏈構(gòu)成。狀態(tài)空間模型建立了可觀測(cè)量與狀態(tài)向量之間的關(guān)系，通過(guò)狀態(tài)空間模型，能將設(shè)備或者系統(tǒng)運(yùn)行過(guò)程中的內(nèi)在狀態(tài)呈現(xiàn)出來(lái)。

圖1描述了系統(tǒng)或者設(shè)備運(yùn)行過(guò)程中可觀測(cè)量與狀態(tài)向量之間的關(guān)系，其中A0，A1，…，An，An+1指的是系統(tǒng)或者設(shè)備運(yùn)行過(guò)程中的狀態(tài)向量，與之對(duì)應(yīng)的B0，B1，…，Bn，Bn+1是系統(tǒng)或者設(shè)備的可觀測(cè)序列。

圖1 狀態(tài)向量與可觀測(cè)序列之間的關(guān)系Fig.1 Relationship between state vector and observable sequence

貝葉斯?fàn)顟B(tài)估計(jì)與預(yù)測(cè)可以用來(lái)對(duì)狀態(tài)空間模型求解。每當(dāng)傳感器觀測(cè)到一個(gè)數(shù)據(jù)時(shí)，可以得到一個(gè)先驗(yàn)分布，通過(guò)貝葉斯公式計(jì)算得到一個(gè)后驗(yàn)分布。當(dāng)觀測(cè)到新數(shù)據(jù)時(shí)，之前產(chǎn)生的后驗(yàn)分布可以作為先驗(yàn)分布，繼續(xù)參與計(jì)算得到新的后驗(yàn)分布，如此循環(huán)往復(fù)，每當(dāng)觀測(cè)到新的數(shù)據(jù)，就可以及時(shí)更新模型參數(shù)，從而實(shí)現(xiàn)對(duì)設(shè)備或者系統(tǒng)的狀態(tài)估計(jì)。

銑頭作為一種復(fù)雜的精密機(jī)械設(shè)備，隨著運(yùn)行時(shí)間的增加，會(huì)因?yàn)殛P(guān)鍵零部件磨損而使得加工精度降低，而關(guān)鍵的零部件磨損無(wú)法通過(guò)相應(yīng)的傳感器測(cè)量出來(lái)，最終表現(xiàn)出來(lái)的是銑頭溫度增加，因此可以通過(guò)測(cè)量溫度，利用建立的狀態(tài)空間模型實(shí)現(xiàn)對(duì)銑頭的壽命預(yù)測(cè)。

銑頭性能退化及實(shí)時(shí)壽命的狀態(tài)空間模型可以通過(guò)式(1)來(lái)實(shí)現(xiàn)[10]：

(1)

其中：gθ(·)為狀態(tài)方程，描述銑頭退化狀態(tài)隨著運(yùn)行時(shí)間的變化關(guān)系；fθ(·)為觀測(cè)方程，指銑頭在運(yùn)行過(guò)程中溫度與其內(nèi)部的關(guān)系；t為銑頭的運(yùn)行時(shí)間；yt為銑頭運(yùn)行到t時(shí)刻的可觀測(cè)序列，即溫度；xt為銑頭運(yùn)行到t時(shí)刻的實(shí)際退化狀態(tài)；wt和vt分別為系統(tǒng)的過(guò)程噪聲和觀測(cè)噪聲。

已知t時(shí)刻的銑頭狀態(tài)后驗(yàn)分布π(xt|y1，t)，根據(jù)貝葉斯公式，可得到任意t+k的狀態(tài)分布：

π(xt+k|y1,n)=

(2)

累積失效概率指設(shè)備或者系統(tǒng)運(yùn)行期間內(nèi)的失效概率。假設(shè)文中研究對(duì)象銑頭設(shè)定的溫度閾值為TR，當(dāng)銑頭運(yùn)行過(guò)程中溫度Tr超過(guò)設(shè)定的閾值時(shí)，即可判定銑頭失效。設(shè)F(t)為銑頭運(yùn)行t時(shí)刻后的累積失效概率，累積失效概率計(jì)算公式[11]為

(3)

利用式(2)得到t+k時(shí)刻的狀態(tài)分布后，根據(jù)式(3)可以得到t+k時(shí)刻的銑頭累積失效概率為

(4)

失效概率密度是對(duì)累積失效概率公式的微分，銑頭在t時(shí)刻的失效概率密度的表達(dá)式為

本次研究納入的4100例孕婦中共計(jì)有26例經(jīng)產(chǎn)后或者引產(chǎn)證實(shí)為胎兒肢體缺陷。產(chǎn)前常規(guī)超聲檢出12例，聯(lián)合檢查檢出17例，連續(xù)順序追蹤超聲檢出25例。連續(xù)順序追蹤超聲檢出效果明顯優(yōu)于其他超聲檢查，P＜0.05，差異具有統(tǒng)計(jì)學(xué)意義，見(jiàn)表1。

(5)

給定銑頭的溫度后，以動(dòng)態(tài)線性增長(zhǎng)模型作為狀態(tài)方程，描述銑頭的退化軌跡[12]，如式(6)所示：

(6)

使用R語(yǔ)言對(duì)建立的銑頭性能退化模型進(jìn)行編程并求解[13]，可以得到任意時(shí)刻陶瓷銑頭的性能退化狀態(tài)、累積失效概率以及失效概率密度。

任何預(yù)測(cè)模型得出的結(jié)果與真實(shí)值都有一定的偏差，在預(yù)測(cè)過(guò)程中也有一定的不確定性，為了消除此種不確定性帶來(lái)的影響，在預(yù)測(cè)模型中引入置信區(qū)間來(lái)增加預(yù)測(cè)結(jié)果的可靠性，文中將置信區(qū)間設(shè)置為95%。

2 結(jié)果與分析

2.1 銑頭退化狀態(tài)估計(jì)與預(yù)測(cè)

由于壽命數(shù)據(jù)具有難以獲取、樣本小以及成本高等特點(diǎn)，因此文中使用的銑頭在運(yùn)行過(guò)程中的溫度數(shù)據(jù)及其獲取方法見(jiàn)文獻(xiàn)[5]。圖2所示為銑頭運(yùn)行過(guò)程中的溫升數(shù)據(jù)。

圖2 銑頭溫升數(shù)據(jù)Fig.2 Milling head temperature rise data

通過(guò)對(duì)銑頭溫升的數(shù)據(jù)進(jìn)行最小二乘擬合，可以得到銑頭壽命公式如式(7)所示：

T(t)=0.408 8×t0.593 9

(7)

其中：t為銑頭實(shí)際運(yùn)行時(shí)間，h；T(t)為銑頭溫升，℃。

假設(shè)當(dāng)溫升達(dá)到13 ℃時(shí)，銑頭壽命結(jié)束。通過(guò)式(7)計(jì)算此時(shí)銑頭理論運(yùn)行時(shí)間為338.6 h。

圖3所示為基于銑頭運(yùn)行前110 h的溫升數(shù)據(jù)得到的銑頭預(yù)測(cè)運(yùn)行時(shí)間。通過(guò)上述建立的狀態(tài)空間模型計(jì)算出銑頭溫升為13 ℃時(shí)所需要的運(yùn)行時(shí)間為191.8 h，與通過(guò)式(7)計(jì)算出的理論運(yùn)行時(shí)間338.6 h差別較大。銑頭運(yùn)行前110 h時(shí)退化狀態(tài)變化率δ110=0.023。圖中陰影面積為根據(jù)式(3)(4)、使用前110 h的溫升數(shù)據(jù)計(jì)算得到的運(yùn)行第220 h時(shí)銑頭的累積失效概率。此時(shí)的陰影面積較小，說(shuō)明銑頭在較短的運(yùn)行時(shí)間時(shí)，失效概率較低，銑頭可靠度較高。

圖3 基于運(yùn)行前110 h銑頭溫升數(shù)據(jù)的退化狀態(tài)估計(jì)與預(yù)測(cè)Fig.3 Degradation state estimation and prediction based on the first 110 h milling head temperature rise data

圖4所示為基于銑頭運(yùn)行前220 h的溫升數(shù)據(jù)得到的銑頭預(yù)測(cè)運(yùn)行時(shí)間。通過(guò)狀態(tài)空間模型計(jì)算出銑頭溫升為13 ℃時(shí)所需要的運(yùn)行時(shí)間為275.3 h，與通過(guò)式(7)計(jì)算出的理論運(yùn)行時(shí)間338.6 h還是有一定的差別，但是與運(yùn)行前110 h溫升數(shù)據(jù)相比，差別進(jìn)一步縮小。銑頭運(yùn)行前220 h時(shí)退化狀態(tài)變化率δ110=0.015，退化狀態(tài)變化率與運(yùn)行前110 h相比，進(jìn)一步減小。圖中陰影面積為根據(jù)式(3)(4)、使用前220 h的溫升數(shù)據(jù)計(jì)算得到的運(yùn)行第310 h時(shí)銑頭的累積失效概率。此時(shí)與前110 h溫升數(shù)據(jù)對(duì)比，陰影面積增大，說(shuō)明銑頭在又運(yùn)行了110 h后，失效概率增大，可靠度降低。

圖4 基于運(yùn)行前220 h銑頭溫升數(shù)據(jù)的退化狀態(tài)估計(jì)與預(yù)測(cè)Fig.4 Degradation state estimation and prediction based on the first 220 h milling head temperature rise data

圖5所示為基于銑頭運(yùn)行前310 h的溫升數(shù)據(jù)得到的銑頭預(yù)測(cè)運(yùn)行時(shí)間。通過(guò)狀態(tài)空間模型計(jì)算出銑頭溫升為13 ℃時(shí)所需要的運(yùn)行時(shí)間為330.5 h，與式(7)計(jì)算出的理論運(yùn)行時(shí)間338.6 h差別進(jìn)一步縮小。可以看出：在置信區(qū)間內(nèi)，失效時(shí)間預(yù)測(cè)已經(jīng)超過(guò)理論運(yùn)行時(shí)間，證明預(yù)測(cè)模型的準(zhǔn)確性。由于失效時(shí)間預(yù)測(cè)已經(jīng)超過(guò)理論運(yùn)行時(shí)間，因此再對(duì)運(yùn)行338.6 h后銑頭的累積失效概率計(jì)算已經(jīng)沒(méi)有實(shí)際意義，這里不再給出累積失效概率預(yù)測(cè)圖。銑頭運(yùn)行前310 h時(shí)退化狀態(tài)變化率δ110=0.008，其退化狀態(tài)變化率進(jìn)一步減小，這也能與圖2銑頭溫升數(shù)據(jù)前期增加較快、后期增加較慢相印證。

圖5 基于運(yùn)行前310 h銑頭溫升數(shù)據(jù)的退化狀態(tài)估計(jì)與預(yù)測(cè)Fig.5 Degradation state estimation and prediction based on the first 310 h milling head temperature rise data

2.2 銑頭失效概率密度預(yù)測(cè)

結(jié)合式(4)(5)計(jì)算可以得到任意時(shí)刻銑頭的失效概率密度分布。文中分別選取銑頭運(yùn)行前110、220、310 h時(shí)的溫升數(shù)據(jù)，計(jì)算對(duì)應(yīng)的失效概率密度分布，如圖6所示。其中豎直垂線是銑頭溫升閾值為13 ℃時(shí)，由式(7)計(jì)算得到的理論運(yùn)行時(shí)間，為338.6 h。

圖6 銑頭失效概率密度預(yù)測(cè)Fig.6 Failure-probability density prediction of milling head

由圖6可以看出：隨著溫升數(shù)據(jù)不斷增多，銑頭失效時(shí)間預(yù)測(cè)越來(lái)越收斂于理論值，說(shuō)明模型預(yù)測(cè)的不確定性越來(lái)越小。這是因?yàn)樨惾~斯能不斷融入新的溫升數(shù)據(jù)，通過(guò)計(jì)算先驗(yàn)分布和后驗(yàn)分布來(lái)及時(shí)更新?tīng)顟B(tài)空間模型參數(shù)。

3 結(jié)論

銑頭在運(yùn)行過(guò)程中，由于關(guān)鍵零部件磨損使得加工精度降低，從而導(dǎo)致其性能狀態(tài)不斷退化，而且退化狀態(tài)不可直接測(cè)量，只能通過(guò)測(cè)定運(yùn)行過(guò)程中銑頭溫升來(lái)判斷退化狀態(tài)；同時(shí)銑頭壽命數(shù)據(jù)具有難獲取、樣本小以及成本高等特點(diǎn)，使得難以實(shí)時(shí)預(yù)測(cè)銑頭壽命。文中以實(shí)際運(yùn)行過(guò)程中銑頭溫升數(shù)據(jù)為基礎(chǔ)，結(jié)合狀態(tài)空間模型與貝葉斯理論，建立了銑頭的實(shí)時(shí)壽命預(yù)測(cè)模型。得到的結(jié)論如下：

(1)建立的銑頭實(shí)時(shí)壽命預(yù)測(cè)模型可以根據(jù)前期銑頭的溫升數(shù)據(jù)，來(lái)預(yù)測(cè)任意時(shí)刻銑頭退化狀態(tài)以及剩余壽命，并且可以根據(jù)最新融合的溫升數(shù)據(jù)及時(shí)更新模型參數(shù)，從而使得預(yù)測(cè)結(jié)果越來(lái)越收斂于理論值。

(2)隨著運(yùn)行時(shí)間的增加，預(yù)測(cè)模型計(jì)算得到的銑頭退化狀態(tài)變化率逐漸減小，這與銑頭剛開(kāi)始投入使用時(shí)溫升增加較快、后期增加較慢這一現(xiàn)象相互印證。

(3)不同運(yùn)行時(shí)間內(nèi)，模型計(jì)算出的銑頭失效概率密度隨運(yùn)行時(shí)間的增加不斷接近真實(shí)值，一定程度上可以為銑頭失效做出判斷。