基于改進有序聚類法的立式加工中心進給系統溫測點優化

李傳珍，李國龍，陶小會，龐 源

（重慶大學機械傳動國家重點實驗室，重慶400044）

隨著智能制造的高速發展，智能化、高精度、高速化的數控機床已成為研發重點［1］。研究表明，在精密加工過程中因熱變形引起的精度誤差占總誤差的40%～70%，因此熱誤差補償是提高機床加工精度的關鍵［2-3］。熱誤差補償的有效性受熱誤差預測模型精度的影響，目前建立精確熱誤差預測模型的難點之一是準確選取溫測點。

針對溫測點優化問題，國內外學者已通過有限元法［4］、神經網絡［5-6］、AVQ（adaptive vector quantization，自適應矢量量化）聚類［7］、模糊聚類法［8］等實現了溫測點優化，有效減少了溫測點的數目，提高了熱誤差預測模型的精度和魯棒性。但是，隨著機床加工精度的提高，對熱誤差預測模型的精度要求也更高，基于傳統單一理論的溫測點優化方法已不能滿足要求。為選取更加準確的溫測點，張偉等［9］采用灰色關聯分析和模糊聚類法實現了臥式加工中心溫測點的優化；Miao等［10］利用模糊聚類法與灰色關聯度分析獲得了各類溫測點的關鍵點，并采用逐步回歸法剔除了非顯著溫測點，完成了高速加工中心主軸溫度測點的優化；叢明等［11］提出了將簡單相關分析、模糊聚類法、灰色關聯度分析和熱誤差預測模型參數統計學分析相結合的溫測點優化方法，并基于這種優化方法確定了臥式加工中心主軸的最佳溫測點；李艷等［12］采用互信息法建立了溫度變量與熱誤差的綜合關聯度矩陣，根據改進模糊聚類法對溫測點進行分類，并結合綜合關聯度矩陣選取數控磨齒機主軸的關鍵溫測點；李逢春等［13］提出一種兼顧歐式距離和相關系數的改進系統聚類法，有效解決了易干擾溫度信號區分問題和溫測點間共線問題，實現了重型落地銑鏜床的溫測點優化；魏弦等［14］采用模糊聚類法和相關性分析確定了龍門機床進給系統的關鍵溫測點，并采用特征提取法獲得熱誤差預測模型的自變量，有效消除了溫測點間的復共線性，確定了最佳溫測點。綜合來看，目前基于多種理論的溫測點優化方法主要涉及溫測點共線性消除和溫測點聚類等，其中：溫測點共線性消除普遍使用逐步回歸法和簡單相關性分析等全局式篩選方法，但這些方法均未考慮溫測點布置的疏密情況，對所有溫度變量籠統地進行篩選，當溫測點布置密集時，上述方法的篩選效果較差；溫測點聚類大多采用模糊聚類法及其相關理論，但模糊聚類分析的閾值大多通過試驗和經驗獲得，且在分類時只根據溫測點溫度樣本的親疏性質或距離而未考慮溫度有序傳遞的特點，這會導致溫測點優化后所建立的熱誤差預測模型的精度和魯棒性較差?？紤]到立式加工中心進給系統的有效進給行程較短，發熱源較集中且溫度有序傳遞，相較于其他機床，該加工中心具有不同的溫度場分布，其溫測點更密集，且溫測點間的耦合效應也更顯著，因此其溫測點優化方法有待深入研究。

針對上述問題，筆者以三軸立式加工中心為研究對象，基于該機床進給系統的發熱特點，提出一種基于改進有序聚類法的溫測點優化方法。首先，該方法改變了傳統溫測點全局式初篩的方式，采用互信息值（即相關性）對機床進給系統各發熱部件的多個溫測點進行局部式初篩，不僅提高了溫測點初篩精度，還有效消除了測點間的共線性。其次，基于機床進給系統溫度有序傳遞的特點，對初篩后的溫測點進行結合統計學檢驗的改進有序聚類分析，從而確定最佳溫測點。最后，根據最佳溫測點建立立式加工中心進給系統的熱誤差預測模型，并與基于傳統有序聚類法和灰色關聯度模糊聚類法建立的熱誤差預測模型進行比較，以驗證利用本文方法建立的熱誤差模型具有更高的精度和魯棒性。

