基于改進多項式響應面的VPMCD方法及其在滾動軸承故障診斷中的應用

楊 宇，潘海洋，李 杰，程軍圣

（湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082）

滾動軸承的故障診斷本質上是一個模式識別的過程。在滾動軸承的故障模式識別中，神經網絡、支持向量機等模式識別方法得到了廣泛的應用。但是它們都存在著一定的局限性，如人工神經網絡具有收斂速度慢、網絡結構難以確定［1］等問題，支持向量機具有核函數及核參數難以確定［2］等問題。除此之外，神經網絡和支持向量機在進行模式識別時都忽略了從原始數據中所提取的特征值之間的相互內在關系。

然而，在機械故障診斷中，所有或部分特征值之間大都具有一定的內在關系，而且這種內在關系在不同的系統或類別（相同的系統在不同的工作狀態下）間具有明顯的不同。因此，可以對各個特征值之間的相互內在關系建立數學模型，對于不同的類別可以得到不同的數學模型，從而可以采用這些數學模型對被測試樣本的特征值進行預測，把預測結果作為分類的依據，進一步進行模式識別。基于此，Raghuraj與Lakshminarayanan提出了一種新的模式識別方法——基于變量預測模型的模式識別（Variable Predictive Model Based Class Discriminate，簡稱VPMCD）方法，同時還將該方法與神經網絡、支持向量機等其它模式識別方法進行了對比，結果驗證了VPMCD方法的有效性和優越性［3］。

在VPMCD方法中，模型的訓練過程也就是變量預測模型的建立過程，在建模過程中，變量預測模型的擬合精度影響著整個模式識別的效果。因此，獲得合理的變量預測模型至關重要。VPMCD是通過特征值之間存在相互內在關系來建立變量預測模型的，但是這種相互關系的具體情況卻難以確定，而且特征值之間相互內在關系的實際預測模型也無法得到。因此，該方法選取了線性（L）模型、線性交互（LI）模型、二次（Q）模型、二次交互（QI）模型四種模型作為實際模型的代理模型，并從中選取最佳代理模型，這四種模型都屬于多項式響應面（Polynomial Response Surface，簡稱PRS）模型［4］。多項式響應面模型是一種常用的響應面（Response Surface Method，簡稱 RSM）［5］模型，也是應用最廣泛的近似模型，該方法計算簡單，但其高階計算量卻很大，且提高的精度有限。另外，響應面法不能夠隨著樣本容量的增大而有效提高其近似精度，這兩點嚴重限制了該方法的使用［6］。究其原因，主要是PRS方法采用最小二乘擬合參數，而最小二乘擬合與真實數據往往存在偏差，并且放棄這些殘差項，從而導致模型近似過程中忽略了殘差項，而大多數情況下殘差項包含了很多重要信息，將其忽略容易使模型擬合誤差偏大。針對這一缺陷，本文將原方法中的多項式響應面法進行了改進，改進的多項式響應面法對PRS法忽略掉的殘差重新進行近似，保留了模型擬合的重要信息，因而比原方法的模型擬合效果更理想。基于此，本文提出了基于改進多項式響應面（Improved Polynomial Response Surface，簡稱IPRS）的VPMCD方法，并將其運用于滾動軸承故障診斷。

1 基于改進多項式響應面的VPMCD方法

1.1 VPMCD方法

以機械故障診斷問題為例，采用p個不同的特征值 X＝［X1，X2，…，Xp］來描述一個故障類別，對于其中的特征值Xi來說，當故障類別不同時，其他的一個或者多個特征值對其影響也會發生變化。因此，特征值Xi與其余的一個或者多個特征值之間存在著一定的函數關系，而這種關系可以是線性的，也可以是非線性的。為了識別滾動軸承的故障模式，需要有能夠描述這些函數關系的數學模型，以便對測試樣本的特征值進行預測，進一步對測試樣本進行分類，這種模型稱為變量預測模型（Variable predictive model，簡稱 VPM）。

為特征值Xi定義的變量預測模型是一個線性或非線性的回歸模型，可以選擇以下四種模型之一：

