風(fēng)力機主軸承故障監(jiān)測方法?

鄭玉巧，魏劍峰，朱 凱，董 博

（蘭州理工大學(xué)機電工程學(xué)院 蘭州，730050）

引言

風(fēng)力機主軸承作為風(fēng)力機傳動系統(tǒng)的核心部件，承受輪轂、葉片等部件帶來的多種變載沖擊作用，極易發(fā)生故障［1-2］。由于風(fēng)力機主軸承在距離地面幾十米甚至上百米高的機艙內(nèi)，一旦發(fā)生故障維護(hù)難度大、成本高，若未及時發(fā)現(xiàn)并解決故障會給風(fēng)力機運行帶來嚴(yán)重影響甚至造成安全事故。因此，實現(xiàn)風(fēng)力機主軸承的故障監(jiān)測，及時發(fā)現(xiàn)故障征兆并開展維護(hù)，對于風(fēng)力機健康運行、降低運維成本具有重要意義。

傳統(tǒng)風(fēng)力機部件故障監(jiān)測方法有振動分析［3-5］、聲發(fā)射監(jiān)測［6-7］和電氣信號監(jiān)測［8-9］等。其中振動監(jiān)測法和聲發(fā)射監(jiān)測法需要安裝大量的傳感器來采集信號，成本昂貴；此外安裝的傳感器一旦發(fā)生故障，不僅導(dǎo)致采集的數(shù)據(jù)可靠性降低，還增加額外的運維成本。相比風(fēng)力機部件的振動信號和聲信號，電氣信號的采集不需要額外安裝傳感器［10］，但電氣信號中所包含的故障信息往往比較微弱，信噪比低，監(jiān)測結(jié)果誤差較大。因此，上述方法難以在風(fēng)力機部件實際監(jiān)測中大范圍推廣。

目前，多數(shù)風(fēng)電場通過安裝數(shù)據(jù)采集與監(jiān)控系統(tǒng)（supervisory control and data acquisition，簡 稱SCADA）采集儲存風(fēng)力機各部件運行參數(shù)及環(huán)境參數(shù)，因此選擇SCADA系統(tǒng)所記錄的風(fēng)力機主軸承運行參數(shù)開展其故障監(jiān)測研究。主軸承故障表現(xiàn)形式為磨損、點蝕或塑性變形等，故障早期主軸承滾過表面損傷點產(chǎn)生的異常摩擦又加劇故障程度，引起主軸承溫度特征的異常波動。鑒于此，本研究以SCADA系統(tǒng)所采集的風(fēng)力機主軸承運行參數(shù)為研究對象，分別采用多元線性回歸（multiple linear regression，簡稱MLR）、灰色模型（grey model，簡稱GM）及支持向量機回歸（support vector machine regression，簡稱SVR）建立主軸承溫度單一預(yù)測模型及它們的組合預(yù)測模型，引入滑動窗口方法實時監(jiān)測溫度預(yù)測殘差的均值和標(biāo)準(zhǔn)差變化情況，通過溫度預(yù)測殘差均值（或標(biāo)準(zhǔn)差）的置信區(qū)間與設(shè)定臨界值的對比來判斷主軸承運行狀態(tài)，從而實現(xiàn)對風(fēng)力機主軸承的故障監(jiān)測。

1 主軸承溫度組合預(yù)測模型的構(gòu)建

1.1 問題描述與數(shù)據(jù)預(yù)處理

風(fēng)力機SCADA系統(tǒng)通常以10 min為采樣周期記錄并傳輸風(fēng)速、輸出功率及主軸承溫度等運行參數(shù)［11］，由于風(fēng)力機SCADA系統(tǒng)輸出的各項運行參數(shù)單位不同，若直接用于建模分析對預(yù)測結(jié)果造成較大誤差，因此需要對原始數(shù)據(jù)進(jìn)行歸一化處理，消除數(shù)據(jù)間的量綱影響。公式如下