1 基于改進有序聚類法的溫測點優化方法

1.1 改進有序聚類法的提出

傳統有序聚類法是：先通過計算類直徑矩陣和誤差函數，獲得溫測點分類；再計算誤差函數曲線中變化較明顯處相鄰聚類數誤差函數的比值，比值較大且通過顯著性檢驗的聚類數為最佳聚類數；最后通過溫度變量與熱誤差之間的簡單相關性分析，選出關鍵溫測點［15］。然而，相鄰聚類數誤差函數的比值僅反映了2個聚類點誤差函數的大小關系，并不能準確反映誤差函數曲線的變化值，且有序聚類法的聚類數一般在12個以下，較少的樣本點受擬合曲線類型的影響較大，因此無法根據誤差函數曲線斜率變化率確定最佳分類。此外，為消除溫測點間的耦合性和共線性，一般先根據溫度變量與熱變形的相關性初步篩選溫測點，再進行分類，這種全局式溫測點初篩方法在溫測點布置稀疏時效果較好，但當溫測點布置密集時，此方法篩選效果較差，無法精確篩選出關鍵溫測點。

立式加工中心屬中小型數控機床，結構緊湊，進給行程較短，進給系統發熱源集中且溫度場分布有序；在熱誤差試驗中溫測點布置得十分密集。因此，結合立式加工中心進給系統的發熱特點，根據互信息值針對性地初步篩選進給系統各發熱部件不同位置的同類溫測點，并針對原有序聚類法最佳聚類數難以確定的問題，提出一種基于多元線性回歸，由F檢驗、t檢驗、判定系數相關性檢驗和殘差比較等統計學分析確定最佳聚類數的改進有序聚類法，并利用該方法確定立式加工中心進給系統最佳溫測點。

1.2 改進有序聚類法的步驟

針對布置了多個溫測點的發熱部件，先由互信息值初步篩選各溫測點，計算初選溫度變量的類直徑矩陣和比較誤差函數，獲得初始分類；再結合統計學分析確定最佳聚類數和最佳溫測點。具體步驟如下所述。

1）計算互信息值。

根據信息論產生的互信息法用信息熵來衡量2個序列變量的相關性［12］。以溫測點溫度變化量為溫度樣本，設立式加工中心進給系統溫測點溫度隨時間的變化量ΔTi(t)=Ti(t)-T（0i表示溫測點編號，i= 1,2,…,n），每個溫測點有m個溫度變量樣本，則各溫測點溫度變量樣本集合ΔTi={ΔTij(t)}(j= 1,2,…,m)，其中ΔTij(t)表示第i個溫測點的第j個溫度變化量。絲杠熱誤差隨時間的變化量Y(t)為該時刻絲杠測量終點熱誤差與絲杠測量起始點熱誤差之差，即Y(t)=Y終(t)-Y0(t)，熱誤差樣本集合Y=Yj(t)(j= 1,2,…,m)，則各溫測點溫度變量ΔTi與熱誤差Y相關性的互信息值為：

式中：p(?)為概率密度函數，p( ΔTij,Yj)是ΔTi與Y的聯合概率密度函數。

根據式（1）可計算得各溫測點溫度變量與熱誤差的互信息值I，I越大表示該溫測點與熱變形的相關性越大。根據相關性大小，對各部件不同位置的同類溫測點（如在螺母上下端面、外殼、內圈布置4個溫測點是為了更準確地獲得螺母的發熱位置，則這4個測點可歸為螺母的同類測點）進行初步篩選，消除溫測點間的耦合性。

2）數據規格化處理。