① 線性模型（L）：

② 線性交互模型（LI）：

③ 二次交互模型（QI）：

④ 二次模型（Q）：

式中：r≤p－1為模型階數。這四種模型都屬于多項式響應面（PRS）模型，以p個特征值為例，選取四種模型中任意一個模型，用特征值Xj（j≠i）對Xi進行預測，都可以得到：

式（5）稱為變量Xi的變量預測模型VPMi。其中，特征值Xi稱為被預測變量；Xj（j≠i）稱為預測變量；e為預測誤差；b0，bj，bjj，bjk為模型參數，可以通過所有訓練樣本的特征值對它們進行參數估計。

對于g類故障分類問題，提取p個特征值X＝［X1，X2，…，Xp］．在分別采用不同故障類別下的訓練樣本數據對預測模型進行訓練后，不同故障類別下的不同特征值就可以分別建立g×p個預測模型VPMki，其中k＝1，2，…，g代表不同的類別，i＝1，2，…，p代表不同的特征值。然后針對測試樣本，提取其特征值，并用g×p個模型VPMki分別對它們進行預測，得到g×p個預測值X～i，以同一類別下所有特征值的預測誤差平方和最小為判別函數，對測試樣本的故障類型和工作狀態進行分類，該方法稱為基于變量預測模型的模式識別（VPMCD）方法。

1.2 多項式響應面法

VPMCD在模型訓練過程中是從四種多項式響應面模型中選擇最佳模型為變量預測模型，在進行模型近似時主要是利用最小二乘擬合的方法，通過對數據的分析，找出相關因素與響應量之間的函數關系，并用這種函數代替原來的函數關系。PRS法通常是高階多項式，基本形式如下：

式中：xi是 m維自變量 x的第 i個分量；β0，βi和 βij是未知參數，將它們按順序構成列向量β，求解PRS模型最重要的就是求解向量β。將樣本點的值代入式（6）中，再利用最小二乘法可以求得其估計值，即 β＝

PRS方法采用的是泰勒展開式的思想，對于函數f（x），當采用常用的二次多項式時，其自變量x在x0處展開后得到：

式中：Δf（x0）為函數 f在 x0點的梯度與 Δx的點積；Hf（x0）為函數f在 x0點的 Hessian矩陣。PRS法根據樣本數據利用最小二乘法計算函數f在某個定義域區間內的梯度及Hessian矩陣，并將項o（‖x‖3）作為殘差而將其忽略。但是，該項往往會比較大，而且包含很多重要信息，忽略掉該項會使近似值和實際值之間存在較大差異。另外，由式（6）可知，如果一個PRS模型具有n個變量，階數為 m，則相應的未知參數為未知參數的數目會隨著階數的增加而增加，所需的計算量也會相應增加，而模型的計算精度卻沒有明顯提高。因此，本文在分析PRS模型缺陷的基礎上，提出了改進的多項式響應面（Improved Polynomial Response Surface，簡稱 IPRS）模型。

1．3 改進的多項式響應面法

相對于PRS方法，IPRS方法在近似過程中并沒有忽略殘差項，而是對殘差重新進行近似。在PRS方法中，設生成的近似曲面為f1，近似曲面在采樣點處和已知結果f存在一定的差別，記為殘差R。IPRS方法則是進一步對殘差進行近似估計，并把殘差的估計結果加入到近似曲面中去。在對殘差進行插值估計的過程中，選擇合適的插值估計方法至關重要，在插值估計方法中，徑向基神經網絡是一種性能良好的前向網絡，具有最佳逼近，能克服局部極小值問題的性能［7－8］，用徑向基神經網絡對殘差進行插值處理可以達到較好的插值效果。

基于此，改進的多項式響應面法運用徑向基神經網絡對殘差進行插值處理，具體的步驟如下：

（1）對初始采樣點生成PRS近似曲面f1；

（2）計算近似曲面f1和實際曲面之間的殘差，設實際曲面為Y，則殘差R＝Y－f1；