其中：xmax為監(jiān)測數(shù)據(jù)中的最大值；xmin為監(jiān)測數(shù)據(jù)中的最小值；xi為原始監(jiān)測值；x*為歸一化處理后的數(shù)據(jù)，x*的取值范圍為［0，1］。

1.2 多元線性回歸（MLR）模型

MLR模型闡述一個因變量依賴多個自變量的變化規(guī)律［12］。風(fēng)力機主軸承溫度的變化視為多個變量參數(shù)影響下的共同作用結(jié)果，因此可將主軸承溫度作為因變量（即建模輸出變量），影響主軸承溫度變化的運行參數(shù)作為自變量（即建模輸入變量），建立主軸承溫度的多元線性回歸預(yù)測模型步驟如下。

設(shè)因變量（即風(fēng)力機主軸承溫度）為y，m個自變量（即影響主軸承溫度變化的變量）分別為x1，x2，…，xm的多元線性回歸方程可表示為

其中：θ0為回歸常數(shù)；θ1，θ2，…，θm為回歸系數(shù)。

假設(shè)在實際問題中已測得n組監(jiān)測數(shù)據(jù)（yi；xi1，，xi2，…，xim）（i=1，2，…，n），則式（2）可表示為

采用最小二乘法［12］可求解得到式（3）的回歸參數(shù)，將所求回歸參數(shù)代入式（2）后，需要對回歸方程進(jìn)行統(tǒng)計量檢驗（F檢驗、T檢驗及多重共線性檢驗），若不滿足檢驗要求則需要重新尋找影響主軸承溫度變化的變量以便重新建模，滿足要求則可直接輸出模型。建模流程如圖1所示。

1.3 灰色預(yù)測模型GM（1，N）

灰色模型（GM）［13］具有建模數(shù)據(jù)少、預(yù)測效果好等優(yōu)點。傳統(tǒng)的GM（1，1）模型是針對單一變量的預(yù)測，然而影響風(fēng)力機主軸承溫度的因素是多變量共同作用的結(jié)果，因此GM（1，N）模型更適合建立風(fēng)力機主軸承溫度預(yù)測模型。

由此可建立GM（1，N）模型為

其中：a為系統(tǒng)發(fā)展系數(shù)，反映預(yù)測值的發(fā)展趨勢；bj為驅(qū)動系數(shù)，反映原始數(shù)據(jù)的內(nèi)在變化。

設(shè)置參數(shù)列μ=（a，b2，b3，…，bm）T，根據(jù)最小二乘法可求解式（6）得到參數(shù)列μ的值。

對GM（1，N）模型（式6）進(jìn)行白化可得

求解式（7）得到時間響應(yīng)函數(shù)，表示為

最終可得累減還原式（即灰色預(yù)測值）為

1.4 支持向量機回歸（SVR）模型

支持向量機回歸的核心思想是建立一個最優(yōu)分類面，使得所有訓(xùn)練樣本離最優(yōu)分類面的誤差最?。?4］。風(fēng)力機主軸承溫度的SVR模型的訓(xùn)練與預(yù)測借助LIBSVM工具箱實現(xiàn)。首先，選定一組風(fēng)力機主軸承SCADA數(shù)據(jù)作為原始數(shù)據(jù)，對其進(jìn)行歸一化處理；然后，根據(jù)工具箱提供的交互檢驗（cross validation，簡稱VC）功能尋求最優(yōu)懲罰參數(shù)c和核函數(shù)參數(shù)g（由函數(shù)SVMcgForRegress實現(xiàn)）；最后，根據(jù)得到的最優(yōu)參數(shù)c，g對SVR模型進(jìn)行預(yù)測。

圖2 SVR算法流程圖Fig.2 Flow chart of SVR algorithm

1.5 組合預(yù)測模型

假設(shè)某預(yù)測問題共有f種單一預(yù)測模型，令yt表示t時刻的實際監(jiān)測值，y?it表示第i種預(yù)測模型在t時刻的預(yù)測值（其中：i=1，2，…，f；t=1，2，…，n），第i種預(yù)測模型在t時刻的相對誤差為eit