為消除不同數據類型間的差異，降低計算誤差，將溫測點溫度變量樣本矩陣X中的元素xab變換為zab，獲得有序樣本矩陣Z，變換公式為：

式中：xab表示溫測點初篩后第b個原始樣本的第a個溫測點溫度變化值；a= 1,2,…,l；b= 1,2,…,m。

3）計算類直徑矩陣D。

有序聚類法采用直徑表示類中的差異程度，類中差異越小，直徑就越小，類中樣本就越集中［15］。設有序聚類樣本中的一個類為{za,za+1,…,zc} (a＜c≤l)，用D(za,zc)表示該類的類直徑，D(za,zc)可用離差平方和表示為：

式中：za為正規化處理后第a個溫測點的溫度變化值，為第a個到第c個溫測點溫度變化值的均值，即

4）計算誤差函數。

將l個溫測點分成k類，各類的誤差函數可表示為各類的直徑之和，將溫度樣本簡記為f，q(k,l)為l個溫測點分成k類的一種分法，則誤差函數ψ[q(k,l)]為：

式中：D(fa,fa- 1 )表示某一子類的離差平方和。

當n和l固定時，ψ[q(k,l)]越小，表示各類的類直徑越小，分類越合理［15］。因此，有序聚類法的最優分類就是選擇一種分法使得ψ[q(k,l)]最小。

5）確定最優分法。

當k=2時，q( 2，l)是所有可能的分類中使誤差函數最小的分法［15］，故有：

對上式進行歸納，即有：

首先尋找合適的a值，使得ψ[q(k,l)]=ψ[q(k-1，fk-1)]+D(fk，fl) 最 小，得 到 第k類Gk={fk，fk+1，…，fl}；接著尋找fk-1，使得ψ[q(k-1，fk-1)]最小，其中

ψ[q(k-1，fk-1)]=ψ[q(k-2，fk-1-1)]+D(fk-1，fk-1)，得到第k-1 類Gk-1={fk-1，fk-1+1，…，fk-1}；依此類推可得到所有的分類G1，G2，…，Gk。

6）確定最佳聚類數。

首先繪制誤差函數隨聚類數k變化的誤差函數曲線，選取該曲線明顯變化處的聚類數作為溫測點數目，然后結合各溫測點的互信息值選取建模所用的溫測點，建立包含多個不同溫測點的熱誤差多元線性回歸模型。對模型先進行F檢驗，判斷熱誤差模型的顯著性；再分析其線性相關性系數R2和最大殘差，R2越大，表明該模型的擬合效果越好，最大殘差值越小，表明該模型的精度越高；最后對回歸系數進行顯著度t檢驗，剔除不顯著的溫測點，進而確定最佳聚類數和最佳溫測點。

1.3 基于改進有序聚類法的溫測點優化流程

根據采集的機床進給系統溫測點的溫度和熱變形數據，基于改進有序聚類法確定最佳溫測點，總體流程如圖1所示。

圖1 基于改進有序聚類法的溫測點優化流程Fig.1 Optimizing process of temperature measurement points based on improved sequential clustering method

2 溫測點優化試驗

2.1 立式加工中心進給系統熱源分析

不同數控機床因結構、工作轉速（進給速度）、有效進給行程存在差異，其溫度場分布不同，則熱誤差建模選取的關鍵溫測點也不同［2］。文獻［14，16-17］分別對龍門加工中心、鏜床和立式加工中心的進給系統進行研究后發現：龍門加工中心各熱源溫度曲線上升趨勢不同且熱平衡時各熱源溫度差較大，鏜床各熱源溫度曲線上升趨勢相似但熱平衡時各熱源溫度差明顯，立式加工中心各熱源溫度曲線分布密集、上升趨勢相似且熱平衡時各熱源溫度差較小，由此可見不同數控機床發熱規律差別較大，溫測點的布置也不同。因此，考慮到立式加工中心進給系統結構緊湊，進給行程較短，發熱源集中且溫度場分布有序的特點，根據文獻［17-18］中的立式加工中心進給系統熱源分析，確定其溫測點布置位置為電機、聯軸器、上下軸承蓋和螺母內外圈。