（3）利用殘差矩陣R進行徑向基神經網絡插值估計，得到插值函數f2；

（4）利用曲面f1和曲面f2分別進行插值計算，得到插值結果y1和y2，并疊加y1和y2，作為最終的插值結果；

計算流程可以參考圖1，這樣，可以得到IPRS方法在n維向量x點處的插值結果fa（x）：

1．4 仿真對比分析

為了驗證IPRS法的擬合效果，利用如下的測試函數對其進行分析：

其中 0．5≤x1，x2≤3．5。

圖1 改進響應面法計算流程圖Fig．1 The calculation flow chart of improved polynomial response surface method

試驗在定義域內隨機對自變量進行取值，為了驗證訓練樣本容量對擬合結果的影響，分別隨機產生30，100，200組訓練樣本，同時在區間內按照線性分布進行插值，插值點數選擇30個，分別采用PRS法和IPRS法對該函數在定義域區間內進行擬合。PRS方法中，與式（8）最接近的模型為二次交互（QI）模型，因此選擇QI模型。在用IPRS法擬合時，需利用徑向基神經網絡對殘差進行處理，本文選擇創建的徑向基網絡，其徑向基神經元數目等于輸入樣本數目，創建該網絡只需要進行一次運算，因而能夠達到較快的訓練速度。然而，創建網絡時需要合理選擇散布常數spread，散布常數控制著網絡輸出的光滑情況，其值越大輸出結果越光滑，但是太大的結果又會造成計算上的困難，經過多次試驗，選擇散布常數spread＝1。用兩種方法擬合后，可以得到各自的近似曲面，圖2～圖4分別是實際圖形以及兩種方法分別在訓練樣本為30、100、200組情況下的插值擬合三維圖。

從三圖中可以明顯看出，在不同訓練樣本容量的情況下，采用IPRS法擬合得到的近似模型相對于采用PRS法擬合得到的近似模型要更加接近于實際模型，隨著采樣樣本數量的增加，PRS法的擬合效果沒有顯著提高，而IPRS法的擬合效果有了明顯的提高。

對于兩種方法的計算結果，可以分別通過以下三種方法來驗證它們的擬合精度：

（1） 預測點的均方根差：MSE＝E［（y＾－y）2］．

（2）為了評價整個模型的精度，可以選擇經驗積累方差標準：其中 m為誤差取樣點總數，在這里m＝30×30；

圖2 訓練樣本為30組情況下的擬合圖形和原始圖形Fig．2 The fitting graphics and original graphics at the training sample group of 30

圖3 訓練樣本為100組情況下的擬合圖形和原始圖形Fig．3 The fitting graphics and original graphics at the training sample group of 100

圖4 訓練樣本為200組情況下的擬合圖形和原始圖形Fig．4 The fitting graphics and original graphics at the training sample group of 200

兩種方法在訓練樣本數目不同的情況下的擬合精度對比如表1所示。從表中可以看出，PRS方法在各驗證方法下的擬合精度都比IPRS方法低，而且隨著訓練樣本數目的增加，其精度提高的也不明顯；相反，IPRS的擬合精度隨著訓練樣本數量的增加有著顯著地提升，這說明IPRS方法在訓練樣本數目較多的情況下比PRS方法的擬合能力更好。因此，將IPRS法代替PRS法而引入VPMCD算法可以達到更理想的模式分類效果。

表1 兩種方法在訓練樣本容量不同的情況下的擬合精度對比Tab．1The contrast of fitting precision of the two methods under different capacity of training sample

1．5 基于改進多項式響應面的VPMCD方法

本文將改進的多項式響應面法代替原始VPMCD中的多項式響應面法，提出了基于改進多項式響應面的VPMCD方法（以下簡稱IPRS－VPMCD），具體的步驟如下所示：

（1）模型訓練過程：

① 對于g類狀態分類問題，共收集n個訓練樣本，每一類狀態樣本數分別為 n1，n2，…，ng。

② 對所有訓練樣本提取特征值X＝［X1，X2，…，Xp］。