由于每種單一預(yù)測模型具有其不足之處，而組合預(yù)測模型將不同類型的單一模型按照一定的權(quán)重比例重新組合可發(fā)揮各種單一模型優(yōu)勢，有效削弱各種單一模型的不足以達(dá)到提高預(yù)測精度的目的。因此，建立基于熵值法［15］的組合模型，步驟如下：

1）求解第i種單一預(yù)測模型在t時刻的相對誤差比重pit

2）求解第i種單一預(yù)測模型與實際值的相對誤差熵值gi

3）求解相對誤差的變異系數(shù)di

4）求解各種單一預(yù)測模型的權(quán)重wi

5）求解組合預(yù)測模型值

1.6 各模型綜合評價指標(biāo)

為實現(xiàn)對各個預(yù)測模型的定量評價，選取均方根誤差（root mean square error，簡稱RMSE）、平均絕對百分誤差（mean absolute percentage error，簡稱MAPE）、平均絕對誤差（mean absolute error，簡稱MAE）及決定系數(shù)（R2）作為評價指標(biāo)，求解公式如下

其 中：yt為 實 際 值為 預(yù) 測 值；RMSE，MAPE及MAE用來檢驗實際值與預(yù)測值的偏差及波動性大小，數(shù)值越小則表明預(yù)測模型精度越高；R2為自變量所解釋的變化量占總變化量的比例，R2取值為［0，1］，R2越接近于1則表明預(yù)測精度越高。

2 殘差統(tǒng)計特性分析及故障監(jiān)測

2.1 滑動窗口殘差統(tǒng)計特性分析

假設(shè)風(fēng)力機主軸承溫度在t時刻的實際監(jiān)測值為yt，預(yù)測模型在t時刻的預(yù)測值為則預(yù)測殘差序 列 為預(yù) 測 殘 差序列的統(tǒng)計特性可由殘差均值及殘差標(biāo)準(zhǔn)差反映。為準(zhǔn)確描述預(yù)測殘差的統(tǒng)計特性（均值和標(biāo)準(zhǔn)差）變化情況，采用滑動窗口方法［16］求解預(yù)測殘差序列的窗口均值及窗口標(biāo)準(zhǔn)差，滑動窗口如圖3所示。

圖3 滑動窗口示意圖Fig.3 Schematic diagram of sliding window

設(shè)定寬度為N的滑動窗口，則窗口內(nèi)的N個殘差的窗口均值μ和窗口標(biāo)準(zhǔn)差σ求解公式為

2.2 主軸承的故障監(jiān)測

主軸承正常工作時，溫度預(yù)測值與實際溫度值接近，其殘差序列分布穩(wěn)定；然而主軸承發(fā)生故障會引起其溫度的異常波動，進(jìn)而帶來預(yù)測殘差統(tǒng)計特性的異常變化。具體表現(xiàn)為以下3種情況：①主軸承溫度預(yù)測殘差的窗口均值μ未發(fā)生較大改變，但其窗口標(biāo)準(zhǔn)差σ變化顯著，殘差范圍增大；②主軸承溫度預(yù)測殘差的窗口標(biāo)準(zhǔn)差σ未發(fā)生較大改變，但其窗口均值μ變化顯著，出現(xiàn)系統(tǒng)偏差；③主軸承溫度預(yù)測殘差的窗口均值μ及標(biāo)準(zhǔn)差σ均變化顯著。

基于上述情況，可根據(jù)預(yù)測殘差序列的統(tǒng)計特性變化來監(jiān)測主軸承運行狀態(tài)，設(shè)定殘差均值和殘差標(biāo)準(zhǔn)差故障臨界值分別為μv和σv。當(dāng)溫度預(yù)測殘差的均值或標(biāo)準(zhǔn)差超過某一臨界值時，系統(tǒng)產(chǎn)生故障報警告知維護(hù)人員進(jìn)行檢修。假設(shè)滑動窗口的溫度殘差序列均值絕對值最大值為μmax，標(biāo)準(zhǔn)差絕對值最大值為σmax，則主軸承溫度故障的臨界標(biāo)準(zhǔn)為