2.2 溫度與熱變形樣本采集試驗

本文的試驗對象為華中八型立式加工中心VMC850，采集其進給系統（以Z向為例）的溫測點溫度與熱變形數據。根據2.1節的熱源分析，在該立式加工中心Z向進給系統的電機、聯軸器、上下軸承蓋和螺母內外圈處共布置13 個溫測點。將貼片式PT100 溫度傳感器安裝在熱源位置處，HIO-1075 溫度采集板采集溫度信號并將它傳至數控系統，利用華中數控專用數據采集軟件SSTT實現在PC（personal computer，個人計算機）端記錄溫度數據并儲存，各傳感器安裝位置如表1所示，13個溫測點的分布如圖2所示。如圖2（b）所示，用RenishawXL-80 激光干涉儀S采集機床Z向進給系統的熱變形數據，激光干涉儀入射鏡與折射鏡安裝在工作臺上，反射鏡安裝在主軸上，激光發射器安裝于機床正前方，PC與發射器相連接。激光干涉儀只能間歇測量，試驗機床Z向進給系統的效行程為600 mm，每隔50 m 設定一個測點，共13個測點；采用正反雙向測量，消除絲杠反向間隙的影響，激光干涉儀測點分布如圖3所示。

表1 傳感器安裝位置Table 1 Installation position of sensors

圖2 Z向進給系統溫測點分布Fig.2 Layout of temperature measurement points for Z-direction feed system

圖3 激光干涉儀測點分布Fig.3 Layout of measurement points of laser interferometer

本文擬研究機床從冷機到熱平衡階段的熱誤差變化規律，采集進給速度為8 000，12 000 和16 000 mm/min時機床熱機階段的溫度與熱誤差數據，共采集48 組熱誤差數據。以進給速度v=8 000 mm/min為例，機床開機后，主軸保持靜止，Z向進給軸以8 000 mm/min的進給速度往復運動，激光干涉儀采樣周期為15 min，整個試驗時長為240 min。濾波平穩處理后，13 個溫測點處溫度隨時間的變化曲線如圖4 所示，由圖可知，t=0—150 min 階段溫度上升迅速，在t=150 min 左右基本達到熱平衡，t=150—180 min 階段溫度略微上升，t=180 min 后溫度基本不再變化。t=0，15，30，45，60，75，90，105，120，135，150 min時熱誤差曲線如圖5所示。

由于各傳感器初始溫度并不是環境溫度，在安裝激光干涉儀時，需要上下移動主軸使激光干涉儀的反射光與折射光重合，不可避免地會造成機床部件升溫。為保證熱誤差預測模型的精度，采用溫度變化量作為溫度變量樣本，即ΔTi(t)=Ti(t)-T0。圖5中的

圖4 溫測點溫度隨時間的變化曲線（v=8 000 mm/min）Fig.4 Variation curve of temperature of temperature measurement points over time（v=8 000 mm/min）

圖5 溫升階段Z向進給系統熱誤差曲線（v=8 000 mm/min）Fig.5 Thermal error curve of Z-direction feed system in temperature rising stage（v=8 000 mm/min）

熱誤差值以絲杠收縮方向為正，伸長方向為負，因試驗時間為冬季，初始測量時絲杠處于收縮狀態。

2.3 溫測點優化結果

根據1.2 節中互信息值計算原理，以t=0—150 min 階段機床進給系統各溫測點溫度變化量ΔTi(t)和熱誤差變化量Y(t)作為樣本。但由于進給系統熱誤差值與進給軸位置有關，為更準確獲取溫度變化量與熱變形間的相關性，考慮到圖5中的熱誤差與進給位置幾乎呈線性關系，采用熱誤差曲線一次擬合斜率Kj(t)(j=1，2，…，13)代替Y(t)進行互信息值計算。將ΔT、Kj(t)代入式（1），計算得到各溫測點的互信息值，如表2所示。