③ 對任意被預測變量 Xi，選擇代理模型（L、LI、QI、Q四種PRS模型之一）、預測變量和模型階數。對于不同的特征值，其預測模型類型、預測變量和模型階數都有可能不同。

④ 令k＝1，對于nk個第k類訓練樣本中的任意一個樣本，分別對每一個特征值Xi建立PRS模型，則對每一個特征值可以建立nk個方程；然后用PRS法獲得Xi的插值結果，再以預測誤差平方和為判別函數，得到預測誤差平方和最小的PRS模型；再利用此模型，得到被預測變量的擬合值和實際值之間的殘差Ri，利用徑向基神經網絡對其進行插值估計，得到插值結果和對應的插值模型，并將其與之前的PRS模型相疊加，所得到的模型即為特征值Xi的預測模型

⑤ 令k＝k＋1，循環步驟④，至k＝g結束。

⑥至此，對所有模型類別下的所有特征值都分別建立了預測模型VPMik，其中 k＝1，2，…，g代表不同類別，i＝1，2，…，p代表不同特征值。

（2）模型分類過程：

① 選擇測試樣本，并提取其特征值 X＝［X1，X2，…，Xp］。

② 對于測試樣本的所有特征值 Xi（i＝1，2，…，p），分別采用 VPMki（k＝1，2，…，g）對其進行預測，得到測試值，其中 k＝1，2，…，g代表不同類別，i＝1，2，…，p代表不同特征值。

③ 計算同一類別下所有特征值的預測誤差平方最小為判別函數對測試樣本進行分類，其中k＝1，2，…，g代表不同類別。當在g個預測誤差平方和值中小時，將測試樣本識別為第k類。

2 基于LCD近似熵和IPRS-VPMCD的滾動軸承故障診斷方法

在滾動軸承故障振動信號特征提取算法中，近似熵在描述信號的復雜性時具有較好的抗噪、抗干擾能力，而且利用較短數據即可以較穩健地估計出信號的近似熵。同時，近似熵能用于隨機過程和確定性過程，其取值大小會隨著隨機過程和確定過程的混合比例不同而不同［9］。因此，近似熵能表征信號的復雜程度和產生新模式的概率，可將其應用于故障診斷領域。而且，文獻［10］指出，在實際計算中，數據長度n為有限值，當近似熵的嵌入維數m＝2，相似容量r＝0．1SDx～0．2SDx（SDx為原始數據 x（i）的標準差）時，熵值對 N的依賴程度最小，具有較合理的統計特性。因此，可以通過近似熵算法對滾動軸承振動信號進行特征提取。

然而滾動軸承振動信號往往表現出非平穩性，若直接進行近似熵計算會影響診斷精度，因此，必須先對原始振動信號進行處理。在非平穩信號的分析中，時頻分析方法能同時提供非平穩信號在時域和頻域的局部化信息，因而得到了廣泛的應用。在時頻分析方法中，局部特征尺度分解（Local Characteristic Scale Decomposition，簡稱LCD）算法［10］是一種新的基于極值點的局部特征尺度參數的自適應、非平穩信號處理方法，它能將信號自適應地分解為一系列瞬時頻率具有物理意義的內稟尺度分量（Intrinsic Scale Component，簡稱ISC）。LCD在分解速度、信號的還原度方面都能夠達到較滿意的效果，因而具有非常大的應用前景。本文將近似熵算法和LCD算法相結合應用于滾動軸承故障診斷，通過LCD方法將滾動軸承振動信號分解為若干個平穩的ISC分量，計算每個ISC分量的近似熵，再利用不同ISC分量中提取的近似熵之間存在相互關系這一特點，采用IPRS-VPMCD方法建立預測模型，從而進行模式分類。

基于LCD近似熵和IPRS-VPMCD的滾動軸承故障診斷方法的具體步驟為：

（1）在一定轉速下以采樣率fs對滾動軸承正常、內圈故障、外圈故障和滾動體故障四種狀態進行采樣，每種狀態采集N組樣本。