式（22）中k1，k2可由現(xiàn)場運維人員根據(jù)經(jīng)驗確定。根據(jù)預(yù)測模型對主軸承溫度進(jìn)行預(yù)測時，會存在各種不確定性因素對輸出結(jié)果造成干擾［17］。為簡化起見，將預(yù)測殘差序列假設(shè)為一組均值和方差未知的正態(tài)分布數(shù)據(jù)，在求解殘差窗口均值和標(biāo)準(zhǔn)差時需要給出置信度為1-α的均值和標(biāo)準(zhǔn)差的置信區(qū)間，表示為

其中：tα/2為t分布的α/2分位點為χ2分布的α/2分位點，其值在給定置信度1-α后可通過查閱t分布臨界值表和χ2分布臨界值表得知。

當(dāng)殘差均值或標(biāo)準(zhǔn)差的置信區(qū)間超過各自臨界值時，發(fā)出風(fēng)力機主軸承故障報警信息。

3 案例分析

3.1 建模變量的選取

根據(jù)文獻(xiàn)［18］可知，關(guān)鍵部件溫度相關(guān)聯(lián)的特征量主要包括環(huán)境溫度、風(fēng)速、轉(zhuǎn)速及功率等，此外由于部件溫度的熱慣性特點使得前一時刻的部件溫度對當(dāng)前時刻部件溫度有直接影響。因此，建立主軸承溫度預(yù)測模型時初步選擇風(fēng)力機輸出功率、風(fēng)速、環(huán)境溫度、主軸轉(zhuǎn)速及上時刻主軸承溫度作為輸入變量，主軸承當(dāng)前時刻溫度（y）作為輸出變量。

選取某風(fēng)力機2019年4月主軸承正常狀態(tài)時的運行參數(shù)進(jìn)行建模，根據(jù)所建模型對5月份主軸承溫度進(jìn)行預(yù)測。按照1.1部分所述步驟對原始SCADA數(shù)據(jù)進(jìn)行歸一化處理，歸一化后的部分?jǐn)?shù)據(jù)如表2所示。

表1 SCADA數(shù)據(jù)的歸一化Tab.1 Normalization of SCADA data

3.2 MLR，GM（1，N），SVR及組合預(yù)測模型的求解

3.2.1 多元線性回歸模型(MLR)

以風(fēng)力機輸出功率（x1）、風(fēng)速（x2）、空氣溫度（x3）、主軸轉(zhuǎn)速（x4）及上時刻主軸承溫度（x5）為自變量，風(fēng)力機主軸承當(dāng)前時刻溫度（y）為因變量建立多元線性回歸方程，根據(jù)1.1部分所述步驟求解得到的MLR模型為

對式（25）進(jìn)行統(tǒng)計量檢驗時發(fā)現(xiàn)輸出功率（x1）所對應(yīng)的方差擴大因子（variance inflation factor，簡稱VIF）為15.641 8，風(fēng) 速（x2）所 對 應(yīng) 的VIF為15.713 4，表明所建模型存在嚴(yán)重的多重共線性（VIF＞10），可能會給回歸結(jié)果帶來較大誤差，因此需要重新修正模型輸入變量，以消除多重共線性。

消除多重共線性的常用方法則是剔除VIF最大值所對應(yīng)的變量，即風(fēng)速（x2）；剔除后的MLR修正模型自變量為風(fēng)力機輸出功率（x1）、空氣溫度（x3）、主軸轉(zhuǎn)速（x4）及上時刻主軸承溫度（x5），根據(jù)1.1部分所述步驟求解得到的MLR修正模型為

據(jù)檢驗，式（26）已消除多重共線性，F(xiàn)檢驗及T檢驗各項指標(biāo)均滿足顯著性條件。利用修正后的MLR模型對5月份主軸承溫度進(jìn)行預(yù)測，部分預(yù)測值及殘差值見表2。

表2 MLR模型部分溫度預(yù)測結(jié)果值Tab.2 Temperature prediction results of MLR model℃

3.2.2 灰色預(yù)測模型GM(1,N)的求解