由上文可知，在進給軸上軸承、螺母和下軸承等重要位置均布置了多個溫測點，它們均能反映相應位置的發熱情況，因此，為消除溫測點間的耦合性，結合互信息值計算結果，選擇各位置處相關性較高的2個溫測點作為主要影響因素，如進給軸上軸處溫測點4，5，螺母處溫測點7，9，下軸承處溫測點10，12，最后得到1，2，4，5，7，9，10，12，13等9個溫測點作為主要影響因素。

表2 各溫測點的互信息值Table 2 Mutual information value of temperature measurement points

根據改進有序聚類分類法，對上述9個溫測點在t=0—150 min 階段的溫度樣本ΔTi(t)升序排列并進行規格化處理，計算直徑矩陣D和最小誤差函數ψ[q(k，l)]。得到的有序聚類結果如表3所示，誤差函數值與聚類數的關系曲線如圖6所示。

表3 溫測點有序聚類結果Table 3 Sequential clustering results of temperature measurement points

圖6 溫測點誤差函數值與聚類數的關系曲線Fig.6 Relation curve of error function value and number of clusters of temperature measurement points

由圖6可以看出，誤差曲線在聚類數為3，4，5時均有明顯轉折，對應建立包含3個測點、4個測點和5個測點的線性回歸模型，多個測點屬于同一聚類時，選擇互信息值大的測點（如分為3類時，測點2，4，5，7，9，10，12屬于一類，選取互信息值最大的測點4，此時包含3個測點的熱誤差多元線性回歸模型的輸入樣本為測點1，4，13的溫度變化值）建立熱誤差模型。由圖5 的熱誤差曲線可以看出，熱誤差與溫度（時間）、進給軸位置有關，溫度一定時熱誤差與進給軸位置近似呈線性變化，因此用E=tanβ(Px-P0)來表示熱誤差與進給軸位置的關系［19］，其中：Px為機床進給軸實時位置，P0為機床進給軸相對零點，tanβ為熱誤差擬合直線斜率?；诙嘣€性回歸建立立式加工中心Z向進給系統溫度變化量與熱誤差曲線一次擬合直線斜率的回歸模型，并對該模型進行統計學分析，結果如表4所示。

分析表4可得如下結論：

1）對不同測點的熱誤差模型進行F檢驗時，各模型的F值均在相應的拒絕域且概率p（即顯著度Sig.）都小于顯著性水平0.05，表明該顯著水平下各溫測點的溫度變化量對熱誤差的線性影響顯著。

2）判定系數R2、標準殘差和最大殘差直觀反映了熱誤差模型的擬合效果與預測精度，R2越大，表明擬合擬合效果越好，標準殘差和最大殘差越小，表明精度越高。從表中可以看出對于包含3個測點和4個測點的熱誤差模型，隨著測點的增多，R2增大，標準殘差與最大殘差減小，但對于包含4個測點和5個測點的熱誤差模型，即使測點增多，但R2、標準殘差和最大殘差均不變，說明測點增多雖然會提高熱誤差模型的精度，但精度的提高卻越來越不明顯。

3）對各模型的偏相關系數進行t檢驗時，包含5個測點的熱誤差模型的ΔT12沒有通過t檢驗，且從相關系數可以看出，ΔT12的貢獻率很小。

綜上，綜合考慮模型精度與魯棒性要求，選取4個測點1，4，9，13作為最佳溫測點。

3 溫測點優化結果驗證

根據改進有序聚類法得到的最佳溫測點，基于多元線性回歸建模方法得到機床Z向進給系統的熱誤差預測模型為：

其中P0由文獻［19］中的均值法計算得到，P0=17.98。