（2）分別對各類別狀態的原始信號進行LCD分解，每個信號得到若干個ISC分量和余量。選擇合適的i個ISC分量，對所選的ISC分量分別提取近似熵作為特征值，組成特征值向量，每種狀態下得到N×i階的特征值矩陣。

（3）對于不同狀態下的N組特征值，選擇n組樣本作為訓練樣本，采用IPRS-VPMCD方法進行訓練得到預測模型；

（4）將剩余的N－n組樣本作為測試樣本，用訓練得到的預測模型對其進行預測并分類，根據 IPRSVPMCD分類器的輸出結果來確定滾動軸承的工作狀態和故障類型。

3 應 用

本文將IPRS-VPMCD方法應用于滾動軸承故障診斷，采用美國西儲大學電氣工程實驗室的滾動軸承實驗數據來對該方法的有效性和優越性進行驗證，所采用的軸承型號、參數和實驗裝置見文獻［11］。采樣頻率為48 kHz，電機負載為 0．746 kW，轉速為 1 772 r／min，狀態類型分別為：正常狀態、外圈故障、內圈故障、滾動體故障。故障點的直徑為0．177 8 mm，故障深度為0．279 4 mm，每種狀態各得到200個樣本。對各樣本的原始信號進行LCD分解，選擇標準偏差法［12］作為終止判據，選擇鏡像對稱延拓方法［13］減少邊界效應，圖5所示的是某一滾動軸承滾動體故障狀態下的信號及其LCD分解后的分量。

由于滾動軸承的故障信息主要集中在高頻段，因此，可選取前四個ISC分量，并對各分量求取近似熵值（算法中選擇 m＝2，r＝0．2SDx），分別標記為 X1，X2，X3，X4。將所得的近似熵值組成特征向量，以此作為分類器的輸入進行模式識別。各類數據可以提取出200組近似熵，每組近似熵為4個。各類別狀態選擇50組數據作為測試樣本，訓練樣本分別選擇30、90、150組，用IPRS-VPMCD方法進行訓練和測試，得到不同訓練樣本容量下的分類結果。同時，將IPRS-VPMCD方法與VPMCD方法在相同情況下的分類結果進行對比分析，以此驗證兩種模式識別方法在不同訓練樣本容量下的模式分類能力。

圖5 滾動體故障狀態下的振動信號及其LCD分解后的分量Fig．5Vibration signal of rolling element bearing fault state and its components after LCD decomposition

在IPRS-VPMCD方法的殘差估計中，選擇徑向基神經網絡時，網絡輸入中散布常數spread的取值不同會對分類結果有一定的影響，需合理選擇，本實驗參照交叉驗證［14］的方法，并通過多次試驗，選擇 spread＝1．1。另外，為了驗證IPRS-VPMCD方法的適用性及優越性，分別與RBF神經網絡、支持向量機和VPMCD比較，首先經過優化選擇，設置RBF神經網絡訓練誤差的平方和為0．05，設置支持向量機的折衷系數為10，核函數為RBF核函數，然后用四種分類器分別進行訓練分類，最后得到不同訓練樣本數目情況下的分類結果。具體情況如表2所示。

表2 四種分類方法的分類時間和識別率對比Tab．2 The comparison of four kinds of classification methods in the classification time and recognition rate

從表中可以看出，隨著樣本容量的增大，四種方法的分類精度都明顯有所提高，提高了3%左右，但是無論樣本容量多少，RBF神經網絡和VPMCD的分類精度都明顯低于支持向量機和IPRS-VPMCD。另外，再從分類時間上看，VPMCD和IPRS-VPMCD的分類時間明顯快于另外兩種，尤其是支持向量機，綜合分析，IPRSVPMCD方法在訓練過程中除了獲得最佳PRS模型之外，對其殘差還進行了插值處理，保留了部分重要信息，因此，IPRS-VPMCD方法在保證分類效率的同時，有效地提高了計算精度，并且分類精度能夠隨訓練樣本容量增加而明顯提高。

4 結 論