基于MLR修正模型所得自變量及因變量，選擇當(dāng)前時刻主軸承溫度、功率、空氣溫度、主軸轉(zhuǎn)速及上時刻主軸承溫度的時間序列為原始數(shù)據(jù)序列其中：當(dāng)前時刻主軸承溫度數(shù)據(jù)序列為系統(tǒng)特征數(shù)據(jù)序列；而功率數(shù)據(jù)序列空氣溫度數(shù)據(jù)序列主軸轉(zhuǎn)速數(shù)據(jù)序列以及上時刻主軸承溫度序列均為相關(guān)因素序列。選取2019年4月份SCADA數(shù)據(jù)依據(jù)1.2所述步驟進(jìn)行建模，求解得到GM（1，N）模型參數(shù)列μ的值為

μ=(2.014 5，-0.000 6，0.062 5，0.146 5，1.972 3)T

將參數(shù)列μ代入式（8），得到時間相應(yīng)函數(shù)，而后根據(jù)式（9）累減可得到GM（1，N）預(yù)測值，部分預(yù)測值及殘差見表3。

表3 GM(1,N)模型部分溫度預(yù)測結(jié)果值Tab.3 Temperature prediction results of GM(1,N)model ℃

3.2.3 支持向量機回歸模型(SVR)的求解

基于MLR修正模型所得自變量及因變量，選擇2019年4月份的當(dāng)前時刻主軸承溫度、功率、空氣溫度、主軸轉(zhuǎn)速及上時刻主軸承溫度等序列數(shù)據(jù)用來訓(xùn)練SVR模型，用訓(xùn)練好的模型來預(yù)測5月份主軸承溫度值。

給定參數(shù)c和g的初始范圍均為［-8，8］，根據(jù)交互檢驗法求解得到SVR模型最優(yōu)參數(shù)c和g為：Bestc=1.414 2，Bestg=0.353 5；根據(jù)所求最優(yōu)參數(shù)c和g對SVR模型進(jìn)行訓(xùn)練，然后對5月份數(shù)據(jù)進(jìn)行預(yù)測，部分預(yù)測值及殘差值見表4。

3.2.4 MLR、GM(1,N)及SVR組合模型的求解

假設(shè)MLR修正模型所得主軸承溫度預(yù)測值為y?1t，GM（1，N）模型所得預(yù)測值為y?2t，SVR模型所得預(yù)測值為則應(yīng)用熵值法確定的各模型最優(yōu)權(quán)重為w1=0.125 8，w2=0.310 5，w3=0.563 7，因此可得組合模型表達(dá)式為根據(jù)所求組合模型對5月份數(shù)據(jù)進(jìn)行預(yù)測，部分預(yù)測值及殘差見表5。

表4 SVR模型部分溫度預(yù)測結(jié)果值Tab.4 Temperature prediction results of SVR model℃

表5 組合模型部分溫度預(yù)測結(jié)果值Tab.5 Temperature prediction results of Combined model ℃

3.3 各模型預(yù)測結(jié)果對比分析

根據(jù)各模型對5月份主軸承溫度數(shù)據(jù)進(jìn)行預(yù)測，實際值與各模型預(yù)測值變化趨勢如圖4所示。

圖4 風(fēng)力機主軸承溫度預(yù)測圖Fig.4 Forecast figure of wind turbine main bearing temperature

圖4 顯示在第48，1 220，1 928及2 206點附近處主軸承溫度較低，查詢主狀態(tài)日志可知是由于風(fēng)力機停機重新啟動所造成。另外從圖中可看到，不同模型所求主軸承溫度預(yù)測值與實際值變化趨勢相近，擬合程度較高；但根據(jù)圖4不易直觀判斷哪種模型預(yù)測值與實際值吻合度更好。為進(jìn)一步對比各模型優(yōu)劣性，進(jìn)而求解各模型綜合評價指標(biāo)值，見表6。

表6 各模型指標(biāo)對比值Tab.6 Comparison results of each model index