采用傳統有序聚類法選取誤差函數曲線非負斜率α(k)最大時的k值為較優聚類數［15］，即有α(3)＞α(4)＞α(5)，建立包含3個測點的熱誤差預測模型為：

采用灰色關聯度模糊聚類法［10］進行溫測點優化，并建立包含4個測點的熱誤差預測模型為：

表4 溫測點的熱誤差模型統計學分析結果Table 4 Statistical analysis results of thermal error model of temperature measurement points

利用上述3種熱誤差預測模型分別對進給速度v=8 000 mm/min的溫升階段中t=60 min與進給速度v=12 000 mm/min的熱平衡階段中t=150 min時的熱誤差進行預測，基于3種模型的熱誤差預測曲線與殘差曲線對比分別如圖7和圖8所示。激光干涉儀的反向位置（600 mm）處因儀器誤差影響測量結果，因此觀測0～550 mm 位置處各熱誤差模型的預測結果。從圖7和圖8發現，相較于利用傳統有序聚類法和灰色關聯度模糊聚類法得到的熱誤差預測曲線，利用改進有序聚類法得到的熱誤差預測曲線在升溫和熱平衡階段都與實際測量曲線更加接近，殘差變化幅度較小。

為進一步驗證3種熱誤差模型的預測效果，分析了溫升階段與熱平衡階段各模型的均方根誤差、最大殘差和平均殘差，結果如表5和表6所示。表中均方根誤差用來描述誤差的總體離散程度，殘差表示回歸模型的擬合值與實測值之間的誤差［13］。

圖7 溫升階段Z 向進給系統熱誤差模型預測結果(v=8 000 mm/min,t=60 min)Fig.7 Prediction results of thermal error model of Z-direction feed system in temperature rising stage (v=8 000 mm/min,t=60 min)

圖8 Z 向進給系統熱平衡階段熱誤差模型預測結果（v=12000 mm/min，t=150 min）Fig.8 Prediction results of thermal error model in temperature stage(v=12 000 mm/min,t=150 min)

表5 溫升階段Z 向進給系統熱誤差模型預測效果(v=8 000 mm/min,t=60 min)Table 5 Prediction effect of thermal error model of Zdirection feed system in temperature rising stage (v=8 000 mm/min, t=60 min) 單位：μm

表6 熱平衡階段Z向進給系統熱誤差模型預測效果(v=12 000 mm/min, t=150 min)Table 6 Prediction effect of thermal error model in temperature balance stage (v=12 000 mm/min,t=150 min) 單位：μm

從表5和表6可以看出，采用改進有序聚類法所建的熱誤差預測模型在溫升階段和平衡階段的均方根誤差低于1.05 μm，最大殘差僅為2.09和1.71 μm，平均殘差最小，低于1 μm，且熱平衡階段熱誤差模型預測精度高于溫升階段?；诟倪M有序聚類法的熱誤差模型在熱平衡階段的精度相較于另外2種模型分別提高了65.7%和39.8%。

4 結 論

1）根據立式加工中心進給系統的發熱特點，由互信息值初步篩選各發熱部件不同位置的同類溫測點，并基于改進有序聚類法選取最佳溫測點，有效解決了多溫測點間共線性的問題，減少了溫測點數量，提高了擬合精度和魯棒性。

2）結合統計學分析與有序聚類法，從顯著性檢驗（F檢驗和t檢驗）、擬合效果和殘差分析等方面確定最佳聚類數與最佳溫測點。相較于傳統有序聚類法，該方法直觀且全面地分析了溫測點數目與回歸模型精度的關系，選取了更精確的最佳聚類數。

3）建模實例表明：基于改進有序聚類法建立的熱誤差模型的均方根誤差、最大殘差、平均殘差分別小于1.05，2.1，1 μm；相較于采用傳統有序聚類法和灰色關聯度模糊聚類法建立的熱誤差預測模型，該熱誤差模型在熱平衡階段的精度分別提高了65.7%和39.8%。