改進的多項式響應面法將多項式相應面法中忽略掉的殘差項重新進行了插值近似，保留了模型擬合過程中的重要信息，因而能夠獲得更加合理的近似模型。通過仿真分析，將多項式響應面法和改進的多項式響應面法進行擬合效果的比較，結果證明了改進多項式響應面法具有更好的擬合效果，而且能夠隨著訓練樣本容量的增加較明顯地提高擬合精度。

將改進的多項式響應面法代替原始VPMCD方法中的多項式響應面法，提出了IPRS-VPMCD方法，并將其應用于滾動軸承故障診斷中。通過實驗，將IPRSVPMCD方法和VPMCD方法在訓練樣本容量不同情況下的分類結果進行對比分析，結果證明了它的優越性。

［1］Wang C C，Kang Y，Shen P C，et al．Applications of fault diagnosis in rotating machinery by using time series analysis with neural network［J］．Expert System with Application，2010，37（2）：1696－1702．

［2］Fei Sheng-wei，Zhang Xiao-bin．Fault diagnosis of power transformer based on support vector machine with genetic algorithm［J］．Expert Systems with Application，2009，36（8）：11352－11357．

［3］Raghuraj R， Lakshminarayanan S． Variable predictive models-A new multivariate classification approach for pattern recognition applications［J］．Pattern recognition，2009，42（1）：7－16．

［4］穆雪峰，姚衛星，余雄慶，等．多學科設計優化中常用代理模型的研究［J］．計算力學學報，2005，36（22）：5－8．MU Xue-feng，YAO Wei-xing，YU Xiong-qing，et al．A survey of surrogate models used in MDO ［J］．Chinese Journal of Computational Mechanics，2005，36（22）：5－8．

［5］Joyce A P，Leung SS．Use of response surface methods and path of steepest ascent to optimize ligand-binding assay sensitivity［J］．Journal of Immunological Methods，2013，392（1－2）：12－23．

［6］潘雷，谷良賢，閻代維．改進響應面法及其近似性能研究［J］．宇航學報，2009，30（2）：806－810．PAN Lei， GU Liang-xian， YAN Dai-wei． Research on improved response surface method and its performance［J］．Journal of Astronautics，2009，30（2）：806－810．

［7］彭磊，劉莉，龍騰．基于動態徑向基函數代理模型的優化策略［J］．機械工程學報，2011，47（7）：164－170．PENG Lei，LIU Li，LONG Teng．Optimization strategy using dynamic radial basis function metamodel［J］．Journal of Mechanical Engineering，2011，47（7）：164－170．

［8］吳月明，王益群，李莉．BP神經網絡與廣義RBF神經網絡在產品壽命分布模型識別中的應用研究［J］．中國機械工程，2006，17（20）：2140－2143．WU Yue-ming，WANG Yi-qun，LI Li．Application study on BP network and generalized RBF network in estimating distribution model of mechanical products［J］． China mechanical engineering，2006，17（20）：2140－2143．

［9］Lee JH，Kim J，Kim H J．Development of enhanced wignerville distribution function［J］．Mechanical Systems and Signal Processing，2001，13（2）：367－398．

［10］Classen T，Mecklenbrauker W．The aliasing problem in discrete-time wigner distribution［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1983，31（5）：1067－1072．

［11］Case Western Reserve University Bearing Data Center．Bearing Data Center Fault Test Data．［EB／OL］．［2009－10－01］． http： ／／www．eecs． case． edu／laboratory／bearing．

［12］Huang N E，Shen Z，Long S R．The empirical mode decomposition and the Hilbert spectrum for nonlinear nonstationary time Series Analysis［J］． Proc． Roy． Soc．London，1998，454：903－995．

［13］Rilling G，Flandrin P，Goncalves P．On empirical mode decomposition and its algorithms［J］．In IEEE-EURASIP Workshop on Nonlinear Signaland Image Processing，Grado，Italy，2003：8－11．

［14］Breiman L．Heuristics of instability and stabilization in model selection［J］．The Annals of Statistics，1996，24（4）：2350－2383．