由表5的各指標(biāo)數(shù)據(jù)對比可知，組合模型的RMSE，MPE及MAE最小，R2最大，表明組合模型偏差及波動性更小，在各模型中具有最高的預(yù)測精度；其中組合模型的R2值分別較其他模型提高0.049 3，0.002 7和0.000 2。因此，選擇在組合預(yù)測模型基礎(chǔ)上引入滑動窗口方法實現(xiàn)對風(fēng)力機主軸承的故障監(jiān)測。此外在單一預(yù)測模型中，SVR模型具有良好的預(yù)測精度，各項指標(biāo)均和組合模型指標(biāo)相近。

3.4 基于滑動窗口的主軸承溫度故障狀態(tài)監(jiān)測

查詢風(fēng)力機主狀態(tài)日志可知，5月份風(fēng)力機主軸承未發(fā)生故障。為驗證方法有效性，人為模擬主軸承因發(fā)生故障而導(dǎo)致主軸承溫度升高的情況。選擇5月份主軸承溫度較高的12～15號3天累計568個溫度點進(jìn)行研究，取滑動窗口寬度N=20，置信水平α=0.05，根據(jù)式（23），（24）求解得到主軸承正常狀態(tài)下預(yù)測殘差均值及標(biāo)準(zhǔn)差的各自95%置信區(qū)間，它們的變化趨勢如圖5所示。

圖5 殘差窗口均值和標(biāo)準(zhǔn)差Fig.5 Residual window mean and standard deviation

圖5 （a）為主軸承溫度殘差窗口均值變化趨勢圖，從圖5（a）可知主軸承溫度殘差窗口均值最大值為μmax=0.223 5（即A點）；圖5（b）為主軸承溫度殘差窗口標(biāo)準(zhǔn)差變化趨勢圖，由圖5（b）可知主軸承溫度殘差窗口標(biāo)準(zhǔn)差最大值為σmax=0.764 6（即B點），由于是人為模擬故障，所以殘差標(biāo)準(zhǔn)差并未發(fā)生變化。故選擇殘差窗口均值展開分析，取k1=2，則窗口均值臨界值根據(jù)式（22）可得μv=±0.447 0。

從第303個溫度點開始人為加入步長為0.015℃的溫度累積偏移量，基于組合模型的主軸承溫度預(yù)測殘差均值變化趨勢如圖6所示。

圖6 模擬故障狀態(tài)的殘差窗口均值Fig.6 Residual window mean value of simulated fault state

圖6 顯示，主軸承正常運行時殘差均值95%置信區(qū)間均未超過臨界值，在C點即第284個窗口（303-20+1）加入溫度偏移后，滑動窗口殘差均值逐漸增大，在D點殘差窗口均值95%置信上界超過臨界值0.447 0，引發(fā)主軸承故障報警。綜上可得，當(dāng)主軸承發(fā)生故障導(dǎo)致其溫度異常時，基于滑動窗口方法的故障監(jiān)測方法可及時發(fā)現(xiàn)故障并報警，有效監(jiān)控主軸承異常狀態(tài)。

4 結(jié)束語

基于風(fēng)力機SCADA數(shù)據(jù)分別建立主軸承溫度的MLR，GM（1，N），SVR 3種單一預(yù)測模型及它們的組合預(yù)測模型；其中組合模型有著更優(yōu)越的預(yù)測性能，其決定系數(shù)分別比其他模型提高0.049 3，0.002 7及0.000 2。此外SVR模型在單一預(yù)測模型中有著更高的預(yù)測精度。在組合預(yù)測模型基礎(chǔ)上引入滑動窗口方法求解分析溫度預(yù)測殘差的統(tǒng)計特性（均值和標(biāo)準(zhǔn)差），當(dāng)窗口均值或標(biāo)準(zhǔn)差的置信區(qū)間超過設(shè)定臨界值時可判斷出主軸承發(fā)生故障，有效監(jiān)測主軸承運行狀態(tài)。筆者所提方法可用于工程實際中風(fēng)力機或其他機械設(shè)備的關(guān)鍵部件溫度預(yù)測及故障監(jiān)